Implications for primordial black holes from cosmological constraints on scalar-induced gravitational wave

Sufficiently large scalar perturbations in the early Universe can create over-dense regions that collapse into primordial black holes (PBH). This process is accompanied by the emission of scalar-induced gravitational waves (SIGW) that behave like an extra radiation component, thus contributing to the relativistic degrees of freedom ($N_{\rm{eff}}$). We show that the cosmological constraints on $N_{\rm{eff}}$ can be used to pose stringent limits on PBHs created from this particular scenario as well as the relevant small-scale curvature perturbation ($\mathcal{P}_{\mathcal{R}}(k)$). We show that the combination of cosmic microwave background (CMB), baryon acoustic oscillation (BAO) and Big-Bang nucleosynthesis (BBN) datasets can exclude supermassive PBHs with peak mass $M_{\bullet} \in [5 \times 10^{5}, 5 \times 10^{10}]\,{\rm M}_{\odot}$ as the major component of dark matter, while the detailed constraints depend on the shape of the PBHs mass distribution. The future CMB mission like CMB-S4 can broaden this constraint window to a much larger range $M_{\bullet} \in [8 \times 10^{-5}, 5 \times 10^{10}]\,{\rm M}_{\odot}$, covering sub-stellar masses. These limits on PBH correspond to a tightened constraint on $\mathcal{P}_{\mathcal{R}}$ on scales of $k \in [10, 10^{22}]\ {\rm{Mpc^{-1}}}$, much smaller than those probed by direct CMB and large-scale structure power spectra.


INTRODUCTION
Large density fluctuation in the early Universe can lead to gravitational collapse of over-dense regions and form primordial black holes (PBHs ;Hawking 1971;Carr & Hawking 1974;. PBHs have been proposed to explain a range of observed black hole conundrums (Carr et al. 2021a), such as the existence of supermassive black holes (Bean & Magueijo 2002) and black hole merging events seen by LIGO and Virgo (Abbott et al. 2019;Jedamzik 2021;Hütsi et al. 2021). PBHs are also considered as a dark matter (DM) candidate (Carr et al. 2016;Carr & Kuhnel 2020 and their abundance has been constrained via various cosmological and astrophysical observations (Khlopov 2010;Tashiro & Sugiyama 2013;Belotsky et al. 2014;Carr et al. 2016;Ali-Haïmoud & Kamionkowski 2017;Wang et al. 2018; 1 Mpc −1 , the curvature perturbation power spectrum P R has been precisely measured by cosmic microwave background (CMB) and large-scale structure (LSS) observations (Hunt & Sarkar 2015;Planck Collaboration et al. 2020), however the small-scale power is still poorly constrained (Inomata & Nakama 2019;Byrnes et al. 2019;Pi & Sasaki 2020;Chen et al. 2021;Yuan & Huang 2021;Kimura et al. 2021;Wang et al. 2022). Thus, the enhanced small scale P R can potentially produce enough PBHs to explain all dark matter without violating the existing observational bounds.
Associated with PBH production, enhanced P R inevitably generates scalar-induced gravitational wave (SIGW) at second order (Carbone & Matarrese 2005;Nakamura 2007;Yuan & Huang 2021;Chen et al. 2021;Karam et al. 2022), which can potentially be observed by gravitational wave (GW) detectors like Taiji (Hu & Wu 2017) and LISA (Amaro-Seoane et al. 2017). After horizon-crossing, SIGW can be considered as dark radiation (DR; Chacko et al. 2015;Serpico 2019;Takahashi & Yamada 2019) in that it free-streams with an energy density redshifting as (1 + z) 4 (Inomata & Nakama 2019;Chen et al. 2021). From a cosmological point of view, the gravitational behavior of SIGW is indistinguishable from a relativistic species such as the massless neutrino. Therefore, just as other forms of dark radiation (e.g. primordial gravitational wave (Meerburg et al. 2015;Aich et al. 2020), axion (Green et al. 2019) and sterile neutrinos (Takahashi & Yamada 2019; Green et al. 2019)), the energy density of SIGW DR can be parameterised by an additional effective degree of freedom where N eff describes the total effects of relativistic species containing both neutrino and SIGW, and N SM eff = 3.046 is the prediction from the standard model (SM) of particle physics (Mangano et al. 2002(Mangano et al. , 2005de Salas & Pastor 2016).
As a measure of cosmic radiation density, a higher N eff can delay radiation-to-matter equality and change the size of the sound horizon (Aich et al. 2020), which can leave distinctive features on CMB anisotropies (Hou et al. 2013;Follin et al. 2015), baryon acoustic oscillations (BAO) and Big-Bang nucleosynthesis (BBN) (Wallisch 2018). Currently the leading ∆N eff constraint is set by measuring the damping tail and the phase shift in the CMB anisotropy spectrum (Hou et al. 2013;Ade et al. 2016;Follin et al. 2015;Wallisch 2018;Aich et al. 2020), which reads ∆N eff < 0.3 at 95% confidence level (C.L.) from the latest Planck CMB data (Planck Collaboration et al. 2020). Future CMB missions are expected to give significant improvements on ∆N eff (Abazajian et al. 2016;Di Valentino et al. 2018;Hanany et al. 2019). For example, the CMB Stage IV (S4) experiment is aiming to constrain ∆N eff < 0.027 at 2σ C.L. (Baumann et al. 2016a;Abazajian et al. 2016;Baumann et al. 2016b;Wallisch 2018).
Here we use the SIGW upper bounds inferred from cosmological constraints on ∆N eff to place limits on both small-scale P R and the associated PBH abundance. The structure of this paper is as follows: Sec. 2 reviews PBH production induced by enhanced curvature perturbation. Sec. 3 discusses the energy density of SIGW and its connection with N eff . Our results are presented in Sec. 4 and we conclude in Sec. 5.

PBH MODEL
In order for inflation to produce density fluctuation required for efficient PBH formation, the primordial curvature perturbation power spectrum P R needs to be boosted to 10 −2 on small scales (k 1 Mpc −1 ) (Josan et al. 2009;Bringmann et al. 2012;Garcia-Bellido & Ruiz Morales 2017;Pi & Sasaki 2020;Green & Kavanagh 2021;Chen et al. 2021;Karam et al. 2022), which can be realised in several inflation theories (Kawasaki et al. 1998;Yokoyama 1998a;Kohri et al. 2013;Garcia-Bellido & Ruiz Morales 2017;Kannike et al. 2017;Cai et al. 2020;Pi & Sasaki 2020). As a good approximation for a wide range of curvature perturbations (Inomata & Nakama 2019;Pi & Sasaki 2020;Chen et al. 2021;Yuan & Huang 2021;Domènech 2021), we adopt a log-normal P R model peaked at the scale k • which was motivated in the Horndeski theory of gravity (Inomata & Nakama 2019;Chen et al. 2021) where A s = 2.1 × 10 −9 and n s = 0.9665 are the amplitude and spectral index of primordial fluctuations fixed at the pivot scale k * = 0.05 Mpc −1 (Planck Collaboration et al. 2020). Model parameters A, k • and σ describe the amplitude, location and width of the enhanced peak. For large scales (k k • ) or for A → 0, P R is reduced to a power-law primordial spectrum.
Once the large scalar fluctuations crossed the Hubble radius during radiation dominated era, the overdensity above a certain threshold (δ c 0.45, see e.g. Musco & Miller 2013;Harada et al. 2013;Carr & Kuhnel 2020 1 ) can gravitationally collapse into a PBH with mass (Carr et al. 2010;Nakama et al. 2017;Özsoy et al. 2018;Chen et al. 2021), where γ is the collapse efficiency, for which we adopt a typical value of γ = 0.2 Carr et al. 2010;Özsoy et al. 2018;Chen et al. 2021). g * (T ) is the total number of effectively massless degrees of freedom (those species with m T ) at horizon-crossing (k = aH), which we adopted SM particle content (Kolb & Turner 1990;Wallisch 2018)).
We use Press-Schechter model (Press & Schechter 1974) to calculate the fraction of Universe's density collapsing into PBHs with mass M Özsoy et al. 2018;Carr & Kuhnel 2020;Chen et al. 2021) where ρ • and ρ cr (z) are the PBH density and the critical density of the Universe at the collapsing time. In deriving the third equality in Eq. (4) we used δ c >σ, which remains valid for all scenarios we explored.σ 2 is the variance of density fluctuations at the scale k −1 (Young et al. 2014;Özsoy et al. 2018;Chen et al. 2021) where k(M ) function is defined via Eq. (3). W (x) is a window function, which we use exp(−x 2 /2) (Özsoy et al. 2018; Chen et al. 2021).
Since PBHs behave as matter, ρ • /ρ cr grows inversely proportional to temperature until matter-radiation equality, and the current distribution of PBH abundance in different masses can be estimated via (see also Young et al. 2014;Inomata et al. 2017;Özsoy et al. 2018;Chen et al. 2021 where f bh ≡ ρ • /ρ c is the fraction of cold dark matter in form of PBHs. For convenience, we describe the distribution calculated from Eq. (6) with a log-normal parameterisation (Dolgov & Silk 1993;Green 2016;Carr et al. 2017;Carr & Kuhnel 2020;Cang et al. 2022) which provides an excellent fit to the actual profile. Model parameters M • and σ • describe peak PBH mass and distribution width respectively, and their values are calculated using a least-square fitting method.

SIGW AS DARK RADIATION
In addition to forming PBHs, scalar perturbations also change the radiation quadruple moment and generate GW at second order, which carries an energy density ρ GW that redshifts as radiation after the horizoncrossing (Inomata & Nakama 2019;Chen et al. 2021), ρ GW ∝ (1+z) 4 . At the present day, the GW energy density parameter Ω GW ≡ (ρ GW /ρ cr ) z=0 follows the distribution (Ando et al. 2018 where Ω r = 9.1 × 10 −5 is present fractional radiation density assuming massless neutrinos. After neutrino decoupling, the cosmic radiation energy density (ρ r ) is a sum of CMB photon (γ), neutrino (ν) and GW densities where Here T γ = 2.728(1+z) K and T ν = (4/11) 1/3 T γ are temperatures of CMB and neutrino respectively. Because the behavior of GW density mimics that of neutrino, Eq. (11) counts the contribution of GW to the effective number of neutrino species as ∆N eff (Eq. (1)). Comparing ρ GW = Ω GW ρ cr (1 + z) 4 with Eq. (11), one can see that where θ indicates the model parameters. Depending on the choice of parameterisation, θ can be either the P R parameters (A, σ, k • ) defined in Eq. (2), or the PBH parameters (f bh , σ • , M • ) defined in Eq. (7). For a given set of P R parameters, the value of Ω GW is directly determined by integrating Eq. (8). To obtain the relation between Ω GW and PBH parameters, we constructed a three dimensional grid in P R parameter space and calculated Ω GW and (f bh , σ • , M • ) for each point on the grid, using the large discrete sample of Ω GW (f bh , σ • , M • ) function obtained in the process, we then built an interpolation function to calculate Ω GW for any (f bh , σ • , M • ) parameter values in between.
To derive our PBH and P R constraints, we use the N eff limits from Planck Collaboration et al.
(2020), which gives N eff = 3.04 ± 0.22 at 95% C.L. for Planck+BAO+BBN data (hereafter PBB) 2 . Using Gaussian statistics, this can be inverted to an 95% C.L. upper bound of ∆N eff < 0.175. We also study the prospective constraints from a future CMB Stage IV experiment, which is expected to constrain ∆N eff < 0.027 at 95% confidence level (Baumann et Through Eq. (12), Eq. (13) is inverted to the upper bounds on SIGW density at 95% C.L.
In the following, we will substitute Eq. (2) into Eq. (8) to calculate Ω GW or numerically using the Ω GW (f bh , σ • , M • ) interpolation function, and then use the upper bound in Eq. (14) to constrain the model parameters. We also notice that because of the positive 2 Note that there are several updated N eff constraints from more recent CMB datasets, e.g. Planck and ACTPol jointly constrains N eff = 2.74 ± 0.17 at 68% C.L. (Aiola et al. 2020), which is lower than standard model N eff = 3.046 by 1.8σ C.L. This might be due to unaccounted systematics in the data. The SPT-3G and Planck data gives N eff = 3.00 ± 0.18 at 68% C.L (Balkenhol et al. 2022), which has larger error than the PBB limits.
correlation between N eff and H 0 , adopting a different value of H 0 may lead to the shift of posterior distribution of N eff . In recent years, measurements from the local distance ladder (Riess et al. 2016(Riess et al. , 2018Riess et al. 2022) give a higher value of H 0 than the CMB (Planck Collaboration et al. 2020; Aiola et al. 2020;Balkenhol et al. 2022) and BAO measurements (Schöneberg et al. 2022), which are inconsistent with each other at almost 5σ C.L. (Verde et al. 2019). However, we notice that there are still a lot of debates on the potential systematics of Cepheids measurement (Freedman et al. 2020;Efstathiou 2020Efstathiou , 2021, which requires more data in the future to clarify the Cepheids measurement. The uncertainty of N eff associated with this on-going debates of anchor galaxy distance certainly exceeds the scope of this paper, but we notice the reader for this potential effect on N eff measurement.

RESULTS
For given σ and k, P R is determined by k • and A parameters. Once we set k • = k · exp[(1 − n s )σ 2 ], which gives a P R spectra that peaks at k (see Eq. (2)), the upper bound on P R (k) can be obtained by fixing A to its upper limit given by Eq. (14). To avoid overlapping with existed large scales constraints (Hunt & Sarkar 2015), we restrict the solution to the range of k > 10 Mpc −1 , which corresponds to M • < 5 × 10 10 M from Eq. (3).
As shown in Fig. 1, in the majority of parameter space we explored, our results show that PBB constrains P R to be ∼ O(10 −1 ). The constraint becomes more stringent for a wider spectral width σ, which is because the peak amplitude of P R is roughly inversely proportional to σ. For a sharp spectra with σ = 0.1, our constraint yields P R 1, whereas the limit tightens to P R 0.1 for σ = 2. In most cases, the linear part in Eq. (2) can be safely ignored, so that P R ∝ A and Ω GW ∝ A 2 , therefore our P R upper limit is proportional to the maximally allowed √ Ω GW . Therefore, compared to PBB, the CMB-S4 experiment will uniformly improve its P R constraint by 60 per cent on all scales.
As with the majority of PBH formation theories (Dolgov & Silk 1993;Yokoyama 1998b;Niemeyer & Jedamzik 1999;Green 2016;Carr et al. 2017;Bellomo et al. 2018;Pi et al. 2018), our PBHs follow an extended distribution and can be well described by the log-normal parameterisation in Eq. (7). Using the Ω GW (f bh , σ bh , M bh ) interpolation function described in previous section, it can be shown numerically that Ω GW increases with f bh , and here we derive our f bh upper bounds by iteratively solving Eq. (14). In Fig. 2 we show the f bh limits for PBHs with distribution widths of σ • = 1 and σ • = 2, along with the existing bounds on monochromatic PBH (assuming all PBHs having same mass) summarised in Carr & Kuhnel (2020) and Serpico et al. (2020). Note that in addition to the inflationinduced scalar perturbation, there are many other PBH formation scenarios (Carr et al. 2021b), e.g. collapse of cosmic strings (Hogan 1984;Hawking 1989;Polnarev & Zembowicz 1991;Hansen et al. 2000) and the bubble collisions (Hawking et al. 1982;Crawford & Schramm 1982;La & Steinhardt 1989), thus caution should be taken while interpreting Fig. 2. Our constraints applies to PBHs formed from enhanced scalar perturbation, whereas the rest of the limits constrain PBHs regardless of their formation mechanism.
For constraint from PBB with σ • = 1, the blue filled region in Fig. 2 shows that it excludes supermassive PBHs with M • ∈ [5 × 10 5 , 5 × 10 10 ] M as the dominant DM component (f bh < 1), spanning over 5 orders of magnitude. The limit covers the range set by X-ray binary (light blue; Inoue & Kusenko 2017) and LSS (brown; Carr & Kuhnel 2020). For σ • = 2, the PBB exclusion window expands to [10 2 , 5 × 10 10 ] M , which constitutes the widest PBH constraints up-to-date. The red and blue dashed lines show the sensitivity of projected CMB-S4 experiment on PBH abundance, which can exclude PBHs in [8 × 10 −5 , 5 × 10 10 ] M while setting the most stringent PBH constraints in a large mass range of [7 × 10 −3 , 5 × 10 10 ] M . Compared to other leading constraints in the supermassive PBH window around [3 × 10 4 , 5 × 10 10 ] M , the limit from projected S4 can be stronger by more than 10 orders of magnitude. Fig. 3 shows the complete constraints for a range of σ • values. For fixed f bh , we find that Ω GW increases with both σ • and M • , therefore our f bh limit is tightened as we increase either σ • or M • .

SUMMARY
Scalar-perturbation-induced PBH formation events emit SIGWs that behave like a relativistic species, so that the measured ∆N eff can yield stringent lim-its on the scalar power spectra P R on small scales 10 < k < 10 22 Mpc −1 , much smaller than the direct CMB power spectra measurement. Using a log-normal parameterisation for P R arised from Horndeski gravity theory (Chen et al. 2021) (also a good approximation for a wide class of perturbation theories), we show that Planck CMB data combined with BBN and BAO datasets (PBB) gives the currently most stringent P R constraints in k ∈ [10 8 , 10 22 ] Mpc −1 . For PBHs with a log-normal width of σ • = 1, PBB exclude supermassive PBHs with peak mass M • ∈ [5 × 10 5 , 5 × 10 10 ] M as the dominant DM component. Future CMB-S4 can probe M • ∈ [8 × 10 −5 , 5 × 10 10 ] M mass window and potentially improve P R limits of PBB by 60%. J.C. and Y.G. acknowledges support from the National Natural Science Foundation of China (12275278), the China Scholarship Council (CSC, No.202104910395) and partially from the Ministry of Science and Technology of China (2020YFC2201601). Y.Z.M. is supported by the National Research Foundation of South Africa under grant No. 120385 and No. 120378, NITheCS program "New Insights into Astrophysics and Cosmology with Theoretical Models confronting Observational Data", and National Natural Science Foundation of China with project 12047503. J.C. thanks for the hospitality at the Cosmology Group of SNS during this project.