Gravitational wave hints black hole remnants as dark matter

In 1987, about a decade and a half years after Hawking and Carr [1, 2] showed that black holes could form in the early Universe, MacGibbon–in a paper published in Nature [3]–already entertained the possibility that cold dark matter (DM) could be made of residues left after black hole evaporation [4, 5]. The idea that black holes might not evaporate completely already appears by Man'ko and Markov in the proceedings of a singular workshop in Moscow in 1981, called 'Seminar on Quantum Gravity' [6] (the proceedings is worth checking). From today's perspective, the possibility of black hole remnants, sometimes also called relics, might not be the most attractive end state of a black hole, albeit it is a plausible one. After all, even if Hawking evaporation leaves something behind, it is not clear how one would directly detect Planck sized objects or prove their existence. Nevertheless, their collective gravitational pull could, e.g. explain a completely invisible DM [3] and their bag-of-gold interior might offer a solution to the black hole information loss paradox [7]; the latter is still under debate [8]. We refer the reader to [9] for a thorough review on the information loss paradox and black hole remnants.

From a theoretical point of view, remnants might be a consequence of quantum gravity. There is some degree of belief that the theory of quantum gravity would be free from singularities, such as those that appear in the interior of a black hole. Within such candidate quantum gravity theories, a regular black hole might well end up in a stable state after evaporation (see e.g. [10]). For instance, regular black hole solutions have been found within loop quantum gravity [11–13], non-commutative geometry [14], limiting curvature models in mimetic gravity [15–17] and generalized uncertainty principles [18]. For more models, see the collection in [9, 19, 20]. Recently, there is also growing interest within asymptotic safe gravity (see [21, 22] for recent reviews). Note that these are static solutions and it is not clear what would occur when considering the initial collapse of matter.

From a cosmologist perspective, black hole remnants appear to be a 'bonus' to the rich phenomenology of primordial black holes (PBHs) for short (curious fact: the usage of the acronym PBH dates back to 1975 [23], but in lower case letters). The PBH scenario is nowadays a very popular topic, as can be seen from the many recent (and thorough) reviews [24–29]. We will be most concerned with tiny ($= M_{\textrm{PBH}}\lt5\times 10^8\,\textrm{g}$) PBHs, since they evaporate much before Big Bang Nucleosynthesis (BBN) [30–32] and may abundantly leave remnants. Though we will not dwell into the details of the formation of such tiny PBHs, they seem to be easily generated towards the end of inflation by preheating instabilities or quantum stochastic effects [33–35].

Tiny PBHs have an interesting early Universe cosmology: they could totally (or partially) reheat the Universe [36–39] and explain the baryon asymmetry of the Universe [40–52]. Tiny PBHs also produce high frequency gravitational waves (GWs) by Hawking evaporation [45, 53–61], PBH binaries [54, 57] and secondary GWs [57, 62–68]. See [57, 65] for a recollection of GWs associated with tiny PBHs. For the original works see Carr [36], Chapline [37], Zeldovich and Starobinsky [42, 43] and Turner [40, 41]. On top of all that, PBH remnants ⁵ could account for a fraction or all of the DM, if they exist. Any overproduction of remnants would then constrain inflationary models, e.g. see [73, 74].

Very interestingly, it has been recently noticed that some GW products of tiny PBHs might be within the range of future experiments, such as CMB-S4 [56, 58–61] and GW detectors such as the Einstein Telescope [57, 62–68] and, perhaps, high frequency GW detectors [49, 61]. In the future, we may find signatures of PBH evaporation in the early Universe. In this note, we add a unique signature to the PBH remnant scenario. We show that density fluctuations due to the initial inhomogeneous distribution of PBHs leads to (induced) GWs within a fixed low frequency range, which enters the LIGO/Virgo band and could be detected in the future by ET. The rest of the note is organized as follows. We briefly review the PBH remnant scenario as DM and the production of low frequency GWs in section 2. We then end with a short discussion in section 3. Most of the details of the formulas in this paper can be found, e.g. in [25, 47, 57, 65, 73]. When needed we use the cosmological parameters of Planck 2018 [75].

In the tiny PBH scenario we have two basic parameters: the initial mass of the PBHs at formation $M_{\textrm{PBH,f}}$ and the initial energy density fraction $\beta = \rho_{\textrm{PBH,f}}/\rho_{\textrm{total}}$ [25]. For simplicity, we will assume that PBHs form by the collapse of primordial fluctuations with a monochromatic PBH mass function. We discuss later the effects of extended mass functions. Under this assumption, $M_{\textrm{PBH,f}}$ and β are related to the Hubble horizon H and the number density of PBHs n_PBH at formation ⁶ . For a fixed $M_{\textrm{PBH,f}}$ and β we have

$$\begin{align} H_{\textrm{f}} = 4\pi\gamma \frac{M_{\textrm{pl}}^2}{M_{\textrm{PBH,f}}}\approx 5\times 10^{13} \,\textrm{GeV}\times \left(\frac{M_{\textrm{PBH,f}}}{1\,\textrm{g}}\right)^{-1}\,, \end{align} \tag{ 1 }$$

where we used that $\gamma\sim0.2$ [25] and that the reduced Planck mass is given by $M_{\textrm{pl}} = (8\pi G)^{-1/2}\approx 4.235\times 10^{18}\,\textrm{GeV}\approx 4.3\times 10^{-6}\,\textrm{g}$. For high energy scale inflation one approximately has $H_{\textrm{f}}\sim 10^{-5}M_{\textrm{pl}}$ and so at least $M_{\textrm{PBH,f}}\gt1\,\textrm{g}$. We also have that ⁷

$$\begin{align} n_{\textrm{PBH,f}} = \frac{\rho_{\textrm{PBH,f}}}{M_{\textrm{PBH,f}}} = \frac{3\beta}{4\pi\gamma}H_{\textrm{f}}^3\approx 10^{-3} \textrm{qm}^{-3}\times \beta\left(\frac{1\,\textrm{g}}{M_{\textrm{PBH,f}}}\right)^{3}\,, \end{align} \tag{ 2 }$$

where we used $\rho_{\textrm{total,f}} = 3H_{\textrm{f}}^2M_{\textrm{pl}}^2$ and equation (1).

After formation, tiny PBHs quickly evaporate. The evaporation time reads [57]

$$\begin{align} t_{\textrm{eva}}\approx \frac{160\alpha}{3.8\pi g_H(T_{\textrm{PBH}})}\frac{M_{\textrm{PBH,f}}^3}{M_{\textrm{pl}}^4}\approx 400\,\textrm{qs}\times \alpha \left(\frac{M_{\textrm{PBH,f}}}{1\,\textrm{g}}\right)^{3}\,, \end{align} \tag{ 3 }$$

where $g_H(T_{\textrm{PBH}})$ is the spin-weighted degrees of freedom and we introduced the parameter α to take into account the effect of the PBH spin. For no spin α = 1, while for a near extremal PBH $\alpha\sim1/2$ [76]. In deriving equation (3) we assumed that the evaporation time is much larger than the formation time and that $g_H(T_{\textrm{PBH}})\approx 108$ [56]. If we compare the evaporation time in equation (3) with the Hubble time at formation, equation (1), in a radiation-dominated Universe,

$$\begin{align} t_{\textrm{f}}\approx\frac{1}{2H_{\textrm{f}}}\approx 10^{-8} \,\textrm{qs}\times \left(\frac{M_{\textrm{PBH,f}}}{1\,\textrm{g}}\right)\,, \end{align} \tag{ 4 }$$

we indeed see that $t_{\textrm{eva}}\gg t_{\textrm{f}}$.

Many of the discussions that follow will depend on the ratio $t_{\textrm{eva}}/t_f$, which tells us how long these tiny PBHs stayed around. So, let us write it explicitly:

$$\begin{align} R_{\textrm{eva, f}}\equiv\frac{t_{\textrm{eva}}}{t_{\textrm{f}}} = \frac{1280\pi\gamma\alpha }{3.8\pi g_H(T_{\textrm{PBH}})}\frac{M_{\textrm{PBH,f}}^2}{M_{\textrm{pl}}^2}\approx 3\times 10^{10}\,\alpha\left(\frac{M_{\textrm{PBH,f}}}{1\,\textrm{g}}\right)^{2}\,. \end{align} \tag{ 5 }$$

For instance, the condition for the PBH dominance is given by

$$\begin{align} \beta > \beta_{\textrm{min}} = 1/\sqrt{R_{\textrm{eva, f}}}\approx 6\times 10^{-6}\alpha^{-1/2}\left(\frac{M_{\textrm{PBH,f}}}{1\,\textrm{g}}\right)^{-1}\,. \end{align} \tag{ 6 }$$

If $\beta\lt\beta_{\textrm{min}}$, PBHs never dominate. The minimum abundance β_min is obtained by requiring that $H_{\textrm{eeq}}\gt H_{\textrm{eva}}$, where H_eeq and H_eva respectively refer to the Hubble parameter at the time of early radiation-PBH equality after PBH formation and the Hubble parameter right after the complete PBH evaporation. We then used that $H_{\textrm{eeq}}/H_{\textrm{f}}\approx\sqrt{2}\beta^2$ and $H_{\textrm{eva}}\sim1/t_{\textrm{eva}}$. The $\sqrt{2}$ in H_eeq appears after using the exact solutions for a radiation-matter Universe (see, e.g. equation (1.81) in Mukhanov's book [77]) which is a good approximation since PBHs evaporate almost instantaneously [57]. We can also find the amount of the expansion of the Universe from PBH formation until PBH evaporation,

$$\begin{align} \frac{a_{\textrm{eva}}}{a_{\textrm{f}}}\approx \left\{ \begin{aligned} &{\sqrt{R_{\textrm{eva,f}}}}\approx {2\times 10^{5}}\,{\alpha^{1/2}}\left(\frac{M_{\textrm{PBH,f}}}{1\,\textrm{g}}\right);\qquad&\beta \lt \beta_{\textrm{min}}\,,\\ &\left(\frac{3\beta^{1/2}R_{\textrm{eva,f}}}{4}\right)^{2/3}\approx {9\times10^{6}}\,{\beta^{1/3}\alpha^{2/3}}\left(\frac{M_{\textrm{PBH,f}}}{1\,\textrm{g}}\right)^{4/3};\qquad&\beta \gt \beta_{\textrm{min}}\,, \end{aligned} \right. \end{align} \tag{ 7 }$$

where, for the case $\beta\gt\beta_{\textrm{min}}$, we used the fact that $H = 2/(3t)$ in a matter dominated Universe. For a radiation dominated Universe, one instead has $H = 1/(2t)$. We see that larger PBH masses lead to larger number of e-folds. The β dependence in the case $\beta\gt\beta_{\textrm{min}}$ can be understood from the fact that the early radiation-PBH equality depends on β.

We can also compute the temperature of the radiation filling the Universe after PBH evaporation. This is given by

$$\begin{align} T_{\textrm{eva}}\approx 2.4 (2.8)\times 10^{10}\,\textrm{GeV}\,\alpha^{-1/2}\left(\frac{M_{\textrm{PBH,f}}}{1\,\textrm{g}}\right)^{-3/2}\left(\frac{g_{s\star}(T_{\textrm{eva}})}{106.75}\right)^{-1/4}\,, \end{align} \tag{ 8 }$$

where the coefficient $2.4(2.8)$ is for the case $\beta\lt\beta_{\textrm{min}}$($\beta\gt\beta_{\textrm{min}}$), $g_{s\star}$ is the effective degrees of freedom for the entropy and we used the Planck 2018 [75] value for H_eq. Later $g_{\star}$ denotes the effective degrees of freedom for the energy density. To evaluate them we use the fitting formulas provided in [78]. Evaporation before BBN, that is $T_{\textrm{eva}}\gt4\,\textrm{MeV}$ [30–32], requires $M\lt5\times 10^8\,\textrm{g}$. That is all we need to understand the PBH remnant scenario and the associated low frequency GWs.

2.1. GWs after PBH domination and evaporation

PBH formation is rather a rare event. Only those Hubble patches with high enough density contrast will collapse; for Gaussian fluctuations this is exponentially unlikely. Thus, the spatial distribution of PBH formation is to a good approximation uniformly random. This leads to a Poisson type spectrum of PBH number density isocurvature fluctuations [62] (see figure 1 for an illustration). Note, however, that the coarse grained, fluid picture for the collection of PBHs breaks down at the mean inter-PBH comoving separation, which is given by

$$\begin{align} k_{\textrm{uv}}\equiv a_{\textrm{f}}\bar d_{\textrm{f}}^{-1} = \left(\frac{4\pi}{3}n_{\textrm{PBH, f}}a_{\textrm{f}}^3\right)^{1/3} = \frac{\beta^{1/3}}{\gamma^{1/3}}a_{\textrm{f}}H_{\textrm{f}}\,. \end{align} \tag{ 9 }$$

Equation (9) gives us the high momenta cut-off for the spectrum of number density fluctuations. This is also where fluctuations are larger and the spectrum of (PBH-isocurvature) induced GWs peaks ⁸ . We consider from now on, in this subsection, only the case where PBHs dominate the Universe, i.e. $\beta\gt\beta_{\textrm{min}}$. Otherwise, the production of induced GWs is very much suppressed and one needs very large isocurvature fluctuations [79].

Zoom In Zoom Out Reset image size

Curiously, when PBHs dominate, there is a β independent relation between k_uv and k_eva, which is given by [63]

$$\begin{align} \frac{k_{\textrm{uv}}}{k_{\textrm{eva}}} = \frac{\beta^{1/3}}{\gamma^{1/3}}\frac{H_{\textrm{f}}}{H_{\textrm{eva}}}\frac{a_{\textrm{f}}}{a_{\textrm{eva}}} = \left( \frac{3R_{\textrm{eva,f}}}{4\gamma}\right)^{1/3}\approx 5000\,\alpha^{1/3}\left(\frac{M_{\textrm{PBH,f}}}{1\,\textrm{g}}\right)^{2/3}\,. \end{align} \tag{ 10 }$$

This allows us to easily find the peak frequency of the GWs today by using that

$$\begin{align} f_{\textrm{eva}} = \frac{k_{\textrm{eva}}}{2\pi a_0}&\approx 3\times 10^{-5}\,\textrm{Hz}\left(\frac{T_{\textrm{eva}}}{1 \textrm{TeV}}\right)\left(\frac{g_{\star}(T_{\textrm{eva}})}{106.75}\right)^{1/2}\left(\frac{g_{s\star}(T_{\textrm{eva}})}{106.75}\right)^{-1/3}\nonumber\\[5pt] &\approx 734 \, \textrm{Hz}\,\alpha^{-1/2}\left(\frac{M_{\textrm{PBH,f}}}{1\,\textrm{g}}\right)^{-3/2}\left(\frac{g_{\star}(T_{\textrm{eva}})}{106.75}\right)^{1/4}\left(\frac{g_{s\star}(T_{\textrm{eva}})}{106.75}\right)^{-1/3}\gamma\,, \end{align} \tag{ 11 }$$

which yields

$$\begin{align} f_{\textrm{uv}} \approx 3.6\times 10^{6}\,\textrm{Hz}\,\alpha^{-1/6}\left(\frac{M_{\textrm{PBH,f}}}{1\,\textrm{g}}\right)^{-5/6}\left(\frac{g_{\star}(T_{\textrm{eva}})}{106.75}\right)^{1/4}\left(\frac{g_{s\star}(T_{\textrm{eva}})}{106.75}\right)^{-1/3}\,. \end{align} \tag{ 12 }$$

The calculation on the dominant contribution to the amplitude of GWs can be found in [63] (see also [62, 80, 81]). Here we only sketch the production of induced GWs and their amplitude after a PBH dominated Universe. Let us first give the result and then explain it. The induced GW spectrum from a PBH dominated Universe evaluated today is given by

$$\begin{align} \Omega_{\textrm{GW},0}h^2 = 1.62\times 10^{-5}\,c_g\,\Omega^{\textrm{peak}}_{\textrm{GW,eva}}\left(\frac{k}{k_{\textrm{uv}}}\right)^{11/3}\Theta(k-k_{\textrm{uv}})\,, \end{align} \tag{ 13 }$$

where

$$\begin{align} c_g\equiv\left(\frac{\Omega_{\textrm{r},0}h^2}{4.18\times 10^{-5}}\right)\left(\frac{g_{\star}(T_{\textrm{eva}})}{106.75}\right)\left(\frac{g_{s\star}(T_{\textrm{eva}})}{106.75}\right)^{-4/3}\,, \end{align} \tag{ 14 }$$

which takes into account the redshift of the GW energy density from evaporation until today as well as the change of relativistic degrees of freedom, $\Omega_{\textrm{r},0}h^2$ is the energy fraction of radiation today from Planck 2018 [75], and

$$\begin{align} \Omega^{\textrm{peak}}_{\textrm{GW,eva}}\approx \frac{\beta ^{16/3} \gamma ^{8/3} }{1536\,\times {2}^{1/3} \sqrt{3} \pi }\left(\frac{k_{\textrm{uv}}}{k_{\textrm{eva}}}\right)^{17/3}\approx\frac{\alpha^{17/9}}{10} \left(\frac{\beta}{10^{-3}}\right)^{16/3}\left(\frac{M_{\textrm{PBH,f}}}{1\,\textrm{g}}\right)^{34/9}\,, \end{align} \tag{ 15 }$$

which gives the peak amplitude right after PBH evaporation.

Now, let us roughly explain the origin of the amplitude equation (15). First of all, the factor $\beta^{16/3}$ is the suppression due to the fact that PBHs are formed during the radiation era and the amplitude of curvature fluctuations decay until PBHs dominate the Universe. The suppression goes as ${\sim}(k_{\textrm{eeq}}/k_{\textrm{uv}})^2\propto \beta^{4/3}$, where $k_{\textrm{eeq}} = a_{\textrm{eeq}}H_{\textrm{eeq}}$ is the mode that enters the horizon at the early radiation-PBH equality. Induced GWs, as a secondary effect, are proportional to the four-point function of curvature fluctuations. Hence, we obtain $(\beta^{4/3})^4 = \beta^{16/3}$. The factor involving $k_{\textrm{uv}}/k_{\textrm{eva}}$ is more interesting. During PBH domination, the curvature perturbation is constant on all scales and density fluctuations grow proportional to the scale factor. On scales corresponding to k_uv, the PBH number density fluctuations at evaporation have grown by a factor $(k_{\textrm{uv}}/k_{\textrm{eva}})^2$, which is very large. And then, PBHs almost suddenly evaporate: huge pressureless density fluctuations are converted into huge radiation fluctuations which generates a huge wake in the velocities of the radiation fluid. For the induced GW spectrum, which is proportional to four gradients of velocities, we get $((k_{\textrm{uv}}/k_{\textrm{eva}})^2)^4 = (k_{\textrm{uv}}/k_{\textrm{eva}})^8$. To go from 8 to $17/3$ one needs to account that only a very narrow window of scalar modes contribute to the integral, that is only those modes close enough to the resonance in the GW kernel, which is proportional to $k_{\textrm{eva}}/k_{\textrm{uv}}$ [57, 64, 80], this gets us to 7. The rest comes from the 'almost' sudden evaporation [57]: very short wavelength modes actually feel the time dependence of Hawking evaporation, which goes as $M_{\textrm{PBH}}(t)\sim (1-t/t_{\textrm{eva}})^{1/3}$, and suppresses curvature fluctuations by a factor $(k_{\textrm{eva}}/k_{\textrm{uv}})^{1/3}$. Thus, we have $8-1-4\times1/3 = 17/3$. We can then use equation (15) to place upper bounds on β from BBN constraints, namely [63]

$$\begin{align} \beta < 10^{-3}\left(\frac{M_{\textrm{PBH,f}}}{1\,\textrm{g}}\right)^{-17/24}\,. \end{align} \tag{ 16 }$$

It is important to emphasize that equation (15) should be understood as a rough order of magnitude estimate, because at some point during PBH domination, the number density fluctuations exceed unity. Furthermore, the amplitude is very sensitive to the width of the PBH mass function. For a log-normal with logarithmic width $\sigma\sim 1$ the amplification is negligible [57]. In any case, let us note that for the parameters of interest, the curvature perturbation and its time derivative are always smaller than unity and well within the perturbative regime.

2.2. PBH remnants as DM

Let us collect all the previous results and assume that PBH evaporation leaves behind Planck relics with mass

$$\begin{align} m_{\textrm{relic}} = r M_{\textrm{pl}}\,, \end{align} \tag{ 17 }$$

where r > 1 is a free parameter. The fact that evaporation stops when $M_{\textrm{PBH}}(t_{\textrm{end}}) = m_{\textrm{relic}}$ does not affect the calculations in section 2.1. That is, $t_{\textrm{end}} = t_{\textrm{eva}}-\delta t$, where $\delta t/t_{\textrm{eva}} = O(t_{\textrm{eva}}M_{\textrm{pl}})$, which is negligibly tiny. Thus, let us use the previous results and require that PBH remnants occupy a fraction f_relic of the total DM today, that is

$$\begin{align} f_{\textrm{relic}}\equiv\frac{\rho_{\textrm{relic}}}{\rho_{\textrm{DM}}}\,. \end{align} \tag{ 18 }$$

Extrapolating backwards from today until evaporation, using that $\rho_{\textrm{DM}}\propto a^{-3}$, we have that

$$\begin{align} \Omega_{\textrm{relic}}\big|_{\textrm{eva}}&\approx7\times 10^{-13}f_{\textrm{relic}}\left(\frac{T_{\textrm{eva}}}{1\textrm{TeV}}\right)^{-1}\left(\frac{g_{s\star}(T_{\textrm{eva}})}{106.75}\right)\left(\frac{g_{\star}(T_{\textrm{eva}})}{106.75}\right)^{-1}\nonumber\\ &\approx 2.9 \,(2.5)\times 10^{-20}f_{\textrm{relic}}\sqrt{\alpha}\left(\frac{M_{\textrm{PBH,f}}}{1\textrm{g}}\right)^{-3/2}\left(\frac{g_{s\star}(T_{\textrm{eva}})}{106.75}\right)\left(\frac{g_{\star}(T_{\textrm{eva}})}{106.75}\right)^{-3/4}\,, \end{align} \tag{ 19 }$$

where again the value between brackets is for $(\beta\gt\beta_{\textrm{min}})$ and we used that $\rho_{\textrm{DM,eq}} = 3 H_{\textrm{eq}}^2M_{\textrm{pl}}^2/2$. If we now extrapolate forwards from the PBH formation, we have that at evaporation

$$\begin{align} \Omega_{\textrm{relic}}\big|_{\textrm{eva}} = \frac{\beta m_{\textrm{relic}}}{M_{\textrm{PBH,f}}}\frac{H_{\textrm{f}}^2}{H_{\textrm{eva}}^2}\frac{a_{\textrm{f}}^3}{a_{\textrm{eva}}^3} \approx\left\{ \begin{aligned} & \dfrac{m_{\textrm{relic}}}{M_{\textrm{PBH,f}}}{\beta}\sqrt{R_{\textrm{eva,f}}}\approx r\beta\sqrt{\alpha}&(\beta < \beta_{\textrm{min}})\\ & \dfrac{m_{\textrm{relic}}}{M_{\textrm{PBH,f}}}&(\beta > \beta_{\textrm{min}}) \end{aligned} \right.\,.\quad \end{align} \tag{ 20 }$$

In this way, we can draw the parameter space where PBH remnants can be a fraction f_relic of DM. In general we find that for a fixed f_relic,

$$\begin{align} {M_{\textrm{PBH,f}}}\approx10^6\mathrm{\,g}\left\{ \begin{aligned} & 2.4\dfrac{r^{2/3}}{f^{\,2/3}_{\textrm{relic}}}\left(\frac{\beta}{10^{-10}}\right)^{2/3}\left(\frac{g_{s\star}(T_{\textrm{eva}})}{106.75}\right)^{-2/3}\left(\frac{g_{\star}(T_{\textrm{eva}})}{106.75}\right)^{1/2}&(\beta < \beta_{\textrm{min}})\\ & 0.5\dfrac{r^{2/5}}{f^{\,2/5}_{\textrm{relic}}\alpha^{1/5}} \left(\frac{g_{s\star}(T_{\textrm{eva}})}{106.75}\right)^{-2/5}\left(\frac{g_{\star}(T_{\textrm{eva}})}{106.75}\right)^{3/10} &(\beta > \beta_{\textrm{min}}) \end{aligned} \right.\,. \end{align} \tag{ 21 }$$

We show the parameter space in figure 2.

Zoom In Zoom Out Reset image size

Note that, if there was a PBH dominated stage in the early Universe ($\beta\gt\beta_{\textrm{min}}$), and if the PBH remnants totally account for the DM ($f_{\textrm{relic}} = 1$), the initial PBH mass is uniquely determined as

$$\begin{align} M_{\textrm{PBH,f}}\approx5 \times 10^{5}\mathrm{\,g}\times \frac{r^{2/5}}{\alpha^{1/5}}\,. \end{align} \tag{ 22 }$$

This correspond to an evaporation temperature of $T_{\textrm{eva}}\approx80\,\textrm{GeV}$. Note that this value is in agreement with [47]. This case corresponds to an induced GW signal with peak at

$$\begin{align} f_{\textrm{uv}}\approx 70\,\textrm{Hz}\times \alpha^{-1/6}\,. \end{align} \tag{ 23 }$$

This frequency falls well within the frequency range of future GW detectors such as LIGO A+, Voyager, Einstein telescope, cosmic explorer and DECIGO. The frequency (23) corresponds to an inter-PBH comoving separation of 600 km, which is also the mean separation between the remnants.

The PBH reheating scenario with PBH remnants has a unique prediction for the peak frequency of the induced GW spectrum. For a fixed PBH mass, the amplitude of the GW spectrum (15) only depends on β and its value today is given by

$$\begin{align} \Omega^{\textrm{peak}}_{\textrm{GW},0}h^2\approx 10^{-11}r^{68/45}\alpha^{17/15}\left(\frac{\beta}{10^{-8}}\right)^{16/3} \,. \end{align} \tag{ 24 }$$

We plot the induced GW spectrum in figure 3. Recall that $1\gt\alpha\gt1/2$ and $r\sim O(1)$, so the predictions (22)–(24) do not depend much on whether PBHs have spin or if the remnant is a bit larger than the Planck length. Although this signal could also be present without remnants, detecting a peak right at this frequency would be a strong indication that the PBH remnants is the DM.

Zoom In Zoom Out Reset image size

The existence of remnants after Hawking evaporation is suggested in some theories of quantum gravity [11, 21, 22]. The remnants could play an important role in the information loss paradox and in cosmology [9]. Here we focused on the possibility that the Universe is filled with the remnants of tiny PBHs which evaporated well before BBN. For some parameter space the PBH remnants could account for all the DM and reheat the Universe [3, 18, 47].

One of the problems of the PBH reheating plus DM remnants scenario was that it seemed almost impossible to probe. In this note, we pointed out that it may lead to a unique prediction for the GW background: a peaked signal at a frequency ∼70 Hz, close to the peak sensitivity of LIGO/Virgo/KAGRA, LIGO A+, voyager Einstein telescope and cosmic explorer. The low frequency tail of the resulting GW spectrum would also be seen in DECIGO but it is out of reach for LISA. While this is not definitive evidence of remnants as DM, finding a peak at such a precise frequency would give a strong indication of the PBH reheating plus remnants scenario. This could be further probed by additional signatures of high frequency GWs and the effective number of species [49, 58, 60, 61, 91, 92].

A remaining issue is to derive a more accurate estimate for the amplitude of PBH isocurvature induced GWs, because the PBH number density fluctuations reach the nonlinear regime close to the final stage of evaporation [57, 63]. However, this requires sophisticated numerical simulations. Another issue is the effect of a finite width in the PBH mass function. While these issues might reduce the amplitude of the induced GW spectrum, the peak frequency would not be significantly affected. Thus the prediction for a GW background peaked at ∼70 Hz seems robust in the scenario of PBH remnants as DM within O(1) factors.

We would like to thank the anonymous referee for the many useful comments. G D thanks A D Rojas, C D Rojas and D Rojas for their support and A Platania for useful correspondence on black holes in asymptotic safe gravity. We also thank P Bari and I Markin for useful comments on the paper. G D is supported by the DFG under the Emmy-Noether program Grant No. DO 2574/1-1, Project Number 496592360. Kavli IPMU is supported by World Premier International Research Center Initiative (WPI), MEXT, Japan. This work was supported in part by the JSPS KAKENHI Grant Nos. 19H01895, 20H04727, and 20H05853.

All data that support the findings of this study are included within the article (and any supplementary files).