Exploring Evaporating Primordial Black Holes with Gravitational Waves

Primordial black holes (PBHs) from the early Universe have been connected with the nature of dark matter and can significantly affect cosmological history. We show that coincidence dark radiation and density fluctuation gravitational wave signatures associated with evaporation of $\lesssim10^9$ g PBHs can be used to explore and obtain important hints about the formation mechanisms of spinning and non-spinning PBHs spanning orders of magnitude in mass-range, which is challenging to do otherwise.

Primordial black holes (PBHs) that could have formed in the early Universe prior to galaxies and stars can constitute a significant fraction of the dark matter (DM), can significantly affect cosmological history and have been associated with a variety of observational signatures (e.g. ). Depending on formation scenario, PBHs can span many orders of magnitude in mass.
Very light PBHs with masses 10 −15 g evaporate through Hawking radiation [42] on time scales shorter than the age of the Universe and thus do not contribute themselves directly to the DM abundance [43,44]. Evaporating PBHs can copiously emit constituents from the whole particle spectrum, including production of particles from the Standard Model extensions and the dark sector. Products of evaporating PBHs have been associated with the observed matter-antimatter asymmetry, DM and dark radiation (DR) [9,[45][46][47][48][49][50][51][52]. PBHs spanning range ∼ 10 −5 − 10 9 g, from the Planck mass to the mass relevant for Big Bang Nucleosynthesis (BBN), remain unconstrained (see e.g. [53]).
PBHs have been often considered as non-rotating (Schwarzschild) [54][55][56]. However, PBHs can be formed with significant spin (Kerr BHs) [24-26, 29, 57-60] as well as acquire spin via accretion [61] or hierarchical mergers [62]. Aside mass and charge, spin constitutes a fundamental conserved BH parameter. BH spin affects Hawking radiation as well as the BH lifetime, with potential consequences for observations [39,[63][64][65][66][67][68]. BH rotation generally increases the emission and favors particle production with larger spin [42,44,69,70]. As shown in Ref. [50] and confirmed by detailed numerical analyses [71,72], spin-2 graviton emission is significantly enhanced in rotating evaporating BHs, which dominates the resulting copious production of DR when PBHs of 10 9 g mass dominate the energy density and that will affect BBN and cosmic microwave background radiation (CMB) observations. Note that the emitted gravi- * domenech@pd.infn.it † volodymyr.takhistov@ipmu.jp ‡ misao.sasaki@ipmu.jp tons correspond to extremely high frequency gravitational waves (GWs), which are undetectable with current or near-future interferometer experiments (e.g. [64,73]). In this work we show that coincidence observations of DR, associated primarily with PBH Hawking radiation graviton production, as well as secondary GWs from PBH density fluctuations allow us to explore and distinguish formation mechanisms of spinning and nonspinning evaporating PBHs, opening an unprecedented window into the unconstrained 10 9 g PBH mass range.
On small scales, the inhomogeneous distribution of PBHs leads to density fluctuations [74]. As light PBHs dominate the Universe, the initial isocurvature fluctuations are converted into curvature perturbations. During the era of PBH domination, density fluctuations grow until PBHs evaporate completely. The PBH evaporation process, which is effectively instantaneous (i.e. sudden reheating approximation), transforms large density fluctuations into radiation and yields large pressure waves [75]. This induces a significant quadrupole moment and, therefore, results in significant GW production [76] (see also [77]).
We now discuss stochastic background GW production from PBH evaporation, including effects of spin, in generality as applicable to a large class of models. Consider a generic (dimensionless) power spectrum in wavenumber k of fluctuations in the gravitational potential Φ during a matter-dominated stage described by where A Φ is the amplitude, n is the spectral index and k UV is an ultraviolet cut-off. This type of spectrum includes the typical nearly scale-invariant spectrum (n = 0) found in models of inflation. As we show later, n ∼ 1 corresponds to PBH density fluctuations.
[Right] GW spectral density induced by early isocurvature fluctuations in the PBH reheating scenario in terms of frequency. The blue and green solid lines correspond to (M BH = 10 6 g, P Φ = 2 × 10 −22 ) and (M BH = 10 8 g, P Φ = 2 × 10 −25 ), respectively. GWs from Kerr PBHs are shown with dashed lines. Note that for this choice of parameters we respectively have keq/k reh ∼ 270, 1000 and, therefore, the Universe reaches complete PBH domination. We also display the power-law integrated sensitivity curves for LISA, DECIGO, ET and Advanced LIGO experiments [80]. In red we show the latest upper bounds on the stochastic GW background from the LIGO/VIRGO/KAGRA collaboration [81].
by multiplying with the appropriate kernel and integrating over the internal momenta. In the approximation of nearly sudden reheating, the contribution to the induced GWs at the peak can be schematically estimated as where q is the internal scalar loop momentum and S Φ (k/k reh ) is a suppression factor to account for the finite width of the transition, so that S Φ → 1 in the instantaneous transition limit. We refer for computation details to Ref. [75,76]. In the last step, we first used the fact that in the sudden reheating limit (associated with k reh ), the conformal time derivative of Φ, Φ , jumps from Φ = 0 in a matter-dominated stage to Φ = 0 in radiation-dominated stage. This yields a significant amplitude of Φ right after reheating and the kernel is essentially given by , where H is the conformal Hubble parameter. The additional factor (k reh /k UV ) originates from the fact that integral's contribution is dominated by a very narrow region, the width of which is set by this ratio. Here, we have also employed that fluctuations on scales much smaller than the typical transition time scale are affected by the finite width of the transition [75]. Due to radiation pressure, such small scale fluctuations are not efficiently produced which results in a suppression of Φ with respect to the exact instantaneous transition parametrised by S Φ (k/k reh ) [75]. We provide an explicit expression for S Φ later. We also assumed that the highest contribution to the integral comes from the high momentum limit. This implies that the power spectrum must not fall too fast at small scales, which leads to consider only spectral indices n < 2. With these assumptions and including the numerical coefficients, the amplitude of the peak GW emission frequency is and so the contribution to stochastic GW background is The subscript "reh" refers to the value of Ω GW right after reheating in the later radiation-dominated stage, where the energy density ratio of GWs and radiation stays constant. Note that by momentum conservation the GW spectrum has the cut-off at k ∼ 2k UV . However, the peak is at k ∼ k UV and then decays rapidly [76,77] (see Fig. 1). Thus, for simplicity, we cut the GW spectrum at k ∼ k UV . We now consider induced GW production in the case of evaporating light PBHs. For simplicity, we assume that PBHs form at a particular time, with a monochromatic mass spectrum, and are randomly distributed across the universe with a given initial number density n PBH,i . This implies that the initial distribution of PBHs follows Pois-son statistics and their (dimensionless) spectrum of fluctuations on small coarse-grained scales is [74] where the small scale cut-off is given in terms of the mean (comoving) distance between PBHs, that is where a i is the scale factor at the initial time. Then the gravitational potential power spectrum of fluctuations can be estimated as where T Φ (k) is a transfer function that connects the PBH formation stage, set in an arbitrary cosmological scenario, to the PBH domination regime. Let us consider the simplest scenario of PBH formation during an earlier radiation-dominated epoch. Then, the transfer function reads as [88,89] T where k eq = a eq H eq is the comoving horizon scale corresponding to radiation-matter (PBH) equality. Note that scales smaller that k eq are suppressed due to the fact that a mode k of Φ decays once it is inside the horizon. Some time after equality, the power spectrum of Eq. (7) is given by In the scenario where PBHs form in a matter-dominated era previous to the radiation-dominated era, there will be an additional scale k eq,eRD that corresponds to the comoving horizon scale when the universe enters the early radiation dominated stage, and an intermediate step needs to be considered in the transfer function above. However, the transfer function remains the same on sufficiently small scales k k eq,eRD , which are the scales of our interest.
In a PBH-dominated stage, the reheating time and duration of the transition are set by BH evaporation. It then follows that (for details see Appendix A) Furthermore, following the results of [77] we have that Inserting equations (9) and (10) to (4), the peak of the GW spectrum (right after evaporation), at k ∼ k UV , is approximately M PBH,f 10 7 g 34/9 .
where ∆N eff is the effective number of neutrino species and where we took into account the temperature dependence of the effective degrees of freedom g * (T ) in the total energy density of radiation. Hence, we can write M PBH,f 10 7 g 34/9 .
We note that g * is not significantly affected by PBH spin (e.g. [72]). The effect of PBH spin 1 will be predominantly through the evaporation lifetime t eva , where for Kerr BH with Kerr spin parameter a * = 0.9999 one has t eva (a * )/t eva (0) ∼ 0.5 [64,66]. To estimate the GW spectral density measured as a stochastic GW background by GW detectors, we just relate the GW density fraction to the total radiation density fraction by where Ω r,0 h 2 ≈ 4.18 × 10 −5 is the present radiation energy density fraction [78]. From Fig. 1 we see that PBH spin has a more dramatic effect on ∆N eff due to DR from PBH decay products compared to induced GWs. On the other hand, the induced GWs can lead to observable signatures even when PBHs are non-rotating. Further, ∆N eff from induced GWs is significantly affected by changes in PBH mass, in contrast to PBH decay emission.
While in the above we have considered a monochromatic PBH spectrum, formation models can often result in PBHs with distribution in spins and masses. While distribution in spins will significantly affect PBH decay signatures, this is not the case for induced GWs. This can be understood as follows, even for the case of Kerr BHs the lifetime will only change by a factor of ∼ 0.5. On the other hand, distribution in masses can significantly modify the PBH lifetime that goes as t eva ∝ M 3 PBH and hence the resulting induced GW signatures. The reason is that distribution in PBH lifetimes will result in a non-instantaneous effective evaporation as we considered above, which translates into more gradual deposits of evaporation radiation and less dramatic oscillations of the perturbations and suppression of GW production [77,91]. A direct detection of stochastic GW counterpart signal (see Fig. 1) would be characteristic of a narrowwidth PBH mass-function.
As shown on the right of Fig. 1, the ∆N eff signature can be testable and correlated further with the stochastic GW signals observed in GW experiments. This allows for further discrimination of the sources of the signals.
In conclusion, we have shown how coincidence signals associated with Hawking evaporation particle production as well as induced GWs can be used to explore and discriminate between different formation mechanisms for spinning and non-spinning evaporationg PBHs. This establishes a concrete methodology for charting the unexplored PBHs with masses 10 9 g. Our results are general and can be readily applied to a wide class of models.

ACKNOWLEDGMENTS
We thank Kazunori Kohri for comments. G.D. as a Fellini fellow is supported by the European Union's Horizon 2020 research and innovation programme under the Marie Sk lodowska-Curie grant agreement No 754496. This work is supported in part by the JSPS KAKENHI Nos. 19H01895 and 20H04727. M.S. and V.T. are also supported by the World Premier International Research Center Initiative (WPI), MEXT, Japan.

Appendix A: Gravitational Wave Production
Here we present the details of the formulation if PBHs form by the collapse of large primordial fluctuations in a radiation dominated universe. Let us assume that primordial fluctuations have a sharp peak at some arbitrary scale k i . This scale enters the horizon at a time when a i H i = k i and PBHs form with a monochromatic mass-spectrum at In the above we assume that the mass enclosed inside the horizon collapses into a black hole with an efficiency parameter γ. In radiation domination, one has γ ∼ 0.2 (see e.g. [27]). Starting with this initial time, PBHs start to evaporate by Hawking radiation. Evaporation: For a BH without spin, the evaporation rate goes as Hence, the evaporation time is If the PBH has spin, the evaporation is more efficient and the rate increases. We parametrise this with a spin dependence of the evaporation time, i.e. t eva (a * ), where a * is the BH spin. Assuming that PBH evaporate long after they dominate, the Hubble parameter at evaporation is given by where g * (T reh ) are the effective degrees of freedom of the energy density of the resulting radiation. The reheating (evaporation) temperature then reads In order to have successful BBN, we roughly need thermal equilibrium at T reh > 4 MeV that for a BH without (with maximal) spin imposes Note that this result depends on effective degrees of freedom at reheating g * (T rh ). Density fluctuations: PBH are produced very rarely. Only the largest primordial fluctuations from the tail of the distribution are above the necessary threshold. Thus, a good assumption is that PBHs are uniformly distributed across the universe with an initial mean PBH separation (assuming monochromatic spectrum) given by where in the last step we introduced β as the initial density fraction of PBHs, i.e. β ≡ ρ PBH,i /ρ rad,i . This uniform distribution leads to a Poisson spectrum for the density fluctuations on coarse grained scales given by [74] δρ with a (comoving) cut-off given by One can the find that the initial (isocurvature) dimensionless power spectrum is given by (A11) From the initial stage to PBH domination: At some moment, the initial fraction of PBH grows to dominate the universe. The time of radiation-PBH equality is approximately a eq /a i = β −1 . (A12) With such relation, we can compute all the relevant scales in terms of β and the evaporation time by Using the transfer function presented in the main text, one finds that at some time inside the PBH domination, the initial isocurvature fluctuations are converted to curvature perturbations with an amplitude of the power spectrum at k = k UV given by The resulting induced GW spectrum has a peak at k ∼ k UV with amplitude given by Ω GW,reh (k UV ) (A15) . This is energy density of GWs at some time right after reheating. However, we must translate this results to the energy density of GWs today or at the time of BBN. To do that we must take into account the expansion of the universe and the change of the effective number of degrees of freedom from reheating to BBN. This is calculated by using that We then find that the density fraction of GWs at BBN is given by Ω GW,BBN ≈ 0.39 g * (T reh ) 106.75 . (A18) To relate the amplitude at BBN to the amplitude measured by GW detectors, we can express it in terms of today's density parameter of radiation, Ω r,0 h 2 . In this case, the spectral density reads Ω GW,0 h 2 ≈ 0.39 Ω r,0 h 2 g * (T reh ) 106.75