Cosmic Conundra Explained by Thermal History and Primordial Black Holes

TL;DR

This paper proposes a unified primordial black hole (PBH) formation model linked to early Universe phase transitions, explaining dark matter, black hole mergers, and cosmic background correlations, with testable predictions for gravitational wave observations.

Contribution

It introduces a multi-modal PBH mass spectrum model that connects cosmic conundra to early Universe phase transitions and predicts specific gravitational wave signatures.

Findings

PBH mass spectrum peaks at specific masses including stellar and supermassive scales.

Recent gravitational wave events match the predicted mass gaps, supporting the primordial origin hypothesis.

The model links PBH formation to early Universe phase transitions, offering new probes of fundamental physics.

Abstract

A universal mechanism may be responsible for several unresolved cosmic conundra. The sudden drop in the pressure of relativistic matter at $W^{\pm}/Z^{0}$ decoupling, the quark--hadron transition and $e^{+}e^{-}$ annihilation enhances the probability of primordial black hole (PBH) formation in the early Universe. Assuming the amplitude of the primordial curvature fluctuations is approximately scale-invariant, this implies a multi-modal PBH mass spectrum with peaks at $10^{-6}$, 1, 30, and $10^{6}\,M_{\odot}$. This suggests a unified PBH scenario which naturally explains the dark matter and recent microlensing observations, the LIGO/Virgo black hole mergers, the correlations in the cosmic infrared and X-ray backgrounds, and the origin of the supermassive black holes in galactic nuclei at high redshift. A distinctive prediction of our model is that LIGO/Virgo should observe black hole mergers in the mass gaps between 2 and $5\,M_{\odot}$ (where no stellar remnants are expected) and above $65\,M_{\odot}$ (where pair-instability supernovae occur) and low-mass-ratios in between. Therefore the recent detection of events GW190425, GW190814 and GW190521 with these features is striking confirmation of our prediction and may indicate a primordial origin for the black holes. In this case, the exponential sensitivity of the PBH abundance to the equation of state would offer a unique probe of the QCD phase transition. The detection of PBHs would also offer a novel way to probe the existence of new particles or phase transitions with energy between $1\,{\rm MeV}$ and $10^{10}\,$GeV.

Tables

Table 1. Table 1: Upper limits (at 95 % percent 95 95\%  C.L.) on the integrated PBH fraction f PBH tot superscript subscript 𝑓 PBH tot f_{\rm PBH}^{\rm tot} from the various probes discussed in the text, calculated for a PBH mass function with n ~ s = 0.96 subscript ~ 𝑛 s 0.96 \tilde{n}_{\mathrm{s}}=0.96 . The third column lists the main assumptions and uncertainties underlying these limits. The PBHs can account for all the dark matter if one uses a realistic estimate of the EROS microlensing constraints (fourth line). Uncertainties in CMB limits for this mass function are discussed in the text.

Probe	$f_{\mathrm{PBH}}^{\mathrm{tot}}$	Assumptions	References$1^{{}^1}$
	$0.08$	isothermal dark matter halo profile, no clustering	EROS Tisserand et al. (2007), OGLE Wyrzykowski et al. (2009)
Microlensing	$0.39$	conservative isothermal profile, no clustering	Green Green (2016)
(LMC/SMC)	$0.16$	realistic profile, no clustering	Calcino et al. Calcino et al. (2018)
	$3.9$	Green Green (2016) model plus $90\%$ of PBHs in clusters $>{10}^3{M}_{\odot }$	appendices
Eridanus II	$2.3$	IMBH of $M=2000M_{\odot }$ at centre	Li et al. Li et al. (2017), Brandt Brandt (2016)
X-ray/radio PS in GC	4.5	accretion parameter $\lambda =0.02$, Maxwellian velocity distribution	Gaggero et al. Gaggero et al. (2017)
Halo wide binaries	$19$	fixed relative velocity, constant dark matter profile	Quinn et al. Quinn et al. (2009)
Planck	$0.17$	conservative case, steady-state accretion, no clustering	Ali-Haïmoud et al. Ali-Haïmoud and Kamionkowski (2017)

Equations

$$\beta(M)\approx{\rm erfc}\!\left[\frac{\delta_{\mathrm{c}}\big{(}w[T(M)]\big{)}}{\sqrt{2}\,\delta_{\rm rms}(M)}\right],$$

$$T\approx 200\sqrt{M_{\odot }/M}\mathrm{MeV}.$$

$${\rho }_{\mathrm{PBH}}(M)\equiv M^2\frac{\mathrm{d}n}{\mathrm{d}M},f_{\mathrm{PBH}}=\frac{{\rho }_{\mathrm{PBH}}(M)}{{\rho }_{\mathrm{CDM}}},$$

$$f_{\mathrm{PBH}}(M)\approx 2.4\beta (M)\sqrt{\frac{M_{\mathrm{eq}}}{M}},$$

$${\delta }_{\mathrm{rms}}(M)=\{\begin{array}{ll}\tilde{A}{\left(\frac{M}{M_{\odot }}\right)}^{(1-{\tilde{n}}_{\mathrm{s}})/4}& \text{(PBH scales)}\\ A{\left(\frac{M}{M_{\odot }}\right)}^{(1-n_{\mathrm{s}})/4}& \text{(CMB scales).}\end{array}$$

$$f_{\mathrm{PBH}}^{\max }={\left(\int \mathrm{d}\text{mspace}-1.0mu\ln M\frac{f_{\mathrm{PBH}}(M)}{f_p^{\mathrm{mon}}(M)}\right)}^{-1},$$

$$\frac{\mathrm{d}{n}_{\mathrm{h}}}{\mathrm{d}\ln M_{\mathrm{h}}}\approx \frac{{\rho }_{\mathrm{m}}}{\sqrt{\pi }{M}_{\mathrm{h}}}{\mathrm{e}}^{-M_{\mathrm{h}}/M_{\ast }},$$

$$\frac{\mathrm{d}{n}_{\mathrm{PBH}}}{\mathrm{d}\ln M}\approx \frac{{\rho }_{\mathrm{m}}{f}_{\mathrm{PBH}}}{M}.$$

$$t_{\mathrm{rel}}\approx \left(\frac{N_{\mathrm{cl}}}{\ln N_{\mathrm{cl}}}\right)\frac{v_{\mathrm{vir}}}{R_{\mathrm{cl}}},$$

$$v_{\mathrm{vir}}\approx \sqrt{\frac{G{M}_{\mathrm{cl}}}{2R_{\mathrm{cl}}}},$$

$$\alpha (\xi )=\frac{4GM(\xi )}{c^2\xi },$$

$$M(\xi )=2\pi {\int }_0^{\xi }\mathrm{d}{\xi }^{\prime }\Sigma ({\xi }^{\prime }){\xi }^{\prime }$$

$$\Sigma (\xi )=\int \mathrm{d}z{\rho }_{\mathrm{PBH}}(\xi ,z)$$

$$\alpha (R_{\mathrm{cl}},M_{\mathrm{cl}})\approx 2\times 10^{-13}\left(\frac{M_{\mathrm{cl}}}{M_{\odot }}\right)\left(\frac{\mathrm{pc}}{R_{\mathrm{cl}}}\right)\approx 2\times 10^{-10}.$$

$$R_{\mathrm{E}}=2\sqrt{G{m}_{\mathrm{PBH}}x(1-x)\frac{D_{\mathrm{cl}}}{c^2}}\sim 10^{-8}\mathrm{pc},$$

$$\frac{{\mathrm{d}}^2\tau }{\mathrm{d}{m}_1\mathrm{d}{m}_2}\propto f_{\mathrm{PBH}}(m_1)f_{\mathrm{PBH}}(m_2)\times \frac{(m_1{m}_2{)}^{2/7}}{(m_1+m_2{)}^{10/7}}.$$

$$R_{\det }=\frac{\sqrt{5}}{24}\frac{(G\mathcal{M}{c}^3{)}^{5/6}}{{\pi }^{2/3}}\times \frac{1}{2.26}{\left[{\int }_{f_{\min }}^{f_{\max }}\mathrm{d}f\frac{f^{-\alpha }}{S_h(f)}\right]}^{1/2},$$

$$f_{\mathrm{ISCO}}(m_1,m_2)=\frac{4400\mathrm{Hz}}{(m_1+m_2)}\times \frac{4m_1{m}_2}{(m_1+m_2{)}^2},$$

$$f_{\mathrm{merge}}(m_1,m_2)=2f_{\mathrm{ISCO}}\times \frac{4m_1{m}_2}{(m_1+m_2{)}^2}.$$

Full text

Dark Matter, Primordial Black Holes, Quark-hadron transition, EW transition, Primordial Nucleosynthesis

Introduction — Primordial black holes (PBHs) in the solar-mass range have attracted a lot of attention since the LIGO/Virgo detection of gravitational waves from coalescing black holes Abbott et al. (2016). The observed merger rate is compatible with what would be expected if PBHs constitute an appreciable fraction, and possibly all, of the cold dark matter (CDM). Moreover, the LIGO/Virgo observations seem to favour mergers with low effective spins, as expected for PBHs but hard to explain for black holes of stellar origin Wysocki et al. (2018). An extended PBH mass function with a peak in the range $1$–$10\,M_{\odot}$ could explain the LIGO/Virgo observations. Based on an argument related to gravitational lensing by PBH clusters, we show that the usual dark-matter constraints from the microlensing of stars, supernovae and quasars in this range can be evaded.

Given the revival of interest in PBHs, one must explain why they have the mass and density required for explaining the LIGO/Virgo events, and why these values are comparable to the mass and density of stars. One approach is to choose an inflationary scenario which produces a peak in the power spectrum of curvature fluctuations at the appropriate scale García-Bellido et al. (1996). The required amplitude of the inhomogeneities must be much larger than that observed on cosmological scales but not too large, so this requires fine-tuning of both the scale and amplitude.

An alternative approach is to assume the power spectrum is smooth (i.e. featureless) but that there is a sudden change in the plasma pressure at a particular cosmological epoch, allowing PBHs to form more easily then. Enhanced gravitational collapse occurs because the critical density fluctuation required for PBH formation ($\delta_{\mathrm{c}}$) decreases when the equation-of-state parameter ($w\equiv p/\rho c^{2}$) is reduced. Since the PBH collapse fraction depends exponentially on $\delta_{\mathrm{c}}$ for Gaussian fluctuations Carr (1975), this can have a strong effect on the fraction of CDM in PBHs. This is particularly important for the Quantum Chromodynamics (QCD) transition at ${\sim}10^{-5}\,\mathrm{s}$, lattice-gauge-theory calculations indicating that the sound-speed decreases by around 30% then Schmid et al. (1999); Laine and Schroder (2006); Bhattacharya et al. (2014); Laine and Meyer (2015); Borsanyi et al. (2016); Byrnes et al. (2018); Saikawa and Shirai (2018).

PBHs formed at the QCD transition would naturally have the Chandrasekhar mass ($1.4\,M_{\odot}$), this also characterising the mass of main-sequence stars, and a collapse fraction of order the cosmic baryon-to-photon ratio (${\sim}10^{-9}$) if PBHs provide most of the dark matter Carr et al. (2019). The latter feature is naturally explained if PBH formation generates a hot outgoing shower of relativistic particles because electroweak baryogenesis can occur very efficiently there and produce a local baryon-to-photon ratio of order unity García-Bellido et al. (2019).

In this paper we point out an interesting consequence of the above scenario, by extending it beyond the QCD scale. As the background temperature decreases from $100$ GeV to $1$ MeV, corresponding to the rest masses of the W and Z bosons, the proton, the pion and the electron, there are four periods at which the sound speed exhibits sudden dips. The proton dip is the biggest (${\sim}30\%$) but the others can also be significant ($5$–10%) because of the exponential dependence of the gravitational collapse probability on the critical curvature fluctuation. These dips produce distinctive features in the PBH mass function at four mass scales in the range $10^{-6}$–$10^{6}\,M_{\odot}$.

An important feature of this scenario is that it predicts the form of the PBH mass distribution very precisely. We show that for a nearly scale-invariant primordial power spectrum, the expected form not only satisfies all the current astrophysical and cosmological constraints, but also allows the PBHs to explain numerous observational conundra: (1) microlensing events towards the Galactic bulge generated by planet-mass objects with about 1% of the CDM density Niikura et al. (2019), well above most expectations for free-floating planets; (2) microlensing of quasars Mediavilla et al. (2017), including ones that are so misaligned with the lensing galaxy that the probability of lensing by a star is very low; (3) the unexpected high number of microlensing events towards the Galactic bulge by dark objects in the mass gap between $2$ and $5\,M_{\odot}$ Wyrzykowski and Mandel (2019), where stellar evolution models fail to form black holes Brown et al. (1999); (4) unexplained correlations in the source-subtracted X-ray and cosmic infrared background fluctuations Kashlinsky et al. (2005); (5) the non-observation of ultra-faint dwarf galaxies below the critical radius of dynamical heating by PBHs Clesse and García-Bellido (2018); (6) the masses, spins and coalescence rates for the black holes found by LIGO/Virgo Abbott et al. (2018), including two recent events with black holes which are probably in the mass gap; (7) the relationship between the mass of a galaxy and that of its central black hole.

Thermal History of the Universe — Reheating at the end of inflation fills the Universe with radiation. In the absence of extensions beyond the Standard Model (SM) of particle physics (eg. with right-handed neutrinos), the Universe remains dominated by relativistic particles with an energy density decreasing as the fourth power of the temperature as the Universe expands. The number of relativistic degrees of freedom remains constant ($g_{*}=106.75$) until around $200$ GeV, when the temperature of the Universe falls to the mass thresholds of SM particles.

As shown in Fig. 1 (upper panel), the first particle to become non-relativistic is the top quark at $T\simeq m_{\mathrm{t}}=172\,$GeV, followed by the Higgs boson at $125\,$GeV, and the $Z$ and $W$ bosons at $92$ and $81\,$GeV, respectively. These particles become non-relativistic at nearly the same time and this induces a significant drop in the number of relativistic degrees of freedom down to $g_{*}=86.75$. There are further changes at the $b$ and $c$ quark and $\tau$-lepton thresholds but these are too small to appear in Fig. 1. Thereafter $g_{*}$ remains approximately constant until the QCD transition at around $200\,$MeV, when protons and neutrons condense out of the free light quarks and gluons. The number of relativistic degrees of freedom then falls abruptly to $g_{*}=17.25$. A little later the pions become non-relativistic and then the muons, giving $g_{*}=10.75$. Thereafter $g_{*}$ remains constant until $e^{+}e^{-}$ annihilation and neutrino decoupling at around $1\,$MeV, when it drops to $g_{*}=3.36$.

Whenever the number of relativistic degrees of freedom suddenly drops, it changes the effective equation of state parameter $w$. As shown in Fig. 1 (lower panel), there are thus four periods in the thermal history of the Universe when $w$ decreases. After each of these, $w$ resumes its relativistic value of $1/3$ but each sudden drop modifies the probability of gravitational collapse of any large curvature fluctuations present at that time. We will see below how these changes in $w$ result in the production of PBHs with different masses and dark-matter fractions.

The above discussion is subject to some uncertainties. Firstly, a proper description of phase transitions requires highly perturbative numerical methods which are not yet fully understood; deviations in $w$ of a few percent are possible Laine and Meyer (2015); Borsanyi et al. (2016); Saikawa and Shirai (2018); Vovchenko et al. (2020). Secondly, before the onset of neutrino oscillations, significant non-zero lepton flavour asymmetries may be present Oldengott and Schwarz (2017), leading to significant changes in the background plasma pressure Middeldorf-Wygas et al. (2020); Vovchenko et al. (2020); Bodeker et al. (2020). In particular, the Universe can pass through a pion condensate phase for large lepton flavour asymmetry Vovchenko et al. (2020). This softens the equation of state more than indicated in Fig. 1 and the QCD transition may even become 1st order Middeldorf-Wygas et al. (2020). However, in this work we assume zero lepton flavour asymmetry and base our results on Ref. Borsanyi et al. (2016).

Primordial Black Hole Formation — There are a plethora of mechanisms for PBH formation. All of them require the generation of large overdensities, specified by the density contrast, $\delta\equiv\delta\rho/\rho$, usually assumed to be of inflationary origin. When overdensities re-enter the Hubble horizon, they collapse if they are larger than some threshold $\delta_{\mathrm{c}}$, which generally depends on the equation of state and density profile. However, there are other (non-inflationary) scenarios for PBH formation, where the inhomogeneities arise from first-order phase transitions, bubble collisions, and the collapse of cosmic strings, necklaces, domain walls or non-standard vacua. Full references can be found in Ref. Carr and Kühnel (2020).

The threshold $\delta_{\mathrm{c}}$ is a function of the equation-of-state parameter $w(T)$, which is shown in Fig. 1, so the thermal history of the Universe can induce pronounced features in the PBH mass function even for a uniform power spectrum. This is because, if the PBHs form from Gaussian inhomogeneities with root-mean-square amplitude $\delta_{\rm rms}$, then the fraction of horizon patches undergoing collapse to PBHs when the temperature of the Universe is $T$ should be Carr (1975)

[TABLE]

where ‘erfc’ is the complementary error function and the temperature is related to the PBH mass by

[TABLE]

This shows that $\beta(M)$ is exponentially sensitive to $w(M)$. Throughout this work, we use the numerical results for $\delta_{\mathrm{c}}$ from Ref. Musco and Miller (2013). We stress that there are theoretical uncertainties in the value of $\delta_{\mathrm{c}}$, associated (for example) with the density profile in the collapsing region Musco et al. (2020); Young et al. (2014); Yoo et al. (2018); Musco (2019); Kuhnel et al. (2016). However, any change in its value can be counterbalanced by an adjustment in $\delta_{\rm rms}$, such that the current PBH dark matter fraction is preserved.

We need an expression for the present fraction of the CDM in PBHs with mass $M$. However, this requires clarification since the density of PBHs with a precise mass is not defined for a continuous mass function. If the number density of PBHs in the mass range $(M,\,M+\mathrm{d}M)$ is $\mathrm{d}n$, then one can define the mass density and dark matter fraction of PBHs with mass ‘around’ $M$ (i.e. in the mass range $M$ to $2M$) by

[TABLE]

where $\rho_{\rm CDM}$ is the CDM density. Integrating this over $M$ gives the total dark matter fraction $f^{\rm tot}_{\rm PBH}$. The present CDM fraction for PBHs with mass around $M$ is then

[TABLE]

where $M_{\rm eq}=2.8\times 10^{17}\,M_{\odot}$ is the horizon mass at matter-radiation equality. The numerical factor is $2\hskip 1.42262pt(1+\Omega_{\mathrm{b}}/\Omega_{\rm CDM})$, with $\Omega_{\rm CDM}=0.245$ and $\Omega_{\mathrm{b}}=0.0456$ being the CDM and baryon density parameters Aghanim et al. (2018).

There are many inflationary models and these predict a variety of shapes for $\delta_{\rm rms}(M)$. Some of them, including single-field models like Higgs inflation Ezquiaga et al. (2018) or two-field models like hybrid inflation Clesse and García-Bellido (2015), produce an extended plateau or dome-like feature in the power spectrum. Instead of focussing on any specific scenario, we will assume a spectrum of the form

[TABLE]

Here the spectral index $n_{\mathrm{s}}$ and amplitude $A$ are taken to have their CMB values Aghanim et al. (2018), $n_{\mathrm{s}}=0.97$ and $A=4\times 10^{-5}$, respectively, while the corresponding small-scale quantities, $\tilde{n}_{\mathrm{s}}$ and $\tilde{A}$, are treated as free phenomenological parameters. Indeed, Eq. (5) can represent any spectrum with an additional broad peak or small-scale enhancement, such as might be generically produced by a second phase of slow-roll inflation. In order to get an integrated PBH abundance of $f_{\rm PBH}^{\rm tot}=1$, the small-scale amplitude has to be $\tilde{A}=0.1487$ if $\tilde{n}_{\mathrm{s}}=n_{\mathrm{s}}$. As discussed in Appendix A1, non-Gaussian and non-linear effects Kalaja et al. (2019) can impact the overall PBH abundance, but one can rescale $\tilde{A}$ to give $f_{\rm PBH}^{\rm tot}=1$ without significantly affecting the mass function.

The ratio of the PBH mass and the horizon mass at re-entry is denoted by $\gamma$ and we assume $\gamma=0.7$ as a benchmark value, following Ref. García-Bellido et al. (2019). The resulting mass function is represented in Fig. 2. It exhibits a dominant peak at $M\simeq 2\,M_{\odot}$ and three additional bumps at $10^{-5}\,M_{\odot}$, $30\,M_{\odot}$ and $10^{6}\,M_{\odot}$, corresponding to transitions in the number of relativistic degrees of freedom predicted by the known thermal history of the Universe.111We also indicate the values of $M$ associated with the three recent LIGO/Virgo events, although this data only became available after our prediction. Jedamizk Jedamzik (1997) first drew attention to the dips at the QCD and electron-annihilation epochs. Byrnes et al. Byrnes et al. (2018) have derived the PBH mass function associated with the QCD transition but this omits the smallest and largest mass bumps. We note that the prominence of the latter depends on the tilt of the power spectrum and our choice of tilt makes it more significant. Byrnes et al. also incorporate a more detailed study of the effects of criticality on the low-mass tail of the PBH mass function. Although the time dependence of the equation of state of the thermal plasma is crucial in determining the PBH mass function, the tail effect is not large; it just lowers and slightly broadens the main peak at around $1\,M_{\odot}$.

Constraints —  In this section, we discuss whether the PBH mass functions shown in Fig. 2, all of which assume $f_{\rm PBH}^{\rm tot}=1$, are compatible with the numerous observational constraints on $f_{\rm PBH}(M)$. There is an underproduction of light PBHs for $\tilde{n}_{\mathrm{s}}\leq 0.955$ and of heavy ones for $\tilde{n}_{\mathrm{s}}\geq 0.965$ but we claim the mass distribution for $\tilde{n}_{\mathrm{s}}\simeq 0.96$ can provide 100% of the dark matter without violating any current reliable constraints, despite some claims to the contrary.

In order of increasing mass, the PBH constraints come from the extragalactic $\gamma$-ray background, microlensing surveys, dynamical effects (such as the heating of ultra-faint dwarf galaxies and their stellar clusters and the disruption of wide binaries), ratio and X-ray point source counts, and CMB anisotropies generated by PBH accretion. The limits are summarised in Ref. Carr et al. (2016) and numerous other papers. Most of these constraints assume a monochromatic PBH mass function (i.e. one with width $\Delta M\sim M$). In the present scenario we predict an extended mass function and cannot simply compare this with the monochromatic constraints. In order to assess the situation, we adopt the approach advocated in Ref. Carr et al. (2017). Assuming that the mass distribution scales linearly with $f_{\rm PBH}^{\rm tot}$, each probe $p$ sets an upper limit

[TABLE]

where $f_{p}^{\rm mon}(M)$ is the limit from probe $p$ for a monochromatic function of mass $M$. We have calculated the value of $f_{\rm PBH}^{\rm max}$ associated with each probe for $\tilde{n}_{\mathrm{s}}=0.96$ but different astrophysical assumptions. These are shown in Table 1 and discussed in more detail in the appendices material.

It is sometimes argued that the EROS/MACHO microlensing limits exclude solar-mass PBHs with $f_{\rm PBH}^{\rm tot}=1$ but this is based on various contentious assumptions (monochromatic mass function, no clustering, isothermal halo profile). In Appendix A2, we show that these limits are evaded in our scenario. This is because the primordial power spectrum is enhanced on small scales, so the corresponding inhomogeneities virialise much more quickly than in the standard scenario, with the PBHs forming compact clusters. The CMB limits of Refs. Ali-Haïmoud and Kamionkowski (2017); Poulin et al. (2017) are still in tension with our model but only for M\mathrel{\raise 1.29167pt\hbox{>\kern-7.5pt\lower 4.30554pt\hbox{\sim}}}10^{3}\,M_{\odot} and the steady-state accretion assumption breaks down for such large masses. The highest mass peak in Fig. 2 also conflicts with the CMB $\mu$-distortion limit Byrnes et al. (2019) (not shown in Fig. 2) but one can avoid this by invoking large non-Gaussianity Nakama et al. (2018). This changes the tails of the fluctuation distribution and typically increases the probability of collapse without changing the shape of the mass function Ezquiaga and García-Bellido (2018) if the distribution above the critical threshold is exponentially suppressed. It is even possible that PBHs formed without an enhancement of the primordial power spectrum due to quantum diffusion Ezquiaga et al. (2020). One can also evade these limits if the transition in the primordial power spectrum indicated by Eq. (5) occurs at sufficiently small scales but one cannot then explain the relation between the galactic halo and central black hole masses.

Observational Conundra — Besides passing the current observational constraints on the form of the CDM, the PBH mass function with $\tilde{n}_{\mathrm{s}}\simeq 0.96$ predicted from the known thermal history of the Universe provides a unified explanation for several other puzzling conundra. We discuss these in order of increasing PBH mass. The status of some of the conundra is still unclear but we include all of them to convey the breadth of predictions.

1. Planetary-Mass Microlenses. Recently Niikura et al. have reported two interesting microlensing results. The first Niikura et al. (2017) comes from observations of M31 using the Subaru telescope, which include one possible detection and place strong constraints on PBHs in the mass range $10^{-10}$ to $10^{-6}\,M_{\odot}$. The constraints are roughly compatible with our model and even the single candidate could be. The second Niikura et al. (2019) uses data from the five-year OGLE survey of 2622 microlensing events in the Galactic bulge Mróz et al. (2017) and has revealed six ultra-short ones attributable to planetary-mass objects between $10^{-6}$ and $10^{-4}\,M_{\odot}$. These would contribute about $1\%$ of the CDM, which is more than expected for free-floating planets van Elteren et al. (2019). This corresponds to the first bump in our predicted PBH mass function and the abundance, when integrated over the mass range probed by OGLE, coincides with our best-fit model with $\tilde{n}_{\mathrm{s}}\simeq 0.96$.

2. Quasar Microlensing. Hawkins has claimed for many years that quasar microlensing data suggest the existence of PBH dark matter Hawkins (1993). He originally argued for Jupiter-mass PBHs but later increased the mass estimate to $0.4\,M_{\odot}$ Hawkins (2006). More recently, the detection of $24$ microlensed quasars Mediavilla et al. (2017) suggests that up to $25\%$ of galactic halos could be in PBHs with mass between $0.05$ and $0.45\,M_{\odot}$ (somewhat below our main peak). These events could also be explained by intervening stars, but in several cases the stellar region of the lensing galaxy is not aligned with the quasar, which suggests a population of subsolar halo objects with $f_{\rm PBH}>0.01$. Indeed, Hawkins has argued that the most plausible microlensers are PBHs, either in galactic halos or distributed along the lines of sight to the quasars Hawkins (2020). For a PBH mass function with $\tilde{n}_{\mathrm{s}}=0.96$, one expects $f_{\rm PBH}\simeq 0.07$ in this mass range. In principle, Ref. Mediavilla et al. (2017) excludes all the dark matter being in the main peak at $2\,M_{\odot}$ but that conclusion can be circumvented if the PBHs are in clusters and we argue in Appendix A2 that only $10\%$ of them should be uniformly distributed.

3. OGLE/GAIA Excess of Dark Lenses in the Galactic Bulge. OGLE has detected around $60$ long-duration microlensing events, of which around $20$ have GAIA parallax measurements which break the mass-distance degeneracy and imply that they are probably black holes Wyrzykowski and Mandel (2019). The event distribution from the posterior likelihood of their masses peaks between $0.8$ and $5\,M_{\odot}$, which overlaps the gap from $2$ to $5\,M_{\odot}$ in which black holes are not expected to form as the endpoint of stellar evolution Brown et al. (1999). Although most of the dark matter is not in this mass gap, this is consistent with the main peak in the PBH mass distribution if $0.6\lesssim\gamma\lesssim 1$.

4. Cosmic Infrared/X-ray Backgrounds. As shown by Kashlinsky Kashlinsky et al. (2005); Kashlinsky (2016), the spatial coherence of the X-ray and infrared source-subtracted backgrounds implies that black holes are required. Although these need not be primordial, the level of the infrared background suggests an overabundance of high-redshift halos and this could be explained if a significant fraction of the CDM comprises PBHs larger than a few solar masses, the Poisson fluctuations in their number density then growing all the way from matter-radiation equality. In these halos, a few stars form and emit infrared radiation, while PBHs emit X-rays due to accretion. It is challenging to find other scenarios that naturally produce such features. The required mass cannot be specified precisely. Ref. Kashlinsky (2016) focuses on $30\,M_{\odot}$ PBHs and the Poissonian power is smaller by an order of magnitude for $3\,M_{\odot}$ PBHs. However, this only reduces the mass of high-$z$ halos by a factor of a few, so the IR background can still be explained with our mass distribution.

5. Ultra-Faint Dwarf Galaxies (UFDGs). For the PBH mass distribution shown in Fig. 2, the critical radius below which CDM-dominated UFDGs would be dynamically unstable is $r_{\mathrm{c}}\sim 10$–$20$ parsecs (depending on the mass of a possible central black hole). The non-detection of galaxies smaller than this critical radius, despite their magnitude being above the detection limit, suggests compact halo objects in the solar-mass range. Moreover, rapid accretion in the densest PBH halos could explain the extreme UFDG mass-to-light ratios observed Clesse and García-Bellido (2018). Recent $N$-body simulations Boldrini et al. (2020) confirm that this mechanism works for PBHs of $25$–$100\,M_{\odot}$ providing they provide at least $1$% of the dark matter.

6. Mass, Spin and Merger Rates for LIGO/Virgo Black Holes. Most of the observed coalesced black holes have effective spins compatible with zero Abbott et al. (2019a). Although the statistical significance of this result is still low, this goes against a stellar binary origin but is a prediction of the PBH scenario. We assume that most of the binaries form at late times rather than primordially and the reasons for this (somewhat controversial assumption) are given in Appendix A3. In this case, the expected rate of PBH mergers is comparable to that observed if PBHs account for a significant fraction of the CDM. With our mass distribution, PBHs in the range $10$–$100\,M_{\odot}$ have $f_{\rm PBH}(M)\sim 0.03$ even if $f_{\rm PBH}^{\rm tot}=1$. Using a method discussed in Appendix A3, we have computed the likelihood distribution of merger events with PBHs of masses $m_{1}$ and $m_{2}$ for $\tilde{n}_{\mathrm{s}}=0.96$. The results are shown in Fig. 3, which shows that the LIGO/Virgo events are in the most likely region. This region might shift to slightly lower masses if one takes into account the difference between the source frame and detector frame mass for the most distant black hole mergers. With the expected number of O3 events, LIGO/Virgo should be able to detect mergers of PBHs larger than $65\,M_{\odot}$ or with a low mass ratio, $q\equiv m_{2}/m_{1}\lesssim 0.1$, these being distinctive predictions of our scenario. The expected PBH merger rate in the solar-mass range, after normalising to the observed rate of $50\,{\rm yr}^{-1}\,{\rm Gpc}^{-3}$ in the larger mass range, is $\tau\approx 10^{3}\,{\rm yr}^{-1}\,{\rm Gpc}^{-3}$ for PBHs between $1$ and $5\,M_{\odot}$, which is below the rate inferred for neutron-star mergers but within range of the current LIGO/Virgo runs. Our scenario could therefore be probed by searching for BH mergers in the $2$–$5\,M_{\odot}$ mass gap or below the Chandrasekhar mass. These could be distinguished from neutron-star mergers using the maximum chirp frequency or non-detection of electromagnetic counterparts. The LIGO/Virgo collaboration has announced the probable detection of two BH mergers (GW190425 and GW190814) with one component in the mass gap Abbott et al. (2020a, b). These populate regions 4 or 5 of Fig. 3 and are consistent with our model. The first event Abbott et al. (2020a) could be a merger of PBHs at the “proton” peak, given that no electromagnetic counterpart was observed and one component has a mass above that expected for a neutron star. The second event corresponds to the “pion” plateau and a recent paper by Jedamzik Jedamzik (2020) also supports this conclusion. Subsequently, LIGO/Virgo has announced the detection of a black hole merger with at least one component in the pair-instability mass gap (GW190521) Abbott et al. (2020c, d), corresponding to region 2 of Fig. 3.

7. IMBHs and SMBHs. Given our mass distribution, we have calculated the number of intermediate-mass and supermassive PBHs for each $10^{12}\,M_{\odot}$ halo. Interestingly, we obtain about one $10^{8}\,M_{\odot}$ PBH per halo and $10$ times as many smaller ones, possibly seeding the formation of a comparable number of dwarf satellites and faint CDM halos. Assuming a standard Press-Schechter halo mass function Press and Schechter (1974),

[TABLE]

where $\rho_{\mathrm{m}}$ is the mean cosmological matter density (both dark and baryonic) and $M_{*}\approx 10^{14}\,M_{\odot}$ is the cut-off halo mass. For a given $M_{\mathrm{h}}$, one can thus identify the corresponding PBH mass that has the same number density,

[TABLE]

This gives a relation $M_{\mathrm{h}}\approx M_{\rm PBH}/f_{\rm PBH}$, corresponding to roughly one IMBH/SMBH per halo of mass $10^{3}\,M_{\rm PBH}$ for our distribution, which is in agreement with observations. Furthermore, a mass distribution with $\tilde{n}_{\mathrm{s}}\approx 0.96$ reproduces the observed relation between the central black hole mass and halo mass Kruijssen and Lutzgendorf (2013), as shown in Fig. 4, but only if $f_{\rm PBH}^{\rm tot}\simeq 1$. A lower (larger) value of the spectral index would imply too many (few) IMBH/SMBHs. Accretion should increase the mass of heavier SMBHs somewhat, and this would make the case $\tilde{n}_{\mathrm{s}}\approx 0.96$ in closer agreement with observations. Although the initial peak is at $10^{6}\,M_{\odot}$, PBHs of this mass would inevitably grow as a result of accretion.

Conclusions — Various cosmic conundra are naturally explained by the PBH mass function expected from the known thermal history of the Universe if $f_{\rm PBH}^{\rm tot}=1$, i.e. if PBHs constitute all of the dark matter. The current LIGO/Virgo run should measure the mass function of coalescing black holes rather precisely and, remarkably, two recent events coincide with the “proton” peak at around $1\,M_{\odot}$, while a third corresponds to the “pion” plateau at around $50\,M_{\odot}$. This is indicated in Fig. 2 and was a prediction of our model. PBHs from the “W/Z” bump are too small to be seen by LIGO/Virgo but they would be detectable by microlensing effects and may indeed have already been found in OGLE data. PBHs from the “electron” bump are too large to be seen by LIGO/Virgo but may be detected by their dynamical effects and may explain the relation between the masses of IMBHs in dwarf spheroidals or SMBHs in galactic nuclei and the masses of the host halos.

It is intriguing that extrapolating the physics of elementary particles back to the early Universe not only resolves the mystery of the dark matter but also addresses so many other cosmic conundra. On the other hand, if firm evidence for PBHs were found, a broad variety of astronomical observations Kashlinsky et al. (2019) could search for extra features in their mass function. This would probe the existence of any new particles thermally coupled to the primordial plasma, independently of their coupling to SM particles, for masses from $1$ MeV to $10^{10}$ GeV (above which PBHs should have evaporated), far beyond the energies accessible by any future particle accelerator.

The exponential sensitivity of the PBH abundance to the equation of state also means that PBHs can be used to probe the characteristics of the cosmic phase transitions at which they form. This is particularly relevant to the detection of gravitational waves from PBHs. For example, we have seen that the LIGO/Virgo results may probe the QCD transition and the presence of lepton flavour asymmetries associated with a pion condensation phase. It is also possible that NANOGrav may have detected a stochastic gravitational wave background and several groups have argued that this could be a 2nd-order background associated with PBH formation Kohri and Terada (2020); De Luca et al. (2020); Vaskonen and Veermäe (2020); Domènech and Pi (2020). If this interpretation of the data were confirmed, this would qualify as another important observational conundrum but we not discuss it further here.

Appendix A Appendices

*A1. Sensitivity to Primordial Non-Gaussianity and Model Dependence — *Non-Gaussian (NG) effects will change the probability of PBH formation and thus their dark matter contribution Bullock and Primack (1997); Kalaja et al. (2019). There are many NG effects that modify the probability of collapse, from changes in the tail of the primordial density contrast distribution function to non-linear effects in the gravitational collapse at PBH formation. We emphasise that our model for the generation of curvature fluctuations on QCD scales is different from the multiple-field (curvaton) model presented in Refs. Carr et al. (2019); García-Bellido et al. (2019). Here we envision an inflation model with two slow-roll phases. For example, Critical Higgs Inflation Ezquiaga et al. (2018); García-Bellido and Ruiz Morales (2017) may induce fluctuations of order $10^{-5}$ on the CMB scale and $0.1$ on the PBH scale. Both can be generated by dynamics consistent with present values of SM parameters and this gives similar spectral tilts ($\tilde{n}_{\mathrm{s}}\sim 0.96$) on two very different scales. The assumption of near-scale-invariance can be relaxed to describe more complex formation mechanisms, but the thermal history will still imprint the PBH mass function in a similar way. These features are therefore universal and would apply for any PBH model.

Note that the $\mathcal{O}(1)$ fluctuations needed for PBH collapse are generic in inflation Hawking (1982). What is unusual is the small amplitude observed in the CMB, which requires some adjustment of parameters, although realised very naturally in Higgs inflation Ezquiaga et al. (2018). Moreover, the NG change in the tail of the PDF of curvature fluctuations due to deviations from slow-roll at the quasi-inflection point in critical Higgs inflation will give rise to an exponential amplification of the probability of collapse to PBHs when those scales re-enter the Hubble scale during the radiation era Ezquiaga and García-Bellido (2018). We have therefore assumed a nearly scale-invariant spectrum characterises this type of single-field model of inflation. A more detailed study will be required to test those scales against the observed mass distribution in the LIGO/Virgo events.

*A2. Microlensing Limits — * The most important uncertainty in the microlensing limits comes from PBH clustering. Rather than being uniformly distributed in the galactic halo, as assumed by the MACHO, EROS and OGLE analyses, PBHs must be clustered to some degree. This applies for all forms of dark matter but the clustering in the PBH case is strengthened for two reasons. First, the Poisson fluctuations due to the discrete nature of PBHs give a dominant contribution to the matter power spectrum below kiloparsec scales Carr and Silk (2018). Our fourth cosmic conundrum relies on this effect. Second, the strongly enhanced and nearly scale-invariant primordial power spectrum ensures that any small-scale inhomogeneity is rapidly driven into the non-linear regime. Even if the impact on the halo mass function is not so significant, overdense regions collapse at much earlier times and this leads to the formation of dark matter clusters (minihalos) which are initially more compact than in the standard cosmological scenario.

The size of these PBH minihalos today is nevertheless limited by dynamical effects, just as for globular clusters. The typical relaxation time for a halo made of compact objects is

[TABLE]

where $N_{\rm cl}$, $R_{\rm cl}$ and $v_{\rm vir}$ are the number of PBHs in the minihalo, its radius and its virial velocity,

[TABLE]

respectively. Assuming one requires $t_{\rm rel}>10^{10}\,$yr for the clusters to survive and a mass $M_{\rm PBH}\approx 2\,M_{\odot}$, corresponding to the proton peak in the PBH mass function, one gets a lower limit $R_{\rm cl}\approx 10\,$pc. The derivation of this limit can be found in Ref. Clesse and García-Bellido (2018). The dependence on the halo mass is only logarithmic and so this critical scale is universal.

A more refined analysis of the dynamics of UFDGs, including the effects of a central IMBH, leads to a similar lower limit on $R_{\rm cl}$ and is supported by UFDG observations (see the fifth cosmic conundrum). For globular clusters, the dynamical hardening of binaries in the core can reduce the critical radius to parsec scales and this might also apply for some PBH clusters. But it is unlikely that smaller PBH clusters exist today, because they would be dynamically unstable.

Given that clustering is unavoidable, the main contribution to a uniform PBH distribution in our own galactic halo would come from PBHs which have been removed from their host cluster by close encounters and tidal interactions (i.e. the slingshot mechanism). Preliminary $N$-body simulations show that the fraction of objects ejected over the age of the Universe is about $10$%. The case of PBHs is not much different and those with a mass comparable to stars should be expelled at a similar rate. For subsolar PBHs, the uniform fraction could be enhanced but they are anyway contributing no more than a few percents of the dark matter, much below the microlensing limit on non-clustered PBHs.

We now explain why PBHs in clusters do not induce detectable microlensing events, thereby evading the usual LMC and SMC microlensing limits, even if PBHs provide all the dark matter. If a PBH cluster, assumed to be spherical for simplicity, is aligned with a star in one of the Magellanic clouds, it will induce strong gravitational lensing of this star, with a deflection angle

[TABLE]

where

[TABLE]

is the mass within a cylinder of radius $\xi$ with axis along the $z$ direction between the observer and the star, and

[TABLE]

is the projected surface density of the PBH cluster. For a cluster of mass $M_{\rm cl}$ and radius $R_{\rm cl}\approx 10\,$pc and $\xi\simeq R_{\rm cl}$, the deflection angle is

[TABLE]

Typically, this is not resolvable but the observed luminosity flux of the star is spread over an Einstein arc. At a cluster distance $D_{\rm cl}\sim\mathcal{O}(10\,{\rm kpc})$, the deflection angle is subtended by a length $L\sim D_{\rm cl}\hskip 1.42262pt\alpha\sim 10^{-9}\hskip 1.42262pt(M_{\rm cl}/\,M_{\odot})$, which for a cluster mass $M_{\rm cl}\gtrsim 10^{3}\,M_{\odot}$ is much larger than the Einstein radius $R_{\rm E}$ of an individual solar-mass PBH in the cluster:

[TABLE]

where $x\equiv D_{\rm cl}/D_{\rm LMC/SMC}$ is the distance ratio between the lens and the source. Therefore the star’s luminosity is only marginally affected by magnification due to a nearby PBH. Indeed, the microlensing by a PBH in a cluster would induce less than $10$% magnification of the star, whereas microlensing surveys only searched for $\mathcal{O}(1)$ magnifications Fleury and García-Bellido (2019). Only if the star is exactly aligned with the center of mass of the cluster, which is very unlikely, or if the cluster mass is below $10^{3}\,M_{\odot}$, which should only represent only a small fraction of the PBH density, would the magnification of the star’s luminosity be appreciable.

This new and rather simple argument provide robust motivation for the microlensing limits presented in line 4 of Table 1 and in Fig. 2. It also applies to the supernova Zumalacarregui and Seljak (2018); García-Bellido et al. (2018) and quasar microlensing limits, since the Einstein radius scales as $R_{\rm E}\propto D_{\rm cl}^{1/2}$ and the arc length as $L\propto D_{\rm cl}$, suppressing the deflection angle even more. Finally, invoking clustering to evade galactic and quasar microlensing limits does not contradict the PBH interpretation of OGLE microlensing events towards the Galactic center, because dense PBH clusters are probably tidally disrupted there, leading to a smoother PBH distribution. Nor does it affect our explanation of the second conundrum (quasar microlensing) since this only requires a homogeneous distribution with $f_{\rm PBH}\sim 0.01$–$0.1$, which is consistent with our clustering scenario.

*A3. Distribution of PBH Merger Detections — *The expected distribution of PBH merger detections as a function of the component masses $m_{1}$ and $m_{2}$ for LIGO/Virgo has been estimated for PBH binaries formed through tidal capture in dark matter halos. Their merging time $\tau$ is given by Clesse and García-Bellido (2017)

[TABLE]

The normalisation is unimportant as long as one is interested in the distribution of BH mergers. The total rate has uncertainties related to the halo mass function and concentration, the minihalo profiles, the PBH velocity distribution etc. Nevertheless, for realistic assumptions one can obtain a merging rate compatible with that inferred by LIGO/Virgo if PBHs constitute an appreciable fraction of the dark matter Bird et al. (2016); Clesse and García-Bellido (2017).

The distribution of mergers as a function of the component masses is less impacted by these astrophysical uncertainties. It is therefore a better discriminant between BHs of primordial and stellar origin. Knowing the detector range as a function of the black hole or neutron star binary chirp mass Chen et al. (2017), one can relate the rate distribution to the probability distribution for merger events. The latter is shown in Fig. 2 in the main text for the O2 observing run, assuming $n_{\mathrm{s}}=0.96$ and a PBH to horizon mass ratio $\gamma=0.7$.

The expected distribution of merger detections is the product of the comoving merging rate and the comoving volume probed by the detectors, $V_{\rm det}=(4\pi/3)\hskip 1.42262ptR_{\rm det}^{3}$, where $R_{\rm det}$ is the luminosity distance of the furthest detectable source with a signal-to-noise ratio exceeding the threshold value of $8$. This ‘astrophysical range’ depends on the detector sensitivity at the different frequencies of the gravitational wave (GW) signal. A simple estimate is given in Ref. Chen et al. (2017) (and references therein), where different effects are analysed with more precise numerical computations. Here we use the simple estimate,

[TABLE]

where $\mathcal{M}\equiv(m_{1}\hskip 1.42262ptm_{2})^{3/5}/(m_{1}+m_{2})^{1/5}$ is the chirp mass, the factor $2.26$ comes from the ratio of the range to the horizon in a Euclidian Universe, and the integral is the contribution to the signal-to-noise ratio for the whole signal duration over the frequency bandwidth of the detector, assumed to have a noise power spectral density $S_{h}$.

For a simple post-Newtonian waveform of the inspiral phase, one has $\alpha=7/3$ until the GW frequency associated with the innermost stable circular orbit (ISCO),

[TABLE]

where the black hole masses are in solar units. However, in our scenario numerous stellar-mass black holes from the peak of the mass function can be captured by heavier black holes from the pion bump thereby producing binaries with mass ratios of order 0.1, so it is crucial to take this effect into account. Typically, it reduces the astrophysical range and hence detection rate for low-mass ratios but there is still a significant probability of detection in the interesting region 5 of Fig. 5 where GW190814 has been observed. We also include the merging phase in our analysis, with $\alpha=2/3$ between $f_{\rm ISCO}$ and the merger frequency

[TABLE]

The last factor can be found in Ref. Lousto and Healy (2020) and is around one for mass ratios close to unity. We take the lower and upper bounds of the integral to be $f_{\rm min}=50\,{\rm Hz}$ and $f_{\rm max}=\min(f_{\rm merge},2000\,{\rm Hz})$. Finally, we account for redshift effects on the comoving rates and waveform by replacing $\mathcal{M}$ by $\mathcal{M}(1+z)$ in Eq. (A9). We have used the detector noise power spectral density $S_{h}$ adopted by the LIGO-Livingstone detector Abbott et al. (2017) in the last month of its second observing run.

This approach still has some inaccuracies. The detector range is wrong by a factor of a few, so this has been corrected by a linear rescaling in order for the range to be in reasonable agreement with Fig. 4 of Ref. Chen et al. (2017) and a binary neutron star range of $90$ Mpc Abbott et al. (2019b). In analysing the third observing run, a more accurate approach may be adopted by carrying out campaigns over the whole run with different sky locations and by using more accurate GW waveforms during the merging and ring-down phases. However, such an analysis is beyond the scope of this paper, in which we estimate the expected event distribution for the second LIGO/Virgo observing, using public data of the detector noise spectral density. The comparison with some events observed in the third observing run is thus only qualitative.

So far the expected event distribution is compatible with observation and its specific features are discussed in the main text. We do not consider the expected event distribution of detections for primordial binaries because it is still unclear how the expected PBH mass function would suppress the merging rate due to early-forming PBH halos seeded by Poisson fluctuations. A more refined analysis of this binary formation channel is needed and this will require $N$-body simulations.

A4. Late versus early binary formation — It is still unclear which binary-formation channel dominates. However, the primordial channel is subject to more uncertainties and there are two reasons for assuming this is disfavoured. First, we are considering a scenario with $f_{\rm PBH}=1$, in which case $N$-body simulations show that rapidly forming clusters lead to a suppression of the merging rate of primordial binaries for both a monochromatic and lognormal mass distribution. The importance of this suppression is still unclear but Ref. Vaskonen and Veermäe (2019) claims it could be up to two orders of magnitude. In our scenario, PBHs in the range $10$–$50\,M_{\odot}$ relevant for LIGO/Virgo account for only a few percent of the dark matter. The two effects combine and the merging rate becomes of order $1$–$10\,$yr-1 Gpc*-3*, consistent with the LIGO/Virgo lower bound. Consequently, the merging rate from captures in halos can dominate. Since the mass dependence of the suppression for primordial binaries is unknown for our mass function, we have only considered the capture channel because the merger rate distribution is subject to fewer uncertainties in this case.

Second, the rate limits from the O1– O2 runs of LIGO/Virgo already disfavour primordial binaries as the dominant channel for our mass distribution if we assume (for simplicity) a suppression rate independent of the mass or even no suppression at all. The merging rates should go roughly as $f(m_{\rm PBH})^{2}/m_{\rm PBH}$ for equal-mass binaries. For unequal-mass binaries, the rates calculated by Kocsis et al. Kocsis et al. (2018) and assumed by Gow et al. Gow et al. (2020) can be used but one comes to the same conclusion: with our bumpy mass function the limits set by LIGO/Virgo on the subsolar scale cannot be evaded for primordial binaries if they also explain some of the larger black hole mergers. We have also reproduced Fig. 2 for the case of primordial binaries in Fig. 5. The expected distribution of merger events in the ($m_{1},\hskip 1.42262ptm_{2}$) plane has a dominant peak at the solar-mass scale and only a subdominant number of events the high-mass region, which is the opposite of what LIGO/Virgo has observed. Finally, we point out the recent claim of Ref. Boehm et al. (2020) that the expected merging rate of primordial binaries is highly suppressed due to the time dependence of the Misner-Sharp mass in the Thakurta metric for the radial density flows.