The impact on distant fly-bys on the rate of binary primordial black hole mergers

By performing Monte Carlo simulations of the evolution of binary primordial black hole (PBH) systems, we estimate the effect of distant encounters with single PBHs upon the coalescence time and merger rate of binary PBHs. We find that, for models where PBHs compose a large fraction of dark matter, $\fPBH\sim 1$, the expected fractional change in coalescence time is negligible, of order $10^{-6}$ for most binaries. For models with significantly lower PBH abundances, $\fPBH\ll 1$, we find that the average change in binary lifetime due to encounters can be as large as $\mathcal{O}(10^{-2})$, with a small number of binaries experiencing an order unity change in lifetime. In the absence of encounters, we also compare the use of an analytic approximation for the coalescence time to numerically evolving the binary system, finding that the analytic approximation results in an order $10\%$ error in the coalescence time. However, when these effects are taken into consideration, there is a negligible change to the calculated merger rate, placing previous constraints on the PBH abundance arising from observed gravitational wave signals from merging binary black holes on a more secure footing.


I. INTRODUCTION
Since the first detection of gravitational waves by the LIGO scientific collaboration on September 14, 2015 from the merger of 2 black holes (BHs) [1], as well as subsequent detections [2][3][4][5][6][7][8] there has been speculation on the origin of these BHs. Did they form through an astrophysical channel, or did the BHs have a primordial origin?
Primordial black holes (PBHs) could form through a variety of mechanisms in the early universe, including the collapse of large density perturbations [62,63], from cosmic strings [64], or from bubble collisions [65]. Part of the motivation for considering PBHs as a candidate for the LIGO BHs is the observed spins of the merging BHs -which can be difficult to explain from astrophysical BHs, but is a natural prediction for PBHs [66][67][68].
The question of whether LIGO had detected PBHs was quickly investigated by Bird et al. [69] and Clesse and García-Bellido [70], and the initial findings were that the observed merger rate matched closely with the merger rate predicted if dark matter (DM) was composed entirely of PBHs. Since then, however, the calculation has been refined, notably including the formation of binary systems in the early universe [71][72][73][74][75][76][77], and the current consensus is that the observed merger rate is too low for PBHs to make up the entirety of DM, implying that, at most, PBHs could compose O(0.1%) of DM, f PBH 0.001.
However, such calculations of the merger rate today typically ignore the effect of other nearby objects on the evolution of binary primordial black holes (BPBHs), although attempts have been made to account for this. Vaskonen and Veermäe [78] considered the disruption of binaries located in haloes undergoing core collapse, which slightly weakened constraints on the PBH abundance coming from the observed merger rate, whilst Raidal et al. [72] used an N -body approach to study binary PBHs in the early universe, finding that initial binaries are likely to be disrupted if the abundance of PBHs is large, f PBH 0.1.
In this paper, we will study the effect of "fly-bys" (other PBHs passing near the binary system) in the late universe. A simple calculation was performed in [79], which found that the chance of single PBHs passing by closely enough to have a significant effect on the coalescence time of BPBH was unlikely, although the approximate calculation employed therein may not be accurate given the highly eccentric nature of BPBHs [58]. Here, we will investigate the cumulative effect of many such fly-bys using more accurate Monte Carlo methods that employ analytic equations to predict the effects of each fly-by.
In particular, we will use the secular approximation, in which the BPBH orbital period is much shorter than the passage time-scale of the perturbing body (e.g., [58,[80][81][82][83][84]). This approximation works well for the overwhelming majority of perturbers in our scenario, as we will show in section III. We combine Monte Carlo sampling of perturbers with the decay of the BPBH orbit due to the emission of GWs, similar in approach to [59] who considered BH binaries of astrophysical origin in globular clusters.
The organisation of this paper is as follows: • Section II will discuss the formation of binaries and initial conditions of BPBHs in the early universe; • Section III discusses the secular regime relevant for distant fly-bys; • Section IV presents an analytic estimate of the change in coalescence time due to the effect of fly-bys; • Section V deals with how the evolution of binary systems is calculated using a Monte Carlo approach; • Section VI presents the results of the Monte Carlo simulations, and finally; • Section VII discusses the conclusions and implications from our investigation.

II. BINARY FORMATION AND INITIAL CONDITIONS
We will be performing a Monte Carlo procedure to model the evolution of binary systems as other PBHs pass nearby, and to this end, we require that the initial distribution of binary systems and their orbital parameters closely match those expected in the early universe. Therefore, in this section we will discuss the formation of binary PBH systems in the early universe, and the calculation of their initial conditions. In order to do this with, we will follow the derivation of initial conditions given in Raidal et al. [72].

A. Primordial black hole mass
We will consider the case that PBHs form from the collapse of large-amplitude overdensities in the early universe. Assuming PBHs form with a relatively narrow mass function (for example, arising from a narrow peak in the primordial power spectrum leading to enhanced PBH formation at those scales), the mass function can be well approximated with a log-normal distribution [99,104] where m is the PBH mass, m c is the mass at which the distribution peaks, and σ m is the width of the distribution. In this paper, we will treat m c and σ m as free parameters, although a best-fit model to the black hole coalescence events observed by LIGO suggests m c ≈ 20 M and σ m ≈ 0.5 [72].
Note that this definition of the mass function is normalised to integrate to unity, dm ψ(m) = 1 (and can therefore be interpreted as the probability distribution function (PDF) for PBH mass).
To describe the total abundance of PBHs, we will use the parameter f PBH , which is the fraction of DM composed of PBHs, where Ω is the density parameter for PBHs or cold dark matter (CDM). The number density n of PBHs can therefore be expressed as wherem = m c exp σ 2 m /2 is the mean PBH mass. We note that, since the PBH density evolves as matter, the (mean) number density expressed in comoving coordinates is constant with respect to time.

B. Initial semi-major axis
In the absence of primordial non-Gaussianity, PBHs are expected to follow a Poissonian spatial distribution. The PDF for the radial distance r 0 from a given PBH to its nearest neighbor at formation is then given by where n is the number density of PBHs at the time of formation. In order for a pair of PBHs to form a binary system, we require that the PBHs are allowed to decouple from the Hubble flow and fall towards each other without being disrupted by other nearby perturbations/PBHs. To achieve this, we include an exclusion zone around the pair, in which no other PBHs are found -which is a factor A larger than the initial separation of the PBH pair. The probability of finding the nearest neighbor in the range r → r + dr, and no other PBHs within the exclusion zone (a sphere of radius Ar) is then We will later want the distribution for the initial separation of binary systems (rather than PBH pairs which may not form binaries), and so the distribution is correctly normalised as where the factor A 3 can be considered a normalisation factor to ensure that dr 0 P bin (r 0 ) = 1, since we want the distribution of r 0 for PBH pairs which form binaries.
The number density of initial pairs forming binary systems, with initial separation in the range r 0 → r 0 + dr 0 and masses in the range m 1,2 → m 1,2 + dm, can be calculated by multiplying the relevant probabilities: where a factor of 1/2 is included to avoid overcounting, and dn(m) = n. This equation is equivalent to equation (2.18) in Raidal et al. [72], with 4πr 2 0 dr 0 = dV (x 0 ) and 4 3 π(Ar 0 ) 3 n =Ñ (y). When PBHs form significantly close to each other relative to the average, the pair can be considered a matter overdensity, which eventually decouples from the Hubble flow and "collapses" to form a binary system. We follow Raidal et al. [72] and define the quantity, where ρ M is the matter energy density. During radiation domination, such matter perturbations collapse when ρ r a −4 = δ b ρ M a −3 , where ρ r is the radiation energy density. An approximate estimate for the decoupling scale factor is then given by where the eq represents matter-radiation equality. The approximate value for the initial semi-major axis r a of the binary when it decouples from the Hubble flow is given by Raidal et al. [72] r a ≈ 0.1a dc r 0 = 0.1a eq r 0 δ b .
An estimate for the initial semi-major axis can therefore be calculated from the PBHs masses and initial separation, for which the PDFs are known. As will be seen, our conclusion is not sensitive to small errors in the initial conditions, and so the approximations made here are considered acceptable.

C. Initial angular momentum
After a PBH pair decouples from the Hubble flow, in the absence of external perturbations/objects modifying the gravitational field around the PBH pair, the pair would fall straight back towards each other and immediately coalesce. However, nearby density perturbations and PBHs can provide a torque to the system, imparting sufficient angular momentum to the pair to prevent a head-on collision and instead form a stable binary system.
When describing the orbits of binary systems, it will be helpful for use the eccentricity e to describe the ellipticity, defined as where r a is the apoapsis, and r p is the periapsis. However, when calculating the distribution of initial conditions, it will be more helpful to describe the dimensionless angular momentum j, related to the eccentricity as Raidal et al. [72] gives an order of magnitude estimate for the initial angular momentum under the assumption that most of the torque is generated by the nearest PBH to the pair. Depending on the exact configuration of nearby density perturbations and PBHs, the actual angular momentum j will vary relative to this value, and the PDF itself varies depending on the expected number of PBHs in the exclusion zone,Ñ . Raidal et al. [72] provides several forms for the PDF depending, depending on the value ofÑ (y). In the limitÑ → 0, the PDF of the initial angular momentum j is given by a power-law: In the limitÑ → ∞, a Gaussian-like PDF is found instead: with σ j given by Here, σ 2 M ≡ Ω M Ω DM δ 2 M is the re-scaled variance of the matter density perturbation, Ω M and Ω DM are the density parameter for matter and DM respectively. We will follow Ali-Haïmoud et al. [76] and take δ 2 M = 0.005 2 when a numerical estimate is required. If the variance of matter perturbations is dominated by the Poisson noise generated by PBHs, the power-law PDF is expected to hold -that is, if f PBH σ M . As can be seen in figures 1 and 2, for PBHs expected to merge around today, that whilst the tails of the distributions can be quite different, the peaks of the distributions for j are quite similar, and so we will proceed by using the power-law distribution to generate initial conditions -which has a negligible impact on our results.
The required orbital parameters can then be sampled by generating random masses, m 1 and m 2 , and initial separation r 0 for the binary PBHs. These values are used to infer an initial semi-major axis r a , and characteristic angular momentum j 0 , which are then used to generate a distribution for the angular momentum, from which a random angular momentum j is drawn, and the eccentricity calculated.

D. Distribution of initial conditions for PBHs merging today
Here we will briefly discuss the typical values for the initial orbital parameters for the binary PBH systems. A full analysis of the entire parameter range will not be considered here (the interested reader can find a more detailed and thorough analysis in Ali-Haïmoud et al. [76] and Raidal et al. [72]) -but a brief summary is included here, because it is useful to consider the typical values which we might expect to find for binary PBHs expected to be merging today, and how these depend on several key parameters.
The first thing to consider is the abundance of PBHs. The higher the number density of PBHs, the closer PBHs will initially be. This means that for large f PBH , we expect to find a smaller initial separation -and this smaller separation means a stronger gravitational binding, meaning an earlier decoupling from the Hubble flow, and an even smaller initial semi-major axis. Therefore, f PBH will have a strong impact on the distribution of the initial semi-major axis.
The mass function will have a similar effect. If the average PBH mass is smaller, this will imply a larger number density of PBHs (assuming the same f PBH ) -and so a higher average PBH mass will typically imply a larger initial semi-major axis.
Here, we make the simple assumption that two PBHs which form close to each other will eventually decouple from the Hubble flow and form a binary, if there are no other PBHs nearby to disrupt this process. This is parameterised by an exclusion radius: we require that there are no other PBHs within a radius Ar 0 (recall that r 0 is the initial separation of the binary). Choosing a higher value for A implies that PBH pairs which form binaries would have a smaller initial separation, and lower initial semi-major axis. A sensible choice for A is likely to be 2 A O(5), and we find that this has a small effect on the initial separation, of order unity. Now, let us turn our attention to the initial angular momentum. After the PBHs form, the PBH pair decouples from the Hubble flow and the PBHs begin to fall back towards each other -before beginning to oscillate around each other. Most of the torque, which provides the initial angular momentum, is expected to originate from the nearest PBH to the fledgling binary -the higher f PBH is, the closer the nearest neigbour is likely to be, and the larger the angular momentum is likely to be.
In addition, the torque is strongest when the nearest neighbour is close relative to the binary separation, before the binary decouples from the Hubble flow -therefore, the longer a system takes to decouple (due to a larger initial separation for example), the higher the total effect of the torque is likely to be. We note that, interestingly, the power-law distribution of j is actually independent of the PBH mass function (although the specific masses of the PBHs in the binary do enter indirectly through j 0 ).
When considering binary systems which are expected to be merging today, a larger initial semimajor axis r a would require a smaller initial angular momentum j and vice-versa (see equation (34) later in the paper). Taken at face value, the simple arguments presented above present a somewhat contradictory picture -a larger f PBH implies r a should be smaller for binaries merging today, whilst also implying j should be smaller. In the end, it is the smaller semi-major axis argument which is more important (owing partly to the fact that a large decrease in r a can be cancelled by a small increase in j to give the same coalescence time).
Depending on the choices of these parameters, typical semi-major axes for PBHs expecting to  (14), or a Gaussian-like distribution, equation (15). We have assumed the following choices for the parameters: masses m 1 = m 2 = m c = 20 M , mass function width σ m = 0.05, and an exclusion zone A = 2 times greater than the initial PBH separation. The abundance of PBHs is described by f PBH , the fraction of dark matter composed of PBHs. Both distributions peak at similar r 0 , but the power-law distribution has a significantly larger tail at small r 0 .  (14), or a Gaussian-like distribution, equation (15). We have used the same parameter choices as previously: masses m 1 = m 2 = m c = 20 M , mass function width σ m = 0.05, and an exclusion zone A = 2 times greater than the initial PBH separation, with the abundance of PBHs is described by f PBH , the fraction of dark matter composed of PBHs. Again, we see that both distributions predict a similar characteristic j, and we see a larger tail for high j for the power-law distribution (corresponding to the low r 0 tail in figure 1). merge today can vary from tens to tens of thousands of AUs, whilst the initial angular momentum can vary from O(10 −3 ) to O(10 −2 ) -meaning extremely high eccentricities, e 0.999. Figure  1 shows the PDF of initial semi-major axes for PBHs expected to merge today (13.7 Gyr after formation of the binary) for different values of f PBH , whilst figure 2 shows the same for the PDF of the initial angular momentum. We have assumed the following choices for the parameters: masses m 1 = m 2 = m c = 20 M , mass function width σ m = 0.05, and an exclusion zone A = 2 times greater than the initial PBH separation.

III. FLY-BYS IN THE SECULAR REGIME
Any object with mass M passing by the PBH binary with mass m ≡ m 1 + m 2 will affect the binary system, potentially breaking it up (e.g., [105][106][107][108][109][110][111][112][113][114][115][116]). In a sub-type of interactions, the third object passes in a wide orbit relative to the binary, which conserves the binary's semimajor axis r a , but induces changes to the angular-momentum and eccentricity vectors. The latter case, known as the secular regime (e.g., [58,80,81,83]), is characterised by the dimensionless quantity (also known as the 'adiabatic ratio') which is the ratio of the perturber's angular speed at periapsis to the binary's mean motion. Here, Q is the perturber's periapsis distance to the binary's center of mass, and E is the eccentricity of the perturber's orbit (E ≥ 1). Note that the perturber eccentricity E can be written in terms of the velocity at infinity v as In our model, encounters are typically highly hyperbolic (E 1). If R 1, this indicates that the secular regime applies, and it is appropriate to average the equations of motion over the binary's orbital phase.
For our fiducial model, the typical adiabatic ratio is where we used that E 1. Since R 1, we can safely assume that the overwhelming majority of perturbers are within the secular regime.
In our Monte Carlo calculations (section V below), we will calculate the effect of the perturbation in the secular approximation using the analytic expressions of [58,84]. Specifically, we include terms of order SA and 2 SA for a given E. Given the excessively large number of individual terms involved and their small values (see Table 1 of [84]), we omit all octupole-order terms associated with 2 SA (the octupole-order terms associated with SA are included).

IV. ESTIMATED EFFECT OF FLY-BYS
In this section, we will provide an analytic estimate of the cumulative effect of nearby PBHs passing near a binary PBH. We will begin by estimating the number density of PBHs residing within a DM halo. For Milky Way-type haloes, the extent of the halo is typically considered to be the region in which the density is 200 times the background matter density of the Universe, where Ω M is the total matter density parameter, and ρ c is the critical density of the universe today. We will take the numerical values, Ω M = 0.315 and ρ c = 1.68 × 10 −23 M au −3 .
For generality, and accounting for the fact that binary PBHs may be found predominantly in PBH clusters with a higher average density, we will consider the halo density to be X halo times greater than the background matter density: wherem PBH is the average PBH mass, and recall that f PBH is the fraction of DM composed of PBHs. For the remainder of this section, we will consider a monochromatic mass function of PBHs, and thus drop the bar notation, such that all PBHs have mass m PBH .
The encounter rate is given, as a function of the encounter radius, by where σ v is the velocity dispersion (which in this section we take to be constant, and all perturbers will move at this speed relative to the binary). The reason for not performing the integral will become apparent soon. Note that we are here ignoring the effect of gravitational focusing, which, in the secular regime being considered, has a negligible effect upon the impact parameter and the distribution of encounters at different radii.
The expected number of encounters is then given by N = Γτ , where τ is the time for which the system is observed, we will use τ = 13 Gyr as a fiducial value: Assuming that all fly-bys are well described by the secular regime, we will use equation (21) above to provide an order of magnitude estimate for the effect of an individual fly-by on the eccentricity of a binary: and note that the effect on the semi-major axis is negligible in the secular regime. Also recall that in this section, we assume all PBHs to have mass m PBH (including the perturbers).
Analytically calculating the total expected change to the eccentricity combined with orbital evolution due to GW emission is complex due to the cumulative nature of the interactions 1 , and so in order to provide a simple order-of-magnitude estimate, we will simply consider the eccentricity to change by a positive SA in each encounter, and ignore the time-evolution of the orbital parameters.
Thus, by equating the impact parameter Q with the radius of an encounter r, an upper limit for the total change in eccentricity (by taking the sum of ∆e from each fly-by, Σ∆e) is given by The integral in the number of encounters now becomes important, because encounters at different radii have a different effect on the eccentricity. We can see that the integral will diverge as r → 0.
This is due to the fact that, as the encounter radius decreases, the chance of an encounter within this radius decreases as r −2 , but the effect increases as r 3 -so even though one does not expect any encounters at small r, the expectation is a large change to the eccentricity from such encounters (although at this point the encounter is no longer secular).
We therefore implement a minimum radius for encounters, such that the expected number of encounters given by equation (25), between 0 and r min , is N = 0.5 -such that most PBHs do not experience a closer encounter than this.
In the Monte Carlo simulations, we require a maximum encounter radius R enc in order to have a finite number of encounters (see section V). However, this is larger, by many orders of magnitude, than r min , and so here we simply take r max → ∞. Accounting for these limits in the integration gives x ≈ 10 −12 r a 10 2 au The coalescence time for a highly eccentric binary PBH system is expected to be proportional to j 7 = (1 − e 2 ) 7/2 . Assuming that the change in eccentricity is small, we can write e = e 0 + x, where the subscript 0 denotes the fiducial eccentricity, and treat x as an expansion variable. To first order in x then, the coalescence time changes by where in the last equality we have made the substitution j 0 = 1 − e 2 0 , and we have neglected the e 0 ≈ 1 term in the numerator. This can be expressed in terms of the model parameters as We can therefore expect that, for the fiducial model, the effect of fly-bys on the encounter rate is negligible. We note that this is only intended as an approximate number, and that many of the parameters are not independent. For example, we saw in section II that the characteristic semi-major axis and eccentricity of the binary orbits are functions of the PBH abundance f PBH and mass m PBH .
We will further investigate the change in lifetime of binary systems by numerically evolving them over time, as well as considering variations from the fiducial model where a significant effect may be seen.

V. BINARY SYSTEM EVOLUTION
We simulate the evolution of a PBH binary after decoupling from the Hubble flow using a Monte Carlo approach. We take into account perturbations from passing PBHs in the secular approximation (see section III), and decay of the orbital energy and angular momentum due to GW emission. Our algorithm consists of the following steps.
• We sample a next perturber given the encounter rate Γ (computed from equation (23)), i.e., the perturber will encounter the binary at a time delay ∆t, where the probability that the time delay exceeds ∆t is given by exp(−Γ∆t). The impact parameter b is sampled from a distribution dN/db ∝ b with 0 < b < R enc . The perturber's mass M is sampled from a lognormal mass function, described in equation (1) • In order for a finite number of encounters to be considered, it is necessary to define a maximum encounter radius R enc , above which we neglect the effect of fly-bys. Since we are considering the effect of an encounter between a binary and a single PBH, we set R enc to a radius around the binary where there is a low probability of finding multiple PBHs: R enc = 0.1 (4πn where n is the number density of PBHs in a DM halo. We do not expect a cut-off in the impact parameter to have a significant effect, as fly-bys with a smaller impact parameter have a larger cumulative effect (see section IV), and we have verified that increasing/decreasing the factor 0.1 by an order of magnitude does not affect the results.
• We apply the effects of the perturber's passage on the binary given M and b using the analytic expressions for ∆e from [58,84]. Here, we assume that the binary's orbital orientation is random (i.e., flat distributions in cos i, ω, and Ω). Evidently, the binary's orientation actually remains fixed between encounters whereas the perturbers plausibly approach from random orientations. However, for the actual computation of ∆e, only the relative orientation matters, so this distinction is unimportant as long as the orientation of the perturbers is isotropic, which is what we assume.
• In-between the current time and the time of the next perturber, we take into account the decay of the orbit due to GW emission by numerically solving the set of ordinary differential equations (ODEs) from [117], i.e., We integrate the above set of equations using the CVODE library [118] in C, and continue until r p = r a (1 − e) < 100 r g , where r g ≡ Gm/c 2 is the binary's gravitational radius. When this condition has been satisfied, we consider the binary to have merged. We note that, by the time of having reached r p = α r g = 100 r g , the binary has mostly circularised; the remaining merger time is [117] T (33)

VI. RESULTS
In this section, we will discuss the results of the simulations of PBH binary evolution. Firstly, we will investigate the merger time calculated by a numeric integration of the equations governing the binary evolution, equation (32), in the absence of encounters which can perturb the system, and then including the effect of perturbers.

A. Numerical evolution of binary primordial black holes
In the limit of high eccentricity, e → 1, an analytic expression for the coalescence time τ A for binary systems is given by [117] τ A = 3 85 Since binary PBHs themselves are typically highly eccentric, this equation has been typically been used to calculate the coalescence time of binary PBH systems (i.e. [72,76]), and is a key component in predicting the merger rate observable today.
In this section, we test the accuracy of this approximation by numerically integrating equation (32), as described in section V, for PBH binaries with random initial conditions to give a numeric value for the coalescence time τ N , and comparing this to the time τ A calculated from equation (34). We then consider the effect this might have on the observable merger rate today, finding that, whilst equation (34) can be inaccurate to O(10%), this is likely to have a negligible effect on the merger rate.
Our initial PBH binaries are randomly generated using the methods described in section II. The relevant parameters which will affect the distribution of initial conditions are: the peak mass of the PBH mass function m c , the width of the mass function σ m , the size of the exclusion zone required for a binary to form A, and the total abundance of PBHs, parameterized by f PBH . The great majority of binary systems merge early during the history of the universe, τ < 1 Gyr. Since such systems are not relevant for the observation of mergers happening today, we limit our selection of binary PBHs to those with a lifetime greater than 1 Gyr.
If PBHs make up the entirety of DM, f PBH = 1, we find that the analytic expressions only matches the numeric results to O(10%), but matches much more closely, to within O(1%) when PBHs are less abundant, f PBH ≈ 0.01. This is likely due to the fact that, whilst all long-lived PBH binaries are highly eccentric, they are significantly more eccentric for low f PBH (see section II for more discussion). We find that changing µ, σ m and A does have an effect, this effect is subdominant to the effect of changing f PBH . We therefore fix µ = 20 M , σ m = 0.5 and A = 2 for the remainder of the discussion in this section.
We define the quantity δ τ as the relative change in τ At first glance, this may be expected to have a significant effect on the merger rate observable today, since the merger rate is expected to decrease over time. To estimate how much the merger rate is affected, we will compare the number of binaries merging over time when the coalescence times is calculated numerically or analytically. To achieve this, we use 100 000 initial binaries (with a lifetime greater than 1 Gyr) and calculate their coalescence time, and sort the coalescence times into 1 Gyr bins.  [72] predicts that the merger rate R follows a power-law with respect to time t, R ∝ t −34/37 , so we fit a power-law to the data points, shown with a red (blue) lines.
We see only a small deviation in the number of PBHs merging at any given time -which is consistent with the random errors caused by the finite sample size. We therefore conclude that, whilst the analytic expression for the merger time can have a significant error (especially for individual binaries), this has a negligible impact on the merger rate predicted -and may therefore be safely ignored. Finally, we turn our attention to the amplitude of the numerical error made when integrating (32). Numerical errors made in the ODE integration can lead to errors in the coalescence time compared to the true solution. It is necessary to quantify these errors when considering the effects of fly-bys as well, in order to be able to properly distinguish the physical effects of fly-bys from unphysical numerical noise. We find that the largest numerical error made in our integrations gives rise to δ τ of order 10 −8 , which, as we will see, is safely several orders of magnitude smaller than the observed signal due to fly-bys.

B. The effect of fly-bys
In this section we will discuss the effect of distant encounters between binary PBH systems and individual PBHs -referred to as fly-bys. We begin with a random sample of PBH binaries, before submitting them to random encounters with other PBHs. To make the task manageable, we make the simple assumption that such binary PBHs are part of a DM halo (and/or PBH cluster) from shortly after the time of formation until today -and neglect any time evolution of such haloes/clusters. A more complete investigation of the evolution of binary systems from the early universe, through structure formation, up until the present epoch will require the use of N -body simulations with PBH DM.
There are a large number of parameters which can affect the typical coalescence times of binaries, and so we will investigate a fiducial model with fixed values for these parameters, before discussing the effect of individual parameters. For the fiducial model, we take the following values for the required parameters: • m c = 20 M , σ m = 0.5, close to the values given by Raidal et al. [72] as a best-fit model to explain the black holes coalescence events observed by LIGO.
• f PBH = 1 (although we note that this is expected to produce a larger frequency of coalescence events than are observed).
• A = 2, throughout we will make the simple assumption that if two PBHs form within a distance x of each other, they will eventually form a binary if there are no other PBHs within a distance 2x. Changing this by order unity has a similar order effect on the distribution of initial semi-major axes, and a negligible effect on our final results.
• X halo = 200, corresponding to the usual definition of DM halos. The great majority of DM in the universe is expected to be found inside halos with a density greater than 200 times the background density of the universe -though we note that binary PBHs may be found in denser sub-haloes, and the effect of mass segregation may mean that the heavier binary systems migrate towards the denser cores of DM halos.
• σ v = 200 km s −1 , corresponding approximately to the velocity dispersion for a Milky Waytype halo calculated using virial theorem. We take a zeroth order estimate for σ, but the actual velocity dispersion within a DM halo is position dependant, and calculating the relative velocity of binaries with respect to nearby PBHs goes beyond the scope of this work.
We will compare the coalescence time calculated with and without encounters from the Monte Carlo simulations of binary PBHs. We will again use δ τ to represent the fractional change in the merger time for PBHs, where τ MC is the coalescence time given by the Monte Carlo simulations including fly-bys, and τ N is the coalescence time predicted from numerically solving equation (32) from the initial conditions of the binary.  As may have been expected from the divergence in the integral in equation (27), the standard deviation is dominated by outliers in the distribution, a small number of binaries experience a change in merger time orders of magnitude larger than typical (up to δ τ = O(0.01) for the binaries considered). Thus, in order to describe the typical effect of fly-bys with a more representative number, and to investigate the effect of changing the model parameters, we will neglect the tails of the distribution, retaining only the central 90% of the data 3 . In this case, we find a significantly smaller standard deviation σ * τ = 4.64 × 10 −6 (and the mean is δ * τ = −5.81 × 10 −8 ).

Changing the model parameters
We now ask the question of how the coalescence time is affected for different parameter choices.
Starting from 10 000 initial binaries in each case, figure 6 shows how σ * τ is affected by changing the model parameters. Here, to increase the sample size, we consider all binaries which merge within we consider the PBH abundance f PBH , seeing that fly-bys have a larger effect for smaller f PBH , due to the increased characteristic semi-major axis of binaries in that case. For illustration purposes, we have plotted linear lines of best fit for m c , σ m , and X halo , and power-law lines of best fit for σ v and f PBH .
shows that σ * τ only has a weak dependence on the mass function 4 . On the other hand, changing the parameters f PBH , X halo or σ v has a relatively strong effect.
Decreasing f PBH might be expected to result in fly-bys having a smaller effect, since there are less PBHs in a DM halo to interact with a given binary. However, as described in section II, a smaller f PBH also implies that binary PBHs will form with significantly larger semi-major axes (and those merging by today also had a higher initial eccentricity) -meaning that encounters have a larger effect on the orbital parameters.
Changing X halo is found to have a significant effect on σ * τ . If we consider that PBH binaries might be found within denser DM haloes (and neglecting the fact that this is likely to change the velocity dispersion), this has the simple effect of increasing the number of encounters that occurwithout changing the nature of the encounters. We find the expected result, therefore, that σ * τ is approximately proportionate to X halo -the effect is not exactly proportionate due to the cut-off in the maximum encounter radius considered.
It might naively be expected that decreasing the velocity dispersion σ v will lead to a smaller effect on the coalescence time -since this implies there will be fewer encounters. However, each fly-by takes a lot longer to occur, and thus has an overall larger impact as σ v decreases (see equation (31)) -and we therefore see that decreasing σ v from 200 km s −1 to 10 km s −1 makes the effect of fly-bys significantly larger.
We now consider the largest value for σ τ which might be obtained for reasonable choices of the model parameters. We have seen that changing the mass function parameters does not strongly affect our results, so we will keep m c = 20 M and σ m = 0.5 fixed. Constraints on PBH abundance arising from GW signals from merging PBHs are of order 10 −3 , so we will take f PBH = 5 × 10 −3 .
For the parameters related to DM haloes, we will consider that binaries might be more likely to be found in denser regions (for example, due to the formation of PBH clusters [119]), where not only is the density higher, but the velocity dispersion is likely to be lower [120]. We will therefore consider X halo = 1000 and σ v = 10 km s −1 . With these choices of parameters, we find σ τ = 8.60 × 10 −2 (and σ * τ = 2.62 × 10 −4 ).

Non-secular encounters in the low f PBH regime
As f PBH becomes small, we find that fly-bys have a larger effect on the coalescence time. This is due to the fact that, while typical semi-major axes of binaries and impact parameters of encounters both increase, the semi-major axis of binaries typically grow by more -resulting in fly-bys having This also means that the chance of non-secular encounters ('strong encounters') increases, i.e. encounters with an adiabatic ratio R 1. Our current formalism is not capable of dealing with such "strong encounters", and when a strong encounter does occur in the Monte Carlo simulation, this is recorded, and evolution of that binary is halted. For f PBH = O(1), the number of binaries experiencing a strong encounter is negligible. However, for f PBH = O(10 −3 ), we find that ∼ 10% of binaries experience a strong encounter -which is actually a larger fraction than the number of binaries merging between 10 − 15 Gyr.
However, the semi-major axes of binaries experiencing a strong encounter are orders of magnitude higher than the semi-major axes of binaries which merge within the simulation time. Figure   7 shows a comparison of the semi-major axes of binaries which merge in the simulation, compared to those which experience a strong encounter (for the parameter choices given at the end of the previous subsection). For binaries which merge are expected to merge in the current lifetime of the universe, the semi-major axes of order 10 3 au, whilst for those experience strong encounters, the semi-major axes are of order 10 6 au.
Neglecting such binaries is therefore unlikely to affect the result and conclusions presented in this paper -although it is conceivable that strong encounters may drive binary systems to merge much earlier, thereby increasing the merger rate observed today. We also note that such wide binaries may have been disrupted in the early universe shortly after formation. We leave further consideration of such binaries for future work.

VII. DISCUSSION
We have considered PBH binaries which form in the early universe, shortly after the formation of the PBHs themselves. Utilising a Monte Carlo approach, the evolution of a large sample of initial binaries are numerically evolved forwards through time in order to determine their coalescence time.
We account for the impact of encounters of fly-bys in DM haloes (whilst binary-binary encounters are rarer and have a negligible effect in the secular limit compared to binary-single encounters [121]).
In section IV we also developed an analytic estimate for the change in coalescence time due to the effect of fly-bys, and in later sections we investigated the effect with a Monte Carlo approach.
We typically find that binary-single encounters have a small effect on the coalescence, changing the lifetime by order 10 −6 , although a small number of binaries experience a much larger effect.
Considering a more extreme model, we find that the typical binary lifetime is unlikely to change by more than 10 −2 .
We note that we have neglected the formation of dense DM haloes ("spikes") around PBHs in models where PBHs only constitute a small fraction of dark matter (i.e. [122,123]). In addition to increasing the effective mass of PBHs, these DM spikes will also affect the in-spiral of binary PBHs and their gravitational waveforms [124]. A PBH (with a DM spike) flying by a binary system is far enough away to be considered as a point mass, and since the important factor is the ratio of perturber mass to binary masses, the change in mass due to the DM spikes is unlikely to have a signficant effect on the outcome of the perturber. Neglecting the the DM spikes surrounding PBHs is therefore expected to have a negligible effect on our conclusions.
We also compared numerically solving the equations governing binary evolution to an analytic approximation used in previous studies, equation (34). By using a numeric method we obtain a more accurate value for the coalescence time than using the analytic expression. We find whilst that the analytic expression generally overestimates the coalescence time for individual systems by around 10%, this is found to have a negligible impact on the merger rate of binaries, and so may safely be neglected.
We therefore conclude that neglecting the impact of binary-single encounters after formation of a binary is unlikely to result in significant error to the coalescence time, placing the constraints on PBH abundance arising from the detected GW signals from merging PBHs on a more secure footing. However, we note that further work may be necessary before this can be stated with certainty -our results show a need to study the evolution of PBH binaries in the early universe, during structure formation, and the dynamics of PBHs within DM halos.
A recent paper by Jedamzik [125] performed a similar analysis to the one presented here, analysing the effect of interactions between a binary PBH and a third by-passing PBH. The results presented there are complimentary to those which we present here, and concern with the evolution of binaries in extremely dense clusters (which are denser by many orders of magnitude than the DM haloes considered here), which form at high redshifts but evaporate by lower redshifts. The conclusion reached in that paper is that many binaries in such clusters are disrupted, and that the merger rate observed today is consistent with PBHs composing the entirety of dark matter.