Lower bound on the primordial black hole merger rate

We derive a lower bound on the merger rate of primordial black hole (PBH) binaries by estimating the maximal fraction of binaries that were perturbed between formation in the early Universe and merger, and computing a conservative merger rate of perturbed binaries. This implies robust constraints on the PBH abundance in the range 1–100M⊙. We further show that LIGO/Virgo design sensitivity has the potential to reach the PBH mass range of 10–10M⊙. The constraint from the merger rate of perturbed binaries is stronger if PBHs are initially spatially clustered.


I. INTRODUCTION
The direct detection of gravitational waves (GWs) from the binary black hole (BH) merger by LIGO [1] marked the dawn of a new era in cosmology. During the first two observing runs LIGO detected nine more BH-BH merger events, latest of which were seen also by the Virgo detector, indicating a merger rate of 9.7-101 Gpc −3 yr −1 for BH binaries with component masses ranging from 7.6M ⊙ to 50.5M ⊙ [2]. Currently possible new events are reported almost on a weekly basis [3]. These events may have originated from astrophysical BHs [4] or from primordial black holes (PBHs) [5][6][7] that, in the most common scenarios, arise from large curvature fluctuations in the early Universe [8,9].
PBH binary formation has received much attention in the last years [5][6][7][10][11][12][13][14][15][16][17][18]. Most PBH binaries whose mergers would be presently observed form already in the radiation dominated Universe [7,13,14]. Due to the random distribution of PBHs before formation of structures, some PBHs may be much closer than their average separation. Such PBHs will form the earliest gravitationally bound systems and produce the initial population of binaries. The distribution of orbital characteristics of this binary population can be estimated analytically and, assuming all binaries remain unperturbed until merger, the merger rate obtained via this mechanism is so large that only a small fraction, Oð0.1%Þ, of dark matter (DM) in PBHs is allowed in the mass range that LIGO/Virgo detectors probe. This implies the strongest constraints on the PBH abundance to date in the mass range 1-100M ⊙ [13,14,17]. Yet, the fate of these binaries remains uncertain. In particular, initial binaries are highly eccentric, thus interactions with surrounding PBHs can significantly increase their coalescence times.
The aim of this paper is to derive a lower bound on the PBH merger rate by considering scenarios where the initial population of binaries is maximally perturbed and by including the contribution from perturbed binaries. These results imply more reliable constraints on the PBH abundance from GW observations. Numerical studies, although confirming the analytical orbital parameter distribution of initial binaries, find that binaries are likely to be perturbed after formation in the case when PBHs make up most of the DM [17]. 1 The binaries can be perturbed by two mechanisms: (1) In case the initial configuration contains a third PBH close to the PBH pair that is expected to form a binary, it is very likely that it collides with the binary. (2) PBHs will form dense N-body systems relatively early, and binaries absorbed by these clusters become more likely to be perturbed.
In [17] we introduced a suppression factor of the merger rate due to disruption by the first mechanism. In this paper we consider the disruption via the second mechanism. Analytic estimates of Ref. [14] indicate that most of the initial PBH binaries are not perturbed between formation 1 We consider disruption (or perturbation) of binaries in a very broad sense: This includes modification of the orbital parameters due to tidal torque, hard collisions with other BH resulting in a binary (including the ones where a BH in initial the binary is swapped), and ionization of the binaries. For hard binaries, only the first two scenarios are probable. and merger in typical DM haloes larger than 10 PBH, assuming that early haloes are disrupted when absorbed by subsequent larger haloes. We show that the dense haloes forming at redshifts z ≳ 100 can significantly perturb their initial binary population in case they survive absorption by their later hosts. Furthermore, self-gravitating systems are not stable, and, due to the gravothermal catastrophe [19,20], the binary may end up in an environment where the PBH density is several orders of magnitude larger than what was used in Ref. [14], for example, when a binary finds itself in the central region of a PBH cluster undergoing a core collapse. This increases the rate of close encounters and the probability of perturbing the binary. To obtain the strongest possible suppression for the merger rate of initial binaries we assume that clusters small enough to experience a core collapse within a Hubble time will perturb all of their initial binaries. Larger haloes are stable and do not perturb their initial binaries due to their smaller densities and higher velocity dispersions.
In this paper we also improve our previous estimate [17] of the merger rate arising from the population of perturbed binaries. Even if most of the early PBH binaries were at least once perturbed, the present merger rate from the resulting population of less eccentric binaries can still reach the observed rate. In that case, binaries contributing to the rate originate from dense 3-body configurations with separations much below the average PBH distance.
By combining the merger rate of perturbed binaries with a maximally suppressed merger rate of initial binaries we derive a lower bound for the merger rate, concluding that it is above the range indicated by the LIGO and Virgo observations in the case that more than 4% of DM is in Oð10M ⊙ Þ PBHs. We finally use our results to obtain robust constraints on PBH abundance. Throughout this paper we use geometric units G N ¼ c ¼ 1.

II. PBH BINARIES FROM THE EARLY UNIVERSE
The coalescence time of a binary due to GW emission is approximately [21] where η ¼ m 1 m 2 =M 2 , M ¼ m 1 þ m 2 denote the mass asymmetry and the total mass of the binary, r a is its semimajor axis, E ¼ M=ð2r a Þ the binding energy per reduced mass, j ≡ ffiffiffiffiffiffiffiffiffiffiffiffi ffi 1 − e 2 p the dimensionless angular momentum, and e is the eccentricity.
Binaries expected to coalesce within a Hubble time are hard. Thus, according to the Heggie-Hills law [22,23], close encounters with other PBHs will, on average, increase their binding energy by an Oð1Þ factor. In particular, these binaries are very unlikely to be ionized. On the other hand, the initial binaries are highly eccentric. The characteristic angular momentum for initial binaries with coalescence time τ is [17] j τ ¼ 0.02f 16 37 PBH ð4ηÞ where t 0 is the age of the Universe. Thus, j very likely increases by more than an order of magnitude leading to an over 7 orders of magnitude increase of the coalescence time. So, when a binary that was initially expected to merge within a Hubble time is perturbed, its coalescence time exceeds the age of the Universe, and thus it will not produce detectable GW signals. We stress that, due to their lower eccentricities, the last effect is insignificant for perturbed binaries.
If the initial population would not be perturbed by the second process, then a fraction f PBH ≥ 10 −3 of PBH DM predicts a merger rate higher than observed by LIGO/Virgo. The differential merger rate at time t in that case is [17] dR np dm 1 dm 2 ≈ 1.6 × 10 6 Gpc 3 yr ψðmÞ is the PBH mass distribution, and S½ψ; f PBH ; M is a suppression factor accounting for the disruption by the first mechanism due to the infall of the PBH close to the binary. For narrow mass functions S ≈ 0.24ð1 þ 2.3σ 2 M =f 2 PBH Þ −21=74 , where σ M ≃ 0.005 is the variance of matter density perturbations at the time the binary was formed. In the rest of the paper we will consider only monochromatic mass functions, ψðm 0 Þ ¼ δðm 0 − mÞ, and denote the PBH mass by m.
When f PBH ≳ 0.1, the fraction of binaries perturbed by the second mechanism, that is, by PBH clusters, has been shown to be relatively high already at z ¼ 1100, indicating that nearly all initial binaries might be perturbed within the age of the Universe [17]. A small fraction of the early binaries may, however, remain unperturbed. So, by (3), the present merger rate from the initial binary population can still exceed the current bounds in the LIGO/Virgo mass range.
This may happen because of several reasons: (a) not all PBH binaries will become bound to PBH clusters early; (b) small dense clusters of PBH are unstable on timescales much shorter than the Hubble time, and thus it is possible that the cluster is dissolved before the binary is perturbed; (c) larger, less dense systems forming the DM haloes of, e.g., dwarf galaxies, must be stable within the Hubble time making disruption unlikely.

III. DISRUPTION OF INITIAL PBH BINARIES IN SMALL HALOES
Let us estimate the probability that a binary will be significantly perturbed via encounters with other PBH. The characteristic timescale for this process is given by where n loc is the local PBH number density, v the perturber velocity, and σ Δj>j τ the cross-section for increasing the angular momentum of the initial binary by an amount comparable to its initial value (2). Note that, by Eq. (1), this corresponds to a 2 order of magnitude increase in the coalescence time. The average change of angular momentum is roughly Δj ≈ m 1=2 r 3=2 a =ðr c bvÞ, where b is the impact parameter and r c is the distance of closest approach [14]. By conservation of angular momentum and energy, they are related as b 2 ¼ r 2 c þ 6mr c =v 2 . The second term dominates in the early Universe where close encounters are more likely due to the lower velocities and higher densities. In this case Consider now the early haloes containing N PBHs. They form approximately when the scale factor is a c ≡ð1þz c Þ −1 ≈ a eq ffiffiffiffi N p =f PBH , where a eq corresponds to matter-radiation equality. Assuming that they are virialized, the velocity dispersion is given by where R is the virial radius and M H ¼ mN=f PBH the mass of the halo. 3 Following Press-Schechter theory, the average energy density of matter in these haloes is ρ ¼ 3M H =4πR 3 ≈ 18π 2 ρ c a −3 , where ρ c is the critical comoving density.
Using v ≈ σ v and a density that is 18π 2 times above the average, we obtain from (4) and (5) that Therefore, when m ¼ Oð10M ⊙ Þ binaries in haloes with N ≲ 2000 are perturbed before today if f PBH ¼ 1. Since, by definition, binaries can be perturbed only in systems with N c ≥ 3 PBH, binaries are unlikely to be perturbed when f PBH ≲ 0.02. These estimates are applicable assuming that the small haloes formed at high redshifts survive, i.e., they maintain their initial density and velocity dispersion. This assumption is, however, easily violated as the haloes can be absorbed into larger structures or expand due to the heating provided by binary-PBH collisions. Both effects reduce the frequency of binary-PBH encounters. In fact, it was estimated that the disruption of initial binaries in DM haloes can be neglected when the earlier smaller haloes are continuously absorbed and disrupted by the subsequent generation of larger haloes [14].
By drawing parallels with globular clusters [25,26], there are two effects that can significantly enhance the disruption probability in PBH haloes: First, since binaries are heavier than a single PBH, they tend to sink towards the center of the halo. Second, during core collapse the central density can increase indefinitely until the collapse is stopped by binary-PBH interactions heating the core. [20]. This will switch on 3-body encounters also in larger haloes in which they were initially unlikely. Thus, to estimate the maximal disruption probability by the second mechanism we assume that all binaries within a cluster that is unstable within a timescale less than a Hubble time are perturbed.
For the timescale of the gravothermal instability we use the characteristic time of core collapse given by t cc ≥ 18t r , where the relaxation time is [27] The Coulomb logarithm is approximately ln Λ≈ lnðN=f PBH Þ. Requiring 18t r < t 0 gives which is consistent with earlier results of [28]. We now find that binaries are perturbed in haloes with N < 5300 if f PBH ¼ 1 and all binaries survive if f PBH ≲ 0.005. We note that binary-PBH collisions can heat the halo and stop the collapse; thus, disruption of initial binaries due to core collapse is relevant for haloes containing 2000 ≲ N ≲ 5300 PBH when f PBH ≈ 1.
Finding the probability for an initial PBH binary to be disrupted thus boils down to finding the fraction of initial binaries in unstable haloes. The distribution of haloes containing N PBHs at redshift z from initially Poisson distributed point masses is approximately [24,29,30] where N Ã ðzÞ is the characteristic number of PBH in a halo at redshift z, which we estimate using analytic results from [24]. These haloes can have substructure which, following [31], we assume to be distributed by the halo mass function (9). The probability of finding a binary in a halo with N PBHs is approximately proportional to p N 4 and the 3 We assume that the fraction of PBHs in haloes matches f PBH . However, it was shown in [24] that this fraction could be larger, especially in the early Universe. This may slightly enhance the disruption of PBH binaries. 4 We remark, that binaries for which the surrounding PBH density was initially larger have, on average, a higher initial angular momentum and must thus have a smaller initial separation if they are to merge within a Hubble time. Therefore, there are fewer initial binaries in dense haloes. Such haloes are likely small. probability of finding a PBH in a subhalo of N PBHs inside a halo of N 0 > N PBHs to p N p N 0 . So, the fraction of nonperturbed binaries is bounded below by where N c ðzÞ defined in Eq. (8) is the smallest number of PBH in haloes or subhaloes that are expected to be stable until redshift z, and the halo distribution is evaluated at the redshift z c at which the haloes with N c PBHs formed. The probabilitiesp N andp N are normalized as Here we assume that all subhaloes survive and that every PBH inside a halo belongs to some subhalo. As substructure can be absorbed or disrupted by the host, this construction will overestimate the probability for the initial binary to be perturbed leading to a conservative merger rate for initial binaries.
In Fig. 1 we show the strongest suppression factor P np ðzÞ of the merger rate. For f PBH ¼ 1, at most 6% of the initial binaries survive unperturbed until today. Moreover, disruption by clusters can be observed at very high redshifts. This is consistent with our N-body simulations [17], according to which only about half the binaries not perturbed by the first mechanism survive until z ≃ 1100 when f PBH ¼ 1.

IV. MERGER RATE OF PERTURBED BINARIES
Even if all early binaries would be perturbed, a large population of PBH binaries capable of contributing to the present merger rate still remains. In the following we will estimate the merger rate of perturbed binaries.
Consider first the initial population of binaries and the resulting merger rate. Following [17] we approximate that all initial pairs are perturbed when there is a third PBH closer than y min , determined by Nðy min Þ ¼ 2, where NðyÞ ≡ 4πn PBH y 3 =3 is the expected number of PBH in the comoving volume sphere of radius y and n PBH is the PBH number density. This comprises 76% of all initial PBH binaries (for f PBH ≳ 1.5σ M ≈ 0.01) and the rate (3) results from the remaining 24%. We stress that the initial binary may be perturbed even if the third PBH is farther than y min by the mechanism discussed in the previous section.
Most of the perturbed PBH binaries whose mergers can be observed today would come from PBH pairs that initially formed an eccentric binary with a very short coalescence time. It is thus necessary to impose that the initial coalescence time (1) is larger than the time required to perturb the binary, τ i ≫ t p . The first close encounter of the initial pair takes place at a ≈ a eq NðxÞ=ð2f PBH Þ, where x denotes the initial comoving separation of the pair. Analogously, the initial encounter of the binary with the third PBH takes place at a ¼ aðt p Þ ≈ a eq NðyÞ=ð3f PBH Þ. This estimate works well when the 3-body system is only weakly coupled to the surrounding PBHs, implying NðyÞ ≲ 1. For larger values of y, it is more likely that the binary is not the closest object to the third PBH. If the 3-body system forms during radiation domination, we obtain For short coalescence times GW emission may be relevant. The energy of eccentric orbits scales as EðtÞ ¼ E 0 ð1 − t=τÞ −2 [21]. Thus, to avoid binaries that emit more than 10% of their initial binding energy in GWs, we consider a population in which the coalescence time of the initial binary satisfies In the matter dominated epoch (12) overestimates t p , and therefore requiring (13) leads to a lower merger rate estimate. The initial coalescence time depends on the initial angular momentum which, assuming a hierarchical collapse, is determined by the tidal forces acting on the pair. The closest PBH generates an angular momentum j 1 ¼ 0.7jsinð2θÞjNðxÞ=NðyÞ, where θ is the angle between the vector joining the PBH pair and the vector joining the center of mass of the pair and the position of the third PBH. The angular momentum generated by tidal forces from all other surrounding PBHs is of the order j 0 ≈ 0.5NðxÞ, where the width of the distribution is σ j ≈ 0.5NðxÞ= ffiffiffiffiffiffiffiffiffiffi NðyÞ p [17]. If NðyÞ≪1, then j 1 ≳σ j ≫j 0 as long as jsinð2θÞj ≳ ffiffiffiffiffiffiffiffiffiffi NðyÞ p . In this case the contribution from other PBHs is negligible and we can use the 3-body approximation. The comoving density of binaries with an initial binding energy in the interval (E, E þ dE) is where is the density of initial configurations specified by x, y, and θ, and is the probability that such initial configurations produce a binary with binding energy E satisfying (13). Θ denotes the step function. The relation between binding energy and initial separation is approximately ; it is independent of the mass of the PBHs and proportional to f 3 PBH . After the initial binaries have interacted with the surrounding PBHs their orbital parameter distribution is approximately where ∂P p ðjÞ=∂j is the angular momentum distribution of the perturbed binaries and KðEjE 0 Þ is the energy distribution of perturbed binaries with initial energy E 0 . The assumption that the angular momentum distribution PðjÞ of perturbed binaries is independent of the initial angular momentum agrees with numerical studies of binary-single PBH collisions [32,33]. The thermal distribution of angular momenta is dP ¼ 2jdj [34]. However, numerical studies find that the j distribution for perturbed binaries in the early Universe is dP ¼ dj [17]. We will consider distributions with γ ∈ ½1; 2. Numerical results of Ref. [17] further imply that where α > 0. Binding energies do not change in the limiting case α → ∞, which gives a lower bound on the merger rate. Being perturbed by the initially closest PBH corresponds to α ≈ 1 [17]. However, binaries may also be perturbed later. Instead of suppressing the merger rate, as in the case of initial binaries, this will lead to an increase in the merger rate because the binaries get harder, i.e., α effectively decreases, while their eccentricity distribution is approximately preserved. The present perturbed binary merger rate can thus be enhanced already when f PBH ≈ 1% as the disruption rate may exceed 50% in that case [see Fig. 1]. We neglect binary-binary collisions, which are likely to ionize the wider binary [35]. Finally, the merger rate of perturbed binaries at time t is When γ > 21=37, this rate grows slower with redshift than the merger rate of unperturbed binaries and scales down faster for decreasing PBH abundance than the rate from unperturbed binaries when γ > 7=24 and f PBH ≳ 0.01. It depends weakly on the lower bound of the initial coalescence time given by (13), i.e., requiring τ i > qt p ðyÞ, we find that R p ∝ q 16ð28−9γÞ 2331 . The present merger rate of perturbed binaries is shown in Fig. 2 for m ¼ 20M ⊙ . The merger rate decreases as a function of α as seen in the upper panel, and the lower panel shows the dependence of the merger rate on f PBH . The merger rate is a decreasing a function of γ, that is, distributions with more eccentric perturbed binaries correspond to a larger merger rate. 5 When f PBH ≈ 1, the merger rate can be dominated by perturbed binaries. For f PBH ≳ 0.3 the merger rate R p exceeds the region observed by the LIGO/Virgo collaboration thus ruling out m ¼ 20M ⊙ PBH DM even when all binaries are perturbed. 5 In [13] the rate R p was found to be about 2 orders of magnitude smaller because it was estimated assuming the extremal case where all perturbed binaries have circular orbits. This corresponds to the limit γ → ∞.

V. DISCUSSION
Consider the modifications due to initial spatial clustering of PBHs [36]. The effect of a nonvanishing 2-point function ξ PBH ðxÞ on R np in typical models of inflationary PBH production is estimated to be small [15]. It should, however, be reconsidered for R p . The spatial correlation of PBH can be included by using NðxÞ ¼ 4π R dxx 2 ð1 þ ξ PBH ðxÞÞ in (15). To find a rough quantitative estimate, we assume a constant 2-point function at scales relevant to binary formation, 1 þ ξ PBH ðxÞ ≈ δ PBH when NðxÞ ≲ 1 [13]. This is equivalent to a local change in the PBH abundance, so its effect on the merger rate amounts to a simple rescaling Thus, R p ∝ δ 144γ 259 þ 10 37 PBH is more sensitive to clustering than R np ∝ δ 16=37 PBH . Moreover, clustering will generally lead to more of the initial binaries being perturbed, i.e., a smaller P np , and thus to a reduction of R np . On the other hand, R p will be enhanced as more frequent 3-body encounters tend to harden the perturbed binaries. Initially clustering therefore enhances R p compared to R np .
We expect our results to hold for narrow mass functions. The extension for wider mass functions is nontrivial due to inherent nonlinearities. Mass segregation will enhance the disruption of heavy initial binaries, especially during core collapse, as it forces them to migrate towards the center of the halo. The opposite holds for light binaries, because, although they are more easily perturbed due to their generally lower binding energy, they tend to migrate away from the center, decreasing the probability of hard collisions. Extended mass functions also shorten the collapse timescale [37], which will make perturbing the binaries more likely. Finally, perturbed binaries are expected to contain PBHs from the heavy end of the mass distribution as the lightest PBH gets ejected in 3-body encounters [35].

VI. A CONSERVATIVE CONSTRAINT ON PBH ABUNDANCE
Constraints arising from the LIGO/Virgo GW measurements for monochromatic PBH mass function are shown in Fig. 3. For these we used the merger rate R ¼ P np R np þ R p , where for the parameters of R p we used the values α ¼ 1 and γ ¼ 2 that yield a smaller rate. Including the rate R p has a relatively small impact on the constraints. The red line shows the 2σ constraint on the PBH abundance corresponding to the case that none of the BH merger events  (2), is neglected. The black dashed contours depict the likelihood fit for a log-normal PBH mass function with width σ ¼ 0.4 on the observed rate and masses. originated from PBH binaries. This is obtained by calculating the expected number of events [17], where T ¼ 165 days is the observation time for LIGO in observation runs O1 and O2, dV c ðzÞ is the differential comoving volume element, ρðm 1 ; m 2 ; zÞ 2 is the signal-tonoise ratio for LIGO observing a merger of BHs of mass m 1 and m 2 at redshift z [38], and ρ c ¼ 8 is the threshold value for detectability of the event [38]. The bounds shown by the red lines correspond to N < 3 which, assuming that the events are Poisson distributed, gives the 95% (or 2σ) upper limit on f PBH . The scenario where all observed BH merger events arise from PBHs takes place in a narrow mass range. For this case we make a maximum likelihood fit of the PBH abundance and mass function (see [17] for details). The 2σ and 3σ contours corresponding to a fit of the observed LIGO/Virgo events to the PBH scenario with a narrow lognormal mass function are shown by the dashed black lines and the best fit by the black dot.

VII. CONCLUSION
In conclusion, we computed the lower bound on the merger rate of PBH binaries by estimating the maximal suppression of the initial PBH binary merger rate. To this aim we assumed that (a) early haloes survive when they are later absorbed by larger structures and (b) haloes small enough to be gravitationally unstable within a the Hubble time will disrupt their binaries with certainty. We remark that omitting assumption (b) would still imply a relatively large suppression as early haloes, due to their larger densities and smaller velocity dispersions, can disrupt their binaries within a timescale much shorter than the age of the Universe. We further gave a conservative estimate of the merger rate from a population of PBH binaries that were formed in the early Universe from dense 3-body systems and found that the resulting merger rate can exceed the one observed by LIGO/Virgo when f PBH ≳ 0.3.
These results put the GW constraints on the PBH abundance on a stronger footing. In particular, scenarios where PBHs make up all DM are ruled out in the range 1-100M ⊙ and LIGO/Virgo design sensitivity has the potential to probe the wide mass range of 10 −2 -10 3 M ⊙ . We find that the PBH scenarios for the LIGO/Virgo GW events could be realized by relatively narrow PBH mass functions centered around 20M ⊙ with the PBH abundance ranging from 0.1%-4%. Our conclusions persist in models where PBHs are initially clustered as, although clustering makes it more likely for initial binaries to be perturbed, it increases the merger rate of perturbed binaries. As our results are based on analytic models of nonlinear PBH structure formation and binary-PBH interactions, for our conclusions to be considered definitive they should be fully and rigorously tested with numerical simulations, which should predict merger rates higher than found here.