key: cord-0647863-k0zgipod
authors: Jedamzik, Karsten
title: Primordial Black Hole Dark Matter and the LIGO/Virgo observations
date: 2020-06-19
journal: nan
DOI: nan
sha: a0d6dbde1455a17501f2183c103e857d0d389c37
doc_id: 647863
cord_uid: k0zgipod

The LIGO/Virgo collaboration have by now detected the mergers of ten black hole binaries via the emission of gravitational radiation. The hypothesis that these black holes have formed during the cosmic QCD epoch and make up all of the cosmic dark matter, has been rejected by many authors reasoning that, among other constraints, primordial black hole (PBH) dark matter would lead to orders of magnitude larger merger rates than observed. We revisit the calculation of the present PBH merger rate. Solar mass PBHs form clusters at fairly high redshifts, which evaporate at lower redshifts. We consider in detail the evolution of binary properties in such clusters due to three-body interactions between the two PBH binary members and a third by-passing PBH, for the first time, by full numerical integration. A Monte-Carlo analysis shows that formerly predicted merger rates are reduced by orders of magnitude due to such interactions. The natural prediction of PBH dark matter formed during the QCD epoch yields a pronounced peak around $1M_{odot}$ with a small mass fraction of PBHs on a shoulder around $30M_{odot}$, dictated by the well-determined equation of state during the QCD epoch. We employ this fact to make a tentative prediction of the merger rate of $sim 30M_{odot}$ PBH binaries, and find it very close to that determined by LIGO/Virgo. Furthermore we show that current LIGO/Virgo limits on the existence of $sim M_{odot}$ binaries do not exclude QCD PBHs to make up all of the cosmic dark matter. Neither do constraints on QCD PBHs from the stochastic gravitational background, pre-recombination accretion, or dwarf galaxies pose a problem. Microlensing constraints on QCD PBHs should be re-investigated. We caution, however, in this numerically challenging problem some possibly relevant effects could not be treated.

The LIGO/Virgo collaboration have by now detected the mergers of ten black hole binaries via the emission of gravitational radiation. The hypothesis that these black holes have formed during the cosmic QCD epoch and make up all of the cosmic dark matter, has been rejected by many authors reasoning that, among other constraints, primordial black hole (PBH) dark matter would lead to orders of magnitude larger merger rates than observed. We revisit the calculation of the present PBH merger rate. Solar mass PBHs form clusters at fairly high redshifts, which evaporate at lower redshifts. We consider in detail the evolution of binary properties in such clusters due to three-body interactions between the two PBH binary members and a third by-passing PBH, for the first time, by full numerical integration. A Monte-Carlo analysis shows that formerly predicted merger rates are reduced by orders of magnitude due to such interactions. The natural prediction of PBH dark matter formed during the QCD epoch yields a pronounced peak around 1M⊙ with a small mass fraction of PBHs on a shoulder around 30M⊙, dictated by the well-determined equation of state during the QCD epoch. We employ this fact to make a tentative prediction of the merger rate of ∼ 30M⊙ PBH binaries, and find it very close to that determined by LIGO/Virgo. Furthermore we show that current LIGO/Virgo limits on the existence of ∼ M⊙ binaries do not exclude QCD PBHs to make up all of the cosmic dark matter. Neither do constraints on QCD PBHs from the stochastic gravitational background, pre-recombination accretion, or dwarf galaxies pose a problem.

Microlensing constraints on QCD PBHs should be re-investigated. We caution, however, in this numerically challenging problem some possibly relevant effects could not be treated.

Five years ago the decades-long efforts of the LIGO collaboration have taken fruit and led to the first unambiguous detection [1] of the gravitational wave signal of an inspiraling massive black hole binary. Since this pivotal moment a larger number of further binary coalescence have been observed partially in collaboration with Virgo [2] . Currently there are 66 detections of which 11 have been analysed. It came somewhat as a surprise that most of these binaries were made of fairly massive ∼ 30M ⊙ black holes (see however [3, 4] ). As star formation scenarios do not necessarily predict this, the theoretical community pondered the question if these black holes could be primordial and constitute the dark matter.

It has been known for long that mildly non-linear, horizon size overdensities could collapse and form an event horizon, i.e. a primordial black hole (PBH, hereafter) [5] [6] [7] (for a review cf. to [8] ). When such collapse occurs during radiation domination in the early Universe, the dynamics is characterized by a competition between self-gravity and pressure forces [9, 10] , and observes the physics of critical phenomena [11] [12] [13] . When pre-existing energy density perturbations, such as believed to emerge from inflationary scenarios, are featureless and almost scale-invariant, as observed in CMBR satellite missions, the equation of state during the PBH formation epoch plays a crucial role. It has been argued [9, 10] that PBH formation during the QCD epoch would be particularly efficient due to a softening of the equation of state. At the time of this realization, the QCD phase transition was believed to be of first order. Fully general relativistic numerical simulations of PBH formation confirmed that PBHs form more easily during the QCD epoch [14] , leading to a pronounced peak of PBHs on the ∼ 1M ⊙ scale. Though the simulations were performed under the assumption of a first order transition, it was argued in [9] , that any softening of the equation state, as for example during higher order transitions, or annihilation phases, such as the e + e − annihilation, would lead to a preferred scale in the PBH mass function. In the case of the latter this mass scale is around ∼ 10 5 − 10 6 M ⊙ . With the advancements of lattice gauge simulations it was possible to derive the zero chemical potential QCD equation of state with high precision [15, 16] . This equation of state, was recently used in approximate analytic calculations to derive the putative PBH mass function [17] [18] [19] . This mass function indeed has a very well developed peak at M ≈ 1M ⊙ and broader shoulder around M ∼ 30M ⊙ . It is tantalizing that this second scale is close to the inferred BH binary masses by the Ligo/Virgo collaboration.

There is a whole suite of observational constraints on the abundances of PBHs in the solar mass range, including the MACHO/Eros Milky Way microlensing constraints [20, 21] , constraints from the CMBR due to PBH accretion [22] [23] [24] [25] [26] , dynamical constraints, such as constraints from the abundances of wide stellar binaries in the galactic halo [27] [28] [29] , constraints from dwarf galaxies [30, 31] (cf. to [32] for a review). These constraints have been recently amended by limits on the PBH mass fraction contributing to dark matter, due to the observed BH binary coalescence rates form LIGO/Virgo data and from the non-observation of a gravitational stochastic wave background due to PBH mergers [33] [34] [35] [36] [37] [38] [39] [40] [41] [42] . The verdict of the community [40, [43] [44] [45] [46] [47] [48] is that solar mass PBHs can not constitute the bulk of the dark matter, mostly because they would greatly surpass the observed LIGO/Virgo rates.

In this paper, we show that this conclusion is with high likelihood incorrect, we believe that QCD formed PBHs may very well constitute all of the dark matter. The findings of this paper, should stimulate intense research on a numerically difficult problem. Limits imposed from the coalescing rate rely on the correctness of the PBH binary semi-major and eccentricity distribution. The initial distribution has been first calculated in [49, 50] subsequently re-calculated by several authors [43, 45, 46] , and its subsequent evolution has been deemed to be too unimportant [40, [44] [45] [46] [47] [48] to circumvent LIGO/Virgo constraints. We will show that this distribution, in fact, dramatically evolves between the first formation of binaries and the present day, allowing for consistency with current data.

PBH dark matter has the characteristics of perfect cold dark matter on large scales, however, there are two important differences when compared to particle cold dark matter. First, PBH dark matter is formed infinitely cold as any peculiar motions of density fluctuations should have inflated away during inflation. When binaries then form from such initial conditions, the distribution function P (a, e) of binary semi-major axis a and eccentricity e is highly peaked at e ≈ 1, reflecting very radial orbits. The binary coalescence time, due to the emission of gravitational waves is given by

and becomes very short when e approaches unity. Any perturbation of the binary would be expected to typically lower e, to get to a more equilibrated, unpeaked steady state in P (a, e). This in turn drastically increases t gw [46, 47, 51, 52] . We will show that simple three-body interactions between the two binary members and a third PBH indeed have this effect. The second difference to particle dark matter pertains to the granularity of PBH dark matter. As individual dark matter particles are fairly massive, PBH dark matter has an additional isocurvature density fluctuation component on small scales simply due to small-number statistics and the associated Poisson noise, i.e the fluctuations of the number of PBHs in a given volume containing N PBHs on average, is δ(N ) = ∆N/N = 1/ √ N . Those isocurvature fluctuations lead to the early formation of PBH clusters at redshifts as high as ∼ 1000. These clusters have been observed in numerical simulations, with the smallest clusters observed in [46] and a full range of clusters observed in [53] . Contrary to prior claims [54, 55] the PBH isocurvature fluctuations have been confirmed to be Poissonian [56, 57] .

Following [53] the growth factor of isocurvature perturbations is given by

where a is scale factor and eq denotes matter-radiation equality. Assuming a spherical top hat collapse model, we may require D(a)δ(N ) ≈ 1.68 to determine the collapse redshift a coll . In this model the final cluster density is given by ρ cl ≈ 178 ρ DM (a coll ) where ρ DM is the average dark matter density. These relations allow us to determine the PBH number density in clusters

and their virial velocity v cl ≈ 0.60 km s

as a function of the number of PBHs in the cluster N cl and PBH mass M pbh . Their properties are shown in Table 1 .

It is well known [58, 59] that dense clusters with a moderate number of cluster members N cl are unstable towards complete evaporation over the life-time of the Universe. Small clusters are therefore only of temporary existence between the initial collapse redshift z coll and the evaporation redshift z evap The approximate time scale for complete evaporation (up to a last remaining binary) of PBH clusters may be estimated [58] by

where t relax is the cluster relaxation time and t cross ≈ R cl /v cl is the cluster crossing time. Here R cl can be determined by the relation (4π/3)n cl R 3 cl = N cl . Both the formation-and evaporation-redshifts, z f orm and z evap , are also shown in Table 1 .

We have considered three-body interactions in these cluster environments.

where v cl th is the virial cluster velocity and v int b is the internal binary velocity, we may not use the impact approximation to derive analytic results, but rather have to resort to full numerical computations. A three-body scattering event is characterized by a variety of parameters, binary member masses M b1 , M b2 perturber mass M p , semi-major axis a, eccentricity e, two inclination angles of the binary plane θ and φ with respect to the perturber velocity direction, the position of the binary members on the orbit given by a parameter ψ, impact parameter at infinity b, as well as the perturber velocity V p in the binary CM frame. The accuracy of individual scattering calculations is confirmed by surveilling energy and angular momentum conservation. Fig.1 and Fig.2 show one example of such three-body scattering. This not atypical scattering event on a hard and very eccentric binary illustrates the complexity of such scatterings. Several times the binary changes partner, binary a and e are constantly modified, until ultimately one PBH escapes to infinity, leaving a binary which is even harder, in agreement with the Heggie-Hills conjecture [60, 61] . However, though the binary is harder, its eccentricity has changed away from e ≈ 1, implying that t gw has increased from an initial 4.48 × 10 −2 Gyr to t gw = 5.0 × 10 10 Gyr.

We have performed Monte-Carlo (hereafter MC) simulations of the effects of three-body scatterings on preexisting binary properties due to their presence in PBH clusters between redshifts z f orm and z evap . The initial binary properties were taken from the distribution calculated in [43, 49] . The MC simulation considers the 10 most important scattering events, where importance is gauged by the smallness of the impact parameter b. The simulation assumes N cl = 100, however, results for N cl = 10 and N cl = 1300 are very similar. For our simulations M b1 = M b2 = M p = M ⊙ are assumed. Cluster frame binary and perturber velocities are given random directions and their magnitude is drawn from a flat distribution ranging from v p = 0 to v p = v cl th for perturbers and v b = 0 and v b = v cl th / √ 2 for binaries. Here the factor of √ 2 takes into account that binaries are twice as heavy. The perturber velocity in the binary CM frame is subsequently obtained by a Galilean transformation. The angles θ, φ, and ψ are randomly drawn from their appropriate distribution functions. Results are shown in Figure 3 . Figure 3 illustrates the following effects. All hard binaries, i.e. binaries with a < ∼ 10 −3 pc are getting harder due to three-body interactions, consistent with the Heggie-Hills conjecture [60, 61] . Eccentricities are drastically lowered due to scatterings. For example, most of the binaries with κ i ∼ 10 −4 , corresponding to eccentricities e i ∼ 0.99995 are scattered to values κ f > ∼ 10 −2 , corresponding to e f < 0.995. Here κ is defined as κ = (1−e 2 ). By this effect t gw increases by a factor > 10 7 . It is seen from the lowest panel in Fig.3 that most binaries with t gw < ∼ 10 13 Gyr get scattered to t gw > ∼ 10 13 Gyr. This is likely due to an equilibration of binary properties due to the complexity of multiple binary-perturber interactions as the one shown in the Fig. 1 and Fig. 2. Fig. 3 illustrates that any conclusion of limits on PBH dark matter relying on the distributions computed in [43, 49] is highly suspect. Many authors [17, 35, 36, 38, 44, 45, 47, 48, 62] have used this distribution to draw conclusions.

Due to PBH three-body interactions, i.e. close encounters of binaries with single PBHs, many of the binaries do get destroyed. Table 1 lists the fraction of binaries which get destroyed in PBH clusters. It also lists the fraction of PBHs which stay in the cluster after destruction. We observe in our MC simulations that only the wider binaries get destroyed, with the hard binaries staying intact. It is also observed that a large fraction of the binaries stay in the cluster. We note here that our criteria for ejection from the cluster is that the post-scattering velocity is larger than the escape velocity, assumed to be v cl esc = 2v cl th [58] . The internal energy of a binary is given by E b = −GM b1 M b2 /2a. As a binary gets destroyed, the cluster kinetic energies of the PBHs are reduced by this amount. This reduction of total cluster kinetic energy will lead to the shrinking of the cluster. Hard binaries get harder during scatterings, i.e. a reduces. In this case, the kinetic energies of perturber and binary are substantially larger after the three-body interaction leading to an expansion of the cluster. It is likely that the first effect dominates, due to the larger number of wide binaries. The probable shrinking of clusters is not taken into account in our study.

We have computed the evolution of the probability per unit time that a binary merges at time t, i.e. dP t (t)/dt (cf. Eq.(13) of [43] ). The problem is numerically chal- [43, 49] and final distributions after processing between τ f orm and τevap in N cl = 10 (green), N cl = 100 (blue), and N cl = 1300 (orange) clusters are shown. The total t dPt(t)/dt after successive membership in N cl = 10, N cl = 100, and N cl = 1300 clusters is shown by the red line. It is seen that dPt(t)/dt at t = 14 Gyr is reduced by a factor ∼ 3 × 10 −5 with respect to its initial value. (cf. to the appendix concerning further details on the production of this figure).

lenging because one needs to obtain the statistics of tails, where already a single scattering is difficult to calculate due to the mismatch of time scales, τ b ≪ τ cl , where τ b is the binary orbital period and τ cl is the typical PBH cluster crossing time. The tail one is interested in are the very few evolutions where a binary gets scattered into, or stays at, the small coalescence time t 0 ≈ 14 Gyr, the current cosmological time. We were not able to solve the complete problem, due to numerical restrictions, but only computed the evolution of dP t (t)/dt for all PBH binaries having t i < 100 Gyr. Results are shown in Fig. 4 . Shown is the change of t dP t (t)/dt from its initial value due to evolution in N cl = 10, 100 and 1300 clusters, during the time interval τ f orm (N cl ) and τ evap (N cl ). Fig. 4 illustrates that PBH clusters of widely differing N cl have in magnitude similar effects on dP t (t)/dt. One may have expected less evolution in larger clusters as their densities are low. However, this effect is almost canceled by the longer evaporation time of large clusters. The number of collisons over a given time interval ∆t is ∆N sc = πb 2 v cl n cl ∆t (6) where b is impact parameter. Treating only the N sc most important events for each cluster, one may show with the help of the equations above that

is almost independent of N cl , with the rhs of Eq. 7 being 1.03,1.00, and 0.81 for N cl = 10, N cl = 100, and I . Shown as a function of the number N cl of PBHs in the cluster, assuming M pbh = 1M⊙, are the fraction of PBH binaries which are destroyed by three-body interactions f d , the fraction of binaries which stay binary and stay in the cluster f bcl ,the fraction of PBH binaries which stay binary and are ejected from the cluster f bej , the fraction of PBHs which after binary disruption stay in the cluster f dcl , the formation redshift of the cluster, the evaporation redshift of the cluster, the cluster density, the cluster virial velocity, and the cluster radius.

N cl = 1300, respectively. In this way binaries receive major perturbations in three different environments, over strongly increasing time intervals and in strongly decreasing densities. We note here that for the initial conditions of the N cl = 100, and N cl = 1300 simulations, strictly speaking we should not use the results of [43, 49] , but rather the distribution of a and e obtained from processing in the smaller cluster of the previous generation. Our simulations are currently not sophisticated enough to do so.

It is important to realize that in the spirit of hierarchical structure formation, N cl = 10 clusters will at lower redshifts be incorporated in N cl = 100 clusters, which in turn, at even lower redshift will become part of a N cl = 1300 cluster. Binaries which had been ejected or evaporated from a smaller cluster will quickly thermalize by three-body interactions in the larger cluster, and may not escape the potential well of the larger cluster, unless they get ejected again by another catastrophic collision. It is therefore that the total reduction of dP t (t)/dt is multiplicative, being the product of reduction rates in all three different size cluster environments. The total t dP t (t)/dt is shown by the red line, illustrating a total factor ∼ 3 × 10 −5 reduction at t 0 ≈ 14 Gyr. This reduction may be directly applied to the observed merger rates as those are given by M = n avg pbh dP t (t 0 )/dt, with n avg pbh the average present day cosmic PBH number density.

When taking a scenario of PBHs formed during the QCD epoch in its most natural form, i.e. without any pre-tuned peak of inflationary perturbations on a particular scale, but only taking into account effects of the softening of the equation of state [9, 10] , most PBHs form on the M ⊙ scale, whereas only a mass fraction of ∼ 10 −2 form on the larger ∼ 30M ⊙ scale [17] [18] [19] . This estimate assumes Gaussian statistics. In this case the relevant limits on PBH dark matter from LIGO/Virgo come from two LIGO/Virgo analysis, limits on M ⊙ PBH dark matter by the non-observation of solar mass PBH binaries [62] and the estimated merger rate of BH-BH binaries on larger mass scales [2, 63] . Using the main result of this paper given in Fig. 4 , we may compute the expected merger rates of pre-existing ∼ M ⊙ PBHs binaries to be M M⊙ ≈ 40 Gpc −3 yr −1 , to be compared to the LIGO/Virgo imposed upper limit of M < M⊙ ≈ 5.2 × 10 3 Gpc −3 yr −1 [62] . For detailed merger rates the reader is referred to our companion paper [64] . In fact, pre-existing binary properties get so efficiently modified by collisions in clus-ters, that the total merging rate is dominated by merging of single PBHs colliding at very low impact parameter [33, 35, 65] in present-day clusters. We derive a rate of M d M⊙ ≈ 640 Gpc −3 yr −1 for N cl = 1300 clusters and M d M⊙ ≈ 160 Gpc −3 yr −1 for N cl = 3000 clusters [64] , well below the LIGO/Virgo limit. Estimating M 30M⊙ is a bit more tricky, as not even the initial distribution is known yet, and the evolution of dP t (t)/dt again is rendered numerically unfeasible due to the small number 3 × 10 −4 of ∼ 30M ⊙ per solar mass PBH. Obviously the evolution of dP t (t)/dt of ∼ 30M ⊙ PBHs has to take into account the evolution of many more solar mass black holes. It would be expected that scatterings of 30M ⊙ PBHs on M ⊙ − M ⊙ or 30M ⊙ − M ⊙ binaries may typically lead to the exchange of the lighter black hole by the heavier one. This may well be the most important contribution to the abundance of 30M ⊙ − 30M ⊙ binaries. In any case, for a tentative estimate, we assume here that as with M ⊙ PBHs the contribution of pre-existing binaries is negligible. When considering the direct merger rate in present day clusters, one has to take into account that ∼ 30 M ⊙ sink to the center of the cluster. We estimate a merger rate of M d 30M⊙ ≈ 12−150 Gpc −3 yr −1 [64] to be compared to the by LIGO/Virgo inferred rate of 53.2 +58.5 −28.8 Gpc −3 yr −1 [63] at 90% confidence level.

These results are rather promising. However, we caution that we were not able to include several possibly relevant effects into the calculations. A question of utmost importance, is which fraction of PBHs never enter into a cluster. Leaning on results of CDM simulations [66] , this fraction could be rather small, but this has to be firmly established. In addition PBHs start infinitely cold which further could reduce this fraction compared to cold dark matter. Further important effects are, the neglected scattering of t i ≫ 10 Gyr binaries into the t f ≈ 10 Gyr band, which could dominate destruction, as there are many more (i.e. ∼ 100) t i > ∼ 10 Gyr binaries than t i < ∼ 10 Gyr, the shrinking of clusters due to binary destruction, and the role of core collapse of clusters [47] . Finally, the inclusion of heavy PBHs into the sea of light ones could also be of some importance. A complete treatment of this problem is currently not straightforward due to a hierarchy of time scales τ b ≪ τ cl ≪ τ H , all of which need to be resolved, where τ H is the Hubble scale. This may only be achieved, if at all, by elaborate adaptive mesh techniques.

It is somewhat unfortunate, that the stellar mass scale and that for QCD era produced PBH dark matter are essentially the same. This fact leads one to easily dismiss a potential discovery of dark matter, by attributing such observations to particularities of star formation. After all, we know stars exist. On the other hand, it is fortunate since gravitational wave detectors optimized to observe binaries on the stellar scale, become at the same time exquisite detectors for PBH dark matter formed during the QCD epoch. In case LIGO/Virgo (in combination with optical observations) detect only one binary including a sub-Chandrasekar mass black hole, M BH < ∼ 1.4M ⊙ , or in a weaker version, one binary including a BH of mass in the by star formation forbidden gap 2M ⊙ < ∼ M BH < ∼ 5M ⊙ , the existence of PBH dark matter may have been established.

In summary, the hypothesis that PBHs in the solar mass range make up the entirety of the dark matter had been seemingly ruled out by a number of observational constraints. However, taking our study at face value, even given the inherent uncertainties, the following statements can be made

• Current observational constraints from the merger rates of ∼ M ⊙ PBHs are evaded due to PBH collisions in clusters studied here.

• Similarly, observational constraints from the stochastic gravitational wave background due to PBH mergers should be evaded by the effect studied here.

• Observational constraints from the CMBR due to accretion or the dynamics of dwarf galaxies are evaded if only a small fraction of the dark matter is in ∼ 30M ⊙ PBHs, as naturally predicted by the QCD equation of state.

• The tentative estimate of the merger rate of ∼ 30M ⊙ PBHs is close to that observed by LIGO/Virgo.

• The mass scale of ∼ 30M ⊙ PBHs is predicted by the QCD equation of state.

Observational constraints from micro-lensing could still pose a problem, but they should be re-investigated in the context of PBH clusters. We believe that the possibility that LIGO/Virgo has detected the dark matter should be taken seriously, and demands further scrutiny.

Acknowledgments. I acknowledge detailed comments by Yacine Ali-Hamoud and Krzysztof Belczynski after first submission of this document. This paper was produced during Covid19 lockdown, allowing for uninterrupted work. I acknowledge my wife Nanci for being patient with me. For the production of Fig. 4 simulations with N cl = 10, 100 and 1300 have been performed. For each cluster size of the order of ∼ 5000 binaries were evolved. Here the N cl = 10 and 100 simulations consider the 10 most important scattering events, whereas the simulation of N cl = 1300 clusters, treat the N sc = 30 most important events. In the latter we had to switch to a larger N sc since otherwise spurious effects of Poisson statistics would arise. Given the long lifetime of such clusters, and the fact that we use Poisson statistics for the time when an event occurs, the probability to have not undergone any scattering before t is given by exp(− N sc t/t 0 ) which can be considerable at t ∼ 3 − 5 Gyr for low N sc , compared to the tail we want to compute. Each scattering event was limited to five minutes of CPU time, and abandoned if taking longer. Thus around 15% of the to-beevolved binaries were abandoned. These cases, typically very hard binaries, would subsequently be re-simulated with N sc = 4 and fifteen minutes of CPU time limit. Here the reduction in N sc helps as it reduces z, and since for small a/z simulations become very difficult. After this second round of integration 0.8%, 0.4%, and 0.05% for N cl = 10, 100, and 1300 clusters, respectively, had to be abandoned. These abandoned runs are included in the figure with their initial values. All of those abandoned runs were very hard binaries a ∼ 10 −5 pc at eccentricities very close to unity, with t gw substantially smaller than the present cosmic time.

Galactic Dynamics