Probing Primordial Black Hole Scenarios with Terrestrial Gravitational Wave Detectors

It is possible that primordial black holes consitute (or consituted) a significant fraction of the energy budget of our universe. Terrestrial gravitational wave detectors offer the opportunity to test the existence of primordial black holes in two different mass ranges, from $10^2\,{\rm g}-10^{16}\,{\rm g}$ to $10^{-6}\,M_\odot-100 \,M_\odot$. The first mass window is open via induced gravitational waves and the second one by gravitational waves from binary mergers. In this review, we outline and explain the different gravitational wave signatures of primordial black holes that may be probed by terrestrial gravitational wave detectors, such as the current LIGO/Virgo/KAGRA and future ones like Einstein Telescope and Cosmic Explorer. We provide rough estimates for the frequency and amplitude of the associated GW background signals. We also discuss complementary probes for these primordial black hole mass ranges.


I. INTRODUCTION
Our universe could contain more black holes than those expected from the result of stellar evolution.Such non-astrophysical black holes, so-called Primordial Black Holes (PBHs1 ), could have formed in the very early universe by an event of cosmological proportions.After formation, a PBH decays via Hawking radiation [2,3], with a lifetime proportional to the cube of the initial mass.A crude estimate then tells us that PBHs heavier than 10 15 g are long-lived and still with us today (and so have a lifetime larger than the age of the universe), while lighter PBHs have already evaporated. 2See, e.g., Refs.[4][5][6][7][8][9] for recent reviews on PBHs.
As surprising as it may seem, PBHs could be (or could have been) an essential component of our universe, according to current observations [6,10].On one hand, long-lived PBHs could explain the dark matter [11][12][13][14][15][16][17], the microlensing events reported by HSC [18] and OGLE [19][20][21], some of TABLE I. Summary of interesting mass ranges for various PBH scenarios, such as PBH reheating, PBH dark matter, HSC and OGLE microlensing events, LVK merger events and seeds of SuperMassive Black Holes (SMBHs).The mass range in the PBH reheating scenario follows from viability: 1 g PBHs form right after GUT scale inflation and 10 9 g PBHs reheat the universe too close to nucleosynthesis.Note that under the null hypothesis of no PBHs, HSC [18] and OGLE [20] respectively probe the mass ranges 10 −11 − 10 −6 M ⊙ and 10 −6 − 10 −3 M ⊙ .Most interestingly, HSC found one candidate event with a mass of about 10 −8 M ⊙ and OGLE found 6 candidates with masses around 10 −4 M ⊙ .A solar mass is M ⊙ ≈ 2 × 10 33 g.
The focus of this review are GW signals of PBHs accessible to terrestrial GW detectors, which happen to include the PBH reheating scenario and PBHs with masses ranging from an earth mass to hundred solar masses.Interestingly, there are tentative hints of solar and sub-solar mass compact objects in the LVK data [106][107][108][109] and solar-mass objects in quasar broad emission lines [110]. 4 There are also the moon-mass and earth-mass microlensing events respectively reported by HSC [18] and OGLE [19,20].For simplicity, although our discussion may apply to more general situations, 5we will specialize to the case of PBH formation from large primordial fluctuations.In that case, the planet to solar mass window has direct connections to the recent Pulsar Timing Arrays (PTAs) results on a possible GW background [112][113][114][115][116][117][118][119][120][121], if interpreted as the induced GW signal  (see also Refs.[160][161][162] for the merger of supermassive PBHs).
This review is organized as follows.In § II, we explain the various GW signals associated with PBHs and provide order of magnitude estimates for the frequency and amplitude.In § III, we focus on the potential of terrestrial GW detectors to probe PBH scenarios.Then, in § IV we briefly discuss other complementary probes to terrestrial GW detectors.We conclude our work with further discussions in § V. Useful references to dig into the details of PBHs from the collapse large primordial fluctuations are Ref.[5] for a general review, Ref. [6] for an extensive review on available constraints, Ref. [8] for the details of PBH formation and Ref. [9] for the connection to inflation.

II. GW SIGNALS OF BLACK HOLES FROM THE EARLY UNIVERSE
Suppose that in the early, radiation dominated universe there is a large positive fluctuation in the spatial curvature of the universe on scales larger than the Hubble radius (sometimes also called Hubble horizon).If such a positive curvature fluctuation is as large as the square of the Hubble parameter when its comoving scale becomes comparable to the Hubble radius, the Hubble volume of that region behaves like a closed universe.And, a closed universe eventually collapses on itself.But, seen from outside of that particular Hubble volume, it is the formation of a PBH.For more details on this nice geometrical picture see Ref. [5].Consequently the mass of the resulting PBH will be proportional to that contained inside the Hubble radius at formation, say H −1 f where H f is the expansion rate of the universe at formation.Explicitly, we have that where γ is the efficiency factor usually taken to be γ ∼ 0.2 [5] and in the last step we wrote the Hubble radius in terms of the temperature T f of the radiation filling the early universe at formation.g ρ (T f ) is the effective number of massless degrees of freedom in the energy density of the primordial plasma.Whenever needed we assume the standard model of particles and use the values of Ref. [163].As a curiosity, note that the Schwarschild radius of the PBH is of the order of the Hubble radius at formation, namely r PBH,f = 2GM PBH,f = γ/H f .The Hubble radius H f , or the temperature T f , is also related to the comoving scale k f of the relevant fluctuation, which in turn is related to the time of generation during inflation and frequency of the induced GWs.We provide an estimate below.Superhorizon curvature fluctuations are, in fact, the initial conditions of the standard cosmological model.They explain the Cosmic Microwave Background (CMB) anisotropies and later give rise to galaxies and other structures we see today in our universe.CMB observations also provide strong evidence that such primordial fluctuations were generated during a period of cosmic inflation [164] (see, e.g., [165,166] for recent brief reviews focused on GWs and other future prospects).The general and simplest prediction of inflation is that of random Gaussian superhorizon curvature fluctuations with an almost scale invariant spectrum, which actually arise from vacuum quantum fluctuations [167][168][169].So, strictly speaking, the universe is filled with fluctuations of all sizes and amplitudes.The issue is that the amplitude of the primordial spectrum measured by the CMB is about 10 −9 [170] and the probability that one fluctuation is large enough to form a PBH is absurdly exponentially suppressed.But, new physics during inflation can enhance the spectrum of fluctuations up to an amplitude of 10 −2 , yielding an interesting fraction of PBHs in the universe. 6nfortunately, we will not dwell into the physics during inflation, which would deserve a whole new chapter.Instead we refer the reader to Ref. [9] (and references therein) for a review of inflation focused on PBHs.
Although there has been substantial progress in recent years in the understanding of when and how PBHs form, the precise condition for the PBH formation is still under study (see e.g., [185][186][187]).Since the PBH formation is not a main topic of this review, for the sake of simplicity, here we consider the PBH formation in the naive Press-Schechter formalism [188].In this formalism, we assume that fractional density fluctuations, δ ≡ δρ/ρ, above a certain threshold, δ th , form PBHs. The abundance of PBHs is then given by the probability to have δ above such a threshold [5], namely, where γ is the efficiency factor, σ M PBH,f is the variance of fluctuations smoothed over the comoving scale R f (basically 1/k f ), and in the last step we assumed that σ M PBH,f ≪ δ th .Note that we assumed that the statistics of the primordial fluctuations are Gaussian.The smoothed variance of fluctuations reads In Eq. (2.3), W (kR f ) is the window function, P R the primordial spectrum of curvature fluctuations R and in the last step we use the relation between density fluctuations and curvature fluctuations on superhorizon scales in a radiation dominated universe [5].The smoothing of the fluctuations is necessary to go from the Fourier space results to real space, focusing only on the relevant scale.From Eqs. (2.2) and (2.3) it is clear that the initail fraction of PBHs, β, is exponentially sensitive to the primordial spectrum P R .For long-lived PBHs, we may relate the initial fraction of PBHs with the fraction of dark matter in the form of PBHs, f PBH .Namely, where g s (T f ) are the effective degrees of freedom in the entropy.Note that the initial fraction β may be initially very small but the fraction of PBHs in the universe grows in time, since the energy density of radiation dilutes faster than that of PBHs.For evaporated PBHs we shall explore their cosmology separately later.As a word of caution, we would like to stress again that the Press-Schechter formalism is the simplest estimate and we believe that, in view of the current state-of-the-art, it should not be used to accurately predict the fraction of PBHs.We refer the reader to Refs.[185][186][187] for reviews on recent advancements in PBH formation criteria.Here we solely mention that the fraction of PBHs is highly dependent on the precise value of the threshold (which also depends on the equation state of the universe [189,190] and the radial profile of the fluctuation [191] but a typical value is δ th ∼ 0.45), the tail of the probability distribution of curvature fluctuations, the critical collapse, the detailed formalism and the window function one chooses.For example, non-Gaussianity of primordial fluctuations drastically changes the amount of PBHs produced [192][193][194][195][196].This effect is particularly important in the induced GW interpretation of the PTA data [124,129,131,135], as assuming Gaussian fluctuations tends to predict too many PBHs (though we stress that there are still uncertainties in the calculations of the exact PBH fraction).Fortunately, as we show later, the GW signals are not sensitive to the non-linear physics of PBH formation.For earlier works on the impact of local non-Gaussianities on the induced GWs see Refs.[155,[197][198][199][200][201][202].See also Ref. [203] for an analysis of LVK data including non-Gaussianities in the induced GWs.

A. GW signatures of PBHs
We proceed to describe the GWs associated with PBHs and provide estimates for their frequency and amplitude.Before doing that though, we would like to note that there are other reviews on GW signatures of PBHs.For instance, there is an extensive review by the LISA cosmology working group [204], with focus on LISA capabilities to test PBH scenarios (see also Ref. [205]).And a pedagogical and intuitive introduction can be found in Ref. [206].Here we aim for a complementary and concise summary focused on terrestrial GW detectors, such as LVK, Einstein Telescope and Cosmic Explorer.For simplicity and analytical viability, we assume a monochromatic (or almost monochromatic) primordial spectrum of curvature fluctuations and PBH mass function.When pertinent, we also discuss and cite works on broad mass functions.We also assume the standard model of particle physics and take appropriate values for the effective degrees of freedom when necessary.To recover the dependence on the effective degrees of freedom we will refer the reader to the relevant references.
We classify the GW signals associated with PBHs into: (i) GWs associated with PBH formation, (ii) GWs associated with PBH reheating, and, (iii) GWs from PBH binary mergers.
The first two signals are induced GWs, while the third are typical GWs from black hole binaries (either resolved or unresolved).At this point, it is important to clarify that, while GWs from PBH binaries (and Hawking evaporation) come from PBHs themselves, that is not the case for induced GWs.Induced GWs [207][208][209][210][211] (see Refs. [212][213][214] for recent reviews), sometimes also called secondary GWs, are a consequence of the cosmological process that led to PBH formation (if associated with PBH formation) or the cosmological process that resulted from the complete PBH evaporation (if associated with PBH reheating).A more concrete explanation for the GWs associated with PBH formation is as follows.In order to form PBHs there must be highly inhomogeneous concentrations of matter involved in the early universe, resulting in density waves and large anisotropic stresses.The latter are responsible for the generation of induced GWs.A similar logic applies after PBH reheating.Before proceeding with some estimates for the frequency and amplitude of the associated GWs, let us clarify the notation in what follows.When dealing with induced GWs, we provide the amplitude of the GW spectrum evaluated at a time when a given GW frequency is sufficiently inside the horizon to be regarded as a proper wave and, therefore, as a radiation fluid.After that epoch the energy density ratio of GWs, defined by Ω GW = ρ GW /ρ total where ρ means energy density (see, e.g., Ref. [215] for the definition of the energy density of GWs in cosmology), is mostly constant in the early universe; it only changes when there is a change in the effective degrees in energy density and entropy.When needed, we take this time to be the epoch of horizon crossing (given by k = aH where k is the comoving wavenumber and a the scale factor), although more realistic estimates suggest it may be at least 2 e-foldings later [216,217].We note however that this does not matter much for the GWs we consider, since in the very early universe the number of the effective degrees of freedom is basically constant.The spectral density of GWs today can then be written as [213,218] where Ω GW, * is the spectral density of GWs well enough inside the Hubble horizon and Ω rad,0 h 2 is the energy density fraction of radiation today given by Planck [170].Note that, for convenience, we drop the subscript * in the estimates below.We provide a summary of our estimates for the frequency and amplitude evaluated today via Eq.(2.5) in Tab.II.
(i) GWs associated with PBH formation.Large enough and rare primordial fluctuations collapse to form PBHs. The other, not large enough, average 7 primordial curvature fluctuations generate density waves with amplitudes damped in time by the primordial plasma.These are the main source of induced GWs. 8 Thus, the largest production of induced GWs occurs at around the time of Hubble horizon crossing of the typical primordial fluctuation, say with comoving wavenumber k p .For an almost monochromatic primordial spectrum, k p determines the frequency at which the spectrum of the GW background peaks, as well as the typical mass of the PBHs by Eq. (2.1) (and using that k f = k p ).The frequency evaluated today is given by 9 where we used Eq.(2.1) to replace k f with M PBH,f since it is more convenient for our discussions.
7 By average we have in mind the root mean square of primordial fluctuations. 8There are also scalar-tensor induced GWs [219][220][221][222] which are often subdominant but might have distinctive features [221]. 9A more detailed formula including the effective degrees of freedom in energy density and entropy is given by The prefactors are important for PBH masses above 10 −3 M⊙ as the corresponding temperature is below 1 GeV.However, this only concerns µHz-nHz frequencies and, therefore, it is not relevant for this review.

Peak frequency (f peak [kHz])
Peak amplitude and spectral index today Reheating I * * (Isocurvature iGWs) 1.7 2 × 10 −7  ).Note that all GW signals have a sharp cut-off for f > f peak , since we assume a monochromatic PBH mass function.The power-law index for f < f peak is given by the value of α.Ω peak GW,0 h 2 is the peak frequency of the spectral density of GWs evaluated today using Eq.(2.5).
The amplitude of the GW spectrum depends on the primordial spectrum of curvature fluctuations.For concreteness, we assume a log-normal spectrum given by where A R is the amplitude and ∆ is the logarithmic width of the spectrum.One recovers the Dirac delta case in the limit when ∆ → 0. The log-normal spectrum allows for nice analytical approximations for the induced GW spectrum, as derived in Ref. [223].For our purposes though, a good enough estimate for the amplitude at the peak is given by where the factor O(10) corresponds to sharp primordial spectrum with ∆ < 0.1.We note, however, that the estimate for the peak amplitude (2.9) is also valid in more general situations.For instance, for a scale invariant primordial spectrum, the GW spectrum is also scale invariant with amplitude ∼ 0.8A 2 R .The main difference is the spectral shape.For a sharp peak, there is a sharp cut-off at f ∼ 2f formation (above which there are no more GWs produced by momentum conservation) and a low frequency tail going as Ω GW ∝ f 2 which transitions to Ω GW ∝ f 3 when far enough from the peak, roughly for f < 2∆ × f formation [223].It is interesting to note that the f 3 is a universal infra-red scaling for GWs from localized sources [224] but induced GWs in radiation domination contain a logarithmic correction [197,225].We also note that these estimates depend on the expansion history of the universe, e.g., on the equation of state of the early universe [226,227].We show the shape of the induced GW spectrum for a log-normal (2.8) with ∆ = 0.1 in blue in Fig. 2. For general semi-analytical formulas of induced GWs see Refs.[213,[226][227][228][229].
(ii) GWs associated with PBH reheating.In two nice papers, Inomata et al. [230,231] studied the induced GWs generated during a transition from pressure-less matter domination to radiation domination, having in mind some models of reheating.They found that sudden transitions enhance the production of induced GWs [230].The main reason is what some of the authors later called the "poltergeist mechanism" [93] (see also Ref. [232] for an application using axions).Essentially what happens is the following: Density fluctuations grow during the matter dominated era and then, suddenly, everything (including those large density fluctuations) are converted into radiation.And, radiation wants to propagate.This creates big sound waves and a loud GW signal.It also gives, in general, a peaked GW signal since smaller scale fluctuations have more time to grow.It turns out that PBH evaporation after PBHs dominate the universe is a good example of an almost sudden transition [93], if one assumes a monochromatic mass function.
PBHs dominate the universe before evaporating if there is a large enough initial fraction of them.This is because their mean energy density redshifts initially as the volume (by the PBH number density conservation), i.e. as a −3 , while the energy density of radiation dominating the universe dilutes as a −4 .Thus, on one hand, PBHs dominate the universe for a/a f > β −1 .On the other hand, the time of evaporation, say t eva , is solely set by the initial PBH mass.And, if PBHs dominate, this time t eva also determines the Hubble parameter at evaporation, since H eva ∼ 1/t eva .From that, one determines that the temperature T eva of radiation filling the universe after PBH evaporation is given by (2.10) Requiring that evaporation occurs much later than domination leads to a lower bound on the initial fraction, namely We can also use that a successful Big Bang Nucleosynthesis (BBN) requires T eva > 4 MeV [233][234][235][236], to place an upper bound to the mass given by M PBH,f < 5 × 10 8 g . (2.12) There are, at least, two type of sources for induced GWs in the PBH reheating scenario.The first one are PBH number density fluctuations, which was first pointed out by Ref. [94].They come from the statistical fluctuations associated with the discreteness of PBHs.And, as PBH appear to a good approximation randomly and uniformly distributed in space, their fluctuations follow a Poisson distribution.The second source are primordial adiabatic curvature fluctuations [93], although one must extrapolate the results from CMB scales down to very small scales.
It is important to note that, in both cases, density fluctuations may enter the non-linear regime (i.e.δρ PBH /ρ PBH > 1), due to the growth of fluctuations in a matter dominated universe.Whether one stops the calculations at the onset of the non-linear regime or not, changes the final amplitude of the induced GWs.However, we note that curvature fluctuations, which are the source of induced GWs, remain always within the validity of perturbation theory.In this review we merely follow the approach of used in the corresponding previous works.In the first case, since the PBH density fluctuations are determined by the presence of PBH themselves, we estimate the GWs using the solutions of linear perturbation theory from the results of [95,96].But, for adiabatic fluctuations, we impose a conservative cut-off in the power spectrum at the comoving scale which becomes nonlinear at reheating as in Ref. [93], which we call k nl-cut .While the first approach may overestimate the GW signal, the second one underestimates it.A most accurate calculation must deal with the non-linear regime but likely requires the use of numerical simulations.There is also the fact that relaxing the monochromatic assumption for the PBH mass suppresses the induced GW signal [93], as PBH evaporation becomes more gradual the broader the mass function is.For recent hybrid N-body and lattice simulations of a gradual transition see Ref. [237].For estimations of the GWs from non-linear structure formation see Ref. [238] and for possible turbulence after evaporation see Ref. [239].See also Ref. [240] for GWs from the structure formation from Yukawa forces in the early universe.
a. GWs associated with PBH reheating I (isocurvature induced GWs): The type of initial conditions for Poisson PBH density fluctuations are called isocurvature fluctuations [241,242].This means that the initial PBH density fluctuations are compensated by equal in amplitude but opposite radiation density fluctuations.Basically, since PBHs originate from the collapse of density fluctuations in the radiation, a hole in the original radiation fluid is filled with a PBH.For a review on isocurvature induced GWs see Ref. [214] (see also Refs.[243] and [63] for GWs and PBHs from dark matter isocurvature and Refs.[244,245] for universal (isocurvature) Gravitational Waves associated with solitons).The spectrum of PBH density fluctuations grows as k 3 and it is largest at the scale corresponding to the mean inter-PBH separation, below which the PBH gas picture is no longer valid.We call the frequency associated with the mean inter-PBH separation f poisson .In terms of the PBH mass this frequency reads [95]  . (2.14) Then, the GW spectral density has a sharp cut-off for f > f possion and goes as f 11/3 for f < f possion .The factors 1/3 in the exponents come from the fact that curvature fluctuations on small scales are suppressed by the decay in the PBH mass, which goes as M PBH ≈ M PBH,f (1 − t/t eva ) 1/3 [93].
We show the GW spectral density in red in Fig. 2. From Eq. (2.14), we may use current BBN constraints [104,246] (see Ref. [247] for a recent review) to place an upper bound on the initial fraction of evaporated PBHs, which reads . (2.15) As far as we are aware, this is the only way to put an upper limit to the initial fraction of PBHs.We note that there are also induced GWs produced during the PBH dominated era [94].But, in the case of a monochromatic PBH mass function, the induced GWs after PBH evaporation constitute the largest contribution to the GW background.
b. GWs associated with PBH reheating II (adiabatic induced GWs): CMB observations measured an almost scale invariant spectrum of curvature fluctuations with amplitude A CMB R ∼ 2 × 10 −9 .If such a spectrum extends to very small scales and PBHs reheat the universe, then it yields an enhanced GW signal.However, if the primordial spectrum extends to arbitrary small scales, density fluctuations on certain scales become non-linear.To avoid such a regime Ref. [93] imposes a cut-off such that no fluctuation enters the non-linear regime up to evaporation.Borrowing the results from Ref. [93], this cut-off is given by Thus there are only density fluctuations with k < k nl-cut .In Eq. (2.16) P Φ (t eva ) is the spectrum of fluctuations of the gravitational potential Φ at evaporation and k eva = a eva H eva is the wavenumber that enters the horizon at evaporation.Translating the cut-off (2.16) into a frequency evaluated today yields 10 (2.19) The GW spectrum then quickly decays for f > f nl-cut and approximately goes as f 7 for f < f nl-cut until it smoothly transitions to an almost scale invariant plateau.The exponential dependence 10 To compare with Ref. [93] one should use that, in their notation, τeq,2/τeq,1 ≈ 100 . We also considered the limit where k nl-cut τeq,1 ≫ 1.And, although the highest GW production happens for k nl-cut τeq,1 ∼ 1, Eq. (2.19) still gives a good order of magnitude estimate. 11In the regime where k nl-cut ≫ keq, with keq being the wavenumber that enters the horizon at the first PBH-radiation equality, we find that f nl-cut /feva ≈ 10 −5 (Ω peak GW,nl−cut ) −3/7 .Then, since we have that the position of f nl-cut is determined by the amplitude of the GW spectrum and the mass of the PBHs.
that appears in both the peak frequency (2.17) and amplitude (2.19) comes from the logarithmic growth of matter fluctuations during the radiation dominated era [93].We show the resulting GW spectrum in green in Fig. 2. For more details on the calculations we refer the reader to Ref. [93].
(iii) GWs from PBH binaries.The last source of GWs associated with PBHs that we consider are GWs from PBH binaries.Most of these PBH binaries form in the early universe (unless the PBH fraction as dark matter is very small around f PBH < 10 −15 ) and they do so via a three body interaction [80].The two nearest PBHs fall towards each other but the third nearest PBH provides enough torque to the system to avoid a head-on collision.Then a very eccentric PBH binary is formed which eventually circularizes.For more details see the review [5] and references therein.
For nearby PBH binaries one may be able to resolve the GW waveform.An approximate estimate for the maximum GW frequency before merging, in the source frame, is given by [215] f max GW,binary ≈ 2.2 kHz This frequency is twice the frequency associated with the innermost stable circular orbit.If the binary is at cosmological distances, the frequency today changes by a factor 1/(1 + z) where z is the redshift.However, if PBH binaries are too far or too weak to be resolved, they contribute to the GW background.For the calculation of the GW background from PBH binaries with a monochromatic PBH mass function see, e.g., Ref. [84] and references therein (for the energy spectrum of binary black holes during the whole inspiral-merger-ringdown phase in the non-spinning limit see Refs.[251,252]).For an example of a broad PBH mass function see Ref. [85].Here we use an analytical estimate of the peak of the GW spectrum from Ref. [253] for a monochromatic PBH mass function, which reads The peak of the GW background from unresolved binaries comes from the nearest binaries and, therefore, the peak frequency is close to the maximum frequency of Eq. (2.21).For f > f max GW,binary the spectrum has a sharp cut-off and for f < f max GW,binary the GW spectrum decays as f 2/3 .We show the GW spectral density in purple in Fig. 2.This completes our list of estimates for the GWs associated with PBHs.

III. TESTABLE PBH MASS RANGE AT TERRESTRIAL GW DETECTORS
With the estimates derived in § II A we are ready to understand which PBH scenarios can be probed by terrestrial GW detectors.We note that our main aim is to recount the potential of terrestrial GW detectors to test PBH scenarios by providing order of magnitude estimates for future studies and searches.We will not dwell into details of how accurately future GW detectors may be able to probe PBH scenarios nor how well they may discern GWs from PBH scenarios from other GW sources.Instead, we refer the interested reader to Refs.[254][255][256][257][258]. As for terrestrial GW detectors, we consider a frequency range roughly from Hz to 10 kHz.When needed, we consider the peak sensitivity of Einstein Telescope which is most sensitive around 100 Hz with Ω GW,0 h 2 ∼ 10 −9 .To be more optimistic, we also consider the power-law integrated sensitivity curve [259] which gives a sensitivity to GW backgrounds of Ω GW,0 h 2 ∼ 10 −13 around 100 Hz, after accumulating several years of data.In the case when the GW peak frequency is higher, we find the parameter space for which the low frequency tail enters a power-law integrated sensitivity curve.
The main message of this and the next sections is summarized in Fig. 1.The estimates of the PBH mass ranges testable at terrestrial GW detectors are also summarized in more detail in Tab.III.At the end of this section, we also show an example of the GW spectral density together with sensitivity curves in Fig. 2. Now, we proceed to describe Fig. 1 based on our discussions in § II A. We start with GWs associated with evaporated PBHs and then turn to GWs associated  III.Summary of the mass ranges testable by terrestrial GWs detectors with the GW signals explained in points (i), (ii) and (iii).Detailed explanations on how to derive such estimates is given in points 1, 2, 3, 4 and 5.When needed, we take the power-law integrated sensitivity of ET, which is most sensitive at 100 Hz with Ω GW,0 h 2 ∼ 10 −13 .
with PBH binaries.

A. Probing evaporated PBHs
As evaporated PBHs cannot be detected directly, we only have the possibility to find their associated induced GW signal.This includes induced GWs associated with PBH formation and to PBH reheating, respectively discussed in points (i) and (ii).Looking at their respective estimates for the frequency, Eqs.(2.7), (2.13) and (2.17), we see that in order to enter the frequency range of terrestrial GW detectors, GWs associated with PBH formation require larger PBH masses than GWs associated with PBH reheating.Inputting some numbers in the estimates we find that: 1. PBHs with masses between 10 8 g − 10 16 g have the peak frequency (2.7) of GWs associated with PBH formation inside the frequency window.One may also be able to probe the amplitude the primordial spectrum almost down to A R ∼ 10 −5 in (2.8) (see Ref. [260] for a detailed analysis).
2. PBHs with masses between 10 6 g − 10 8 g lead to GWs associated with PBH formation detectable only through its low frequency tail.However, one needs A R > 10 −2 in (2.8).
3. PBHs with masses between 10 3 g − 10 8 g can be probed via isocurvature induced GWs associated with PBH reheating, since the peak frequency (2.13) is inside the observable frequency range.From Eqs. (2.13) and (2.14) we find that one could probe an initial PBH fraction from β > 6 × 10 −9 , where we used the peak sensitivity of the power-law integrated sensitivity curve of ET.
4. PBHs with masses between 10 2 g−10 4 g yield an adiabatic induced GWs from PBH reheating accessible to terrestrial GW detectors.Such mass range comes from using Eq.(2.17) and requiring the maximum amplitude possible in the GW spectrum (2.19). 13And once the GW amplitude is fixed, the peak frequency (2.17)only depends on the PBH mass.We note though that relaxing the non-linear cut-off would enhance the amplitude of the GWs and broaden the parameter space.
For Point 1, it is interesting to note that, although the fraction of PBHs in the mass range 10 9 g − 10 17 g is tightly constrained [6,[260][261][262], one may still be able to probe the induced GWs associated with their formation.Points 2, 3 and 4 show that we will be able to probe the existence of very small PBHs in the early universe, although they evaporated well before BBN.For Points 3 and 4, it would be intriguing to derive more accurate estimates using numerical simulations.
6. PBHs with masses between 10 −6 M ⊙ − 100 M ⊙ may yield a GW background signal from unresolved mergers.The mass range from 10 −6 M ⊙ − 10 −2 M ⊙ can be probed with the low frequency tail of the GW spectrum (2.22).In deriving this mass range we assumed f PBH ∼ 10 −2 which is consistent with current observations.A lower fraction of PBH as dark matter would yield a lower signal and a smaller mass range.
We also note that it may be possible to test the mass range from 10 −6 M ⊙ − 10 −2 M ⊙ via continuous GWs [255,263] (see also Ref. [264]), if the PBH fraction is not too small, around f PBH ∼ O(10 −2 ).For the mass ranges where there is overlap with astrophysical black holes, namely for M PBH,f > M ⊙ , one must carry out population analyses [87,258,[265][266][267], study the statistical nature of the GW background [257], or search for GW background anisotropies [268,269], in order to distinguish PBHs from astrophysical BHs.Interestingly, earth-to-solar mass PBHs may be tested by other complementary means such as microlensing or low frequency GWs, which we discuss in more detail in the next section.

IV. COMPLEMENTARY PROBES
This review focused on the role of terrestrial GW detectors in testing PBH scenarios.However, there are other promising ways to test PBHs and complement the information from terrestrial GW detectors.For instance, we may use Big Bang Nucleosynthesis (BBN) predictions and CMB observations, microlensing of electromagnetic waves and GW detectors in other frequency ranges.We list and describe them below.
(a) BBN & CMB.BBN predictions [246,247,276,277] as well as CMB observations [216,278] are sensitive to the presence of additional relativistic particles (sometimes also called dark radiation).Constraints from BBN and CMB are then usually parametrized with an effective number 2. Spectral density of GWs vs frequency in the range relevant for terrestrial GW detectors.We show one example for each GW signal associated with PBH discussed in this review (see points (i), (ii) and (iii)).In solid blue we show the GW spectrum of GWs associated with PBH formation.We considered M PBH,f ≈ 10 9 g and a log-normal primordial spectrum (2.8) with A R = 0.1 and ∆ = 0.1.We computed the GW spectrum using SIGWfast [270].In solid red we show the GW spectrum of isocurvature GWs associated with PBH reheating for M PBH,f ≈ 10 6 g and β ≈ 2 × 10 −8 .In solid green we show the adiabatic GWs associated wwith PBH reheating for M PBH,f ≈ 10 2 g and β ≈ 3 × 10 −5 , which was kindly provided by Keisuke Inomata.We believe that although the green line does not enter the observable window, the non-linear cut-off imposed by Ref. [93] largely underestimates the GW signal.In solid purple we show the low frequency tail of the GW background from PBH binaries with M PBH,f ≈ 10 −4 M ⊙ and f PBH ≈ 10 −2 [84].We also show the power-law integrated sensitivity curves [259] for Einstein Telescope (ET), Cosmic Explorer (CE), Voyager and LIGO A+ experiments (see Refs. [271][272][273][274] for the sensitivity curves).In light blue we plot the upper bounds from the LVK collaboration [275].The blue dashed line shows the current constraint from BBN [104,246,247].
of additional relativistic species, denoted by ∆N eff .Current limits from BBN [104] and CMB [170] respectively give ∆N eff ≲ 0.5 and ∆N eff ≲ 0.3.Future CMB experiments, such as CMB-S4 [279] might reach ∆N eff ≲ 0.02.The crucial point is that GWs with frequencies f ≳ 10 −10 Hz and f ≳ 10 −15 Hz may be considered as a dark radiation fluid, respectively, at the time of BBN and CMB [216].Thus, BBN and CMB provide an integrated constraint on the total spectral density of GWs above these frequencies, which roughly yields Ω GW,0 h 2 ≲ 10 −6 [216,278] (see also Ref. [280] sec.4.1 for a summary with a nice explanation).
Regarding the GWs associated with PBHs, BBN and CMB open the possibility to test a considerable amount of high frequency GWs.This is relevant for the GWs associated with PBHs reheating for 1 g < M PBH,f < 10 3 g as well as the GWs associated with PBH formation for 1 g < M PBH,f < 10 6 g, since they are not accessible to terrestrial GW detectors.Most interestingly, BBN and CMB might provide additional information on the PBH reheating scenario.GWs from Hawking evaporation of spinning PBHs is within the reach of future CMB-S4 experiments [102,104,105,281,282] (see also Ref. [283] for PBH evaporation with large extra dimensions).Such additional signatures might help in discerning the formation mechanism of PBHs for M PBH,f > 10 3 g [96].
(b) Microlensing.Electromagnetic waves, e.g. from stars, travel through the dark matter halos and, if dark matter is composed by PBHs, one expects a certain amount of lensing events depending on f PBH [284] (see also sec.3.1 of [5] for a detailed explanation of microlensing).Current bounds from the absence of microlensing events set a constraint of about f PBH ≲ O(10 −2 ) for a PBH mass range between 10 −10 M ⊙ − 10 M ⊙ [6]. 14However, most interesting are the microlensing candidate events reported by HSC [18] and OGLE [20] respectively with masses about 10 −8 M ⊙ and 10 −4 M ⊙ .We note that 10 −4 M ⊙ PBHs might also be linked to the reported PTA signal [126].For another example, see Ref. [286] where a broad PBH mass function may explain the reported Pulsar Timing Array (PTA) signal as well as PBHs as dark matter, and can be tested by microlensing observations.(c) PTAs.Another complementary window to terrestrial GW detectors are nHz GWs which may be probed by Pulsar Timing Arrays (PTAs).From the estimate of the peak frequency of GWs associated with PBH formation, Eq. (2.7), we see such GWs fall in the PTA range for O(10M ⊙ ) ≳ M PBH,f ≳ O(10 −3 M ⊙ ).If we consider the low frequency tail of the GW spectrum then the range might be extended down to M PBH,f ≳ 10 −9 M ⊙ with future SKA sensitivity [287].Most interesting though are the current results from PTA data [112][113][114][115][116][117][118][119][120][121], which seems to suggest M PBHs ∼ 10 −4 M ⊙ − 10 −3 M ⊙ , if interpreted as an induced GW signal associated with PBH formation.Such a mass range is very interesting as it has implications for the microlensing events reported by OGLE and it may have a detectable GW background from unresolved PBH binaries, as pointed out by Ref. [126].µHz GW detectors like µ-Ares [288] (see also Refs.[289][290][291]) would extend the range of PTAs and provide more evidence for the PBH interpretation and extend the testable PBH mass range.
(d) Space-based GW detectors.Future GW detectors such as LISA [292,293], Taiji [294], TianQin [295] and DECIGO [296,297], will bridge PTAs and µHz GW detectors with terrestrial GW detectors.This will offer the opportunity to test the low frequency tail of GW signals associated with PBHs that enter the terrestrial GW detector's window as well as, of course, to extend the testable PBH mass ranges.Details on the capabilities of LISA to test PBH scenarios can be found in the review by the LISA cosmology working group [204].
(e) MHz-GHz GW detectors.GWs associated with PBH formation and reheating of light PBHs as well as the mergers of planet-mass PBHs (see e.g.Refs.[298,299] for detailed studies) are the sources of high frequency GWs.Even GWs from Hawking evaporation in the PBH reheating scenario might be testable by MHz-GHz GW detectors [105,283].We note that although there are interesting events detected at current MHz GW detectors [300] and the frequency could be explained by the merger of planet-mass PBHs, it seems an extremely unlikely explanation given the current sensitivity [253] (see also Ref. [301]).Nevertheless, an improved sensitivity in future high frequency GW detectors will present an exciting window to further test PBH scenarios [302][303][304][305]. And, they will complement any possible signals seen at terrestrial GW detectors.
(f ) Lensing of GWs.Another interesting probe to PBH scenarios using GWs is the recently proposed lensing of GWs [306][307][308].Lensed GWs might be sensitive to dark matter halos substructure due to frequency dependent wave optic effects, which could probe PBHs as a fraction dark matter in the M ⊙ − 10 5 M ⊙ range [306].It would be interesting to investigate in which circumstances one might probe lighter PBH masses.

V. DISCUSSION AND CONCLUSIONS
Cosmic events that produced PBHs shook the spacetime, resulting in ripples that we see today as GWs.Such induced GWs are also produced if PBHs reheat the universe.The first type of induced GWs, that we called GWs associated with PBH formation, can test the presence of PBHs with masses 10 6 g−10 16 g in current and future terrestrial GW detectors, such as LVK, Einstein Telescope and Cosmic Explorer.The second type of induced GWs, here referred to as GWs associated with PBH reheating, offer means to probe even lighter PBHs with masses from 10 2 g to 10 8 g.Thus, terrestrial GW detectors have the potential to find GW signals associated with evaporated PBHs, otherwise unexplorable.
Interestingly, we may also be able to find hints of black hole remnants as dark matter [60], if PBH evaporation leaves Planck mass remnants behind as suggested by some quantum gravity theories [309][310][311][312][313] (see also Refs.[314][315][316][317]).If so, there is a unique initial PBH mass, that is M ∼ 5 × 10 5 g, that can reheat the universe and its remnants be the dark matter, with a sharp prediction for the frequency of induced GWs peaked at 100Hz [60].
In addition to evaporated PBHs, black hole binaries of long-lived PBHs, those with masses between 10 −6 M ⊙ − 10 2 M ⊙ , also produce GW signals at reach of terrestrial GW detectors.Most interesting for the PBH scenario is the possibility of finding evidence for planet-mass to sub-solar mass black holes, as their origin can only be primordial.And, there are potential candidates for such small PBHs in the LVK data [106][107][108][109], for sub-solar mass PBHs, and in the microlensing data of OGLE [19,20], for earth-mass PBHs (HSC [18] also reported one candidate even from a moon-mass object).One may search for such small PBHs with continuous GWs [255,263] and the GW background from unresolved binaries [84][85][86].The latter most often peaks at MHz-GHz frequencies, which may be accessible to high frequency GW detectors [302][303][304][305].
In a different direction, it would also interesting to fully explore how GW background anisotropies (see, e.g., Ref. [318][319][320]) might help in gathering evidence for the presence of PBHs [268,269,321,322] and, perhaps, even the presence of primordial non-Gaussianity [135,323,324].Moreover, as discussed in the context of MHz GW detectors by Ref. [301], high frequency GWs, e.g. from the mergers of planet-mass PBHs, might leave a GW memory trail at terrestrial GW detectors [325].
We end this review by collecting the various tables and illustrations that we hope will be helpful to the interested reader.Relevant estimates for the frequency and amplitude of the GWs are summarized in Tab.II.In Tab.III we provide details on the testable PBH mass range corresponding to the various GW signals discussed in § II A, points (i), (ii), (iii).We also illustrated the different possibilities and other complementary probes using GWs at other frequencies in Fig. 1.And in Fig. 2 we display the GW spectral shapes of the GW signals associated with PBHs for several different scnarios.
GW peak frequency only depends on the initial PBH mass.The amplitude of GWs at the peak frequency is estimated to be[95]

3 . 11 β 10 − 5 M
(2.17)    In this case the frequency depends on all the model parameters.The GW spectrum peaks at f nl-cut with amplitude 11 Ω peak GW,nl-cut ≈ 5 × 10 −

0. 1 MFIG. 1 .
FIG.1.Illustration of the PBH mass ranges that may be accessible to terrestrial GW detectors together with other complementary GW signals.Describing from left to right: First, we have PBHs evaporated before BBN (M PBH,f ∼ 10 2 g − 10 8 g) which can be detected from GWs associated with PBH reheating as well as high frequency GWs associated with formation and GWs from Hawking evaporation using CMB and BBN.Second, we have evaporated PBHs (M PBH,f ∼ 10 6 g − 10 16 g) that can be detected via GWs associated with PBH formation.Third and fourth, we respectively have GWs from unresolved (M PBH,f ∼ 10 −6 M ⊙ −100 M ⊙ ) and resolved (M PBH,f ∼ 0.1 M ⊙ − 10 2 M ⊙ ) PBH binaries.Their GWs associated with PBH formation are low frequency GWs within the range of PTAs.Other PBH mass ranges are accessible to space-based GW detectors, such as LISA and Taiji, PTAs as well as high frequency GW detectors.

M ≈ 10 2 gM ≈ 10 6 gM ≈ 10 9 g M ≈ 10 − 4 M
FIG.2.Spectral density of GWs vs frequency in the range relevant for terrestrial GW detectors.We show one example for each GW signal associated with PBH discussed in this review (see points (i), (ii) and (iii)).In solid blue we show the GW spectrum of GWs associated with PBH formation.We considered M PBH,f ≈ 10 9 g and a log-normal primordial spectrum (2.8) with A R = 0.1 and ∆ = 0.1.We computed the GW spectrum using SIGWfast[270].In solid red we show the GW spectrum of isocurvature GWs associated with PBH reheating for M PBH,f ≈ 10 6 g and β ≈ 2 × 10 −8 .In solid green we show the adiabatic GWs associated wwith PBH reheating for M PBH,f ≈ 10 2 g and β ≈ 3 × 10 −5 , which was kindly provided by Keisuke Inomata.We believe that although the green line does not enter the observable window, the non-linear cut-off imposed by Ref.[93] largely underestimates the GW signal.In solid purple we show the low frequency tail of the GW background from PBH binaries with M PBH,f ≈ 10 −4 M ⊙ and f PBH ≈ 10 −2[84].We also show the power-law integrated sensitivity curves[259] for Einstein Telescope (ET), Cosmic Explorer (CE), Voyager and LIGO A+ experiments (see Refs.[271][272][273][274] for the sensitivity curves).In light blue we plot the upper bounds from the LVK collaboration[275].The blue dashed line shows the current constraint from BBN[104,246,247].

TABLE II
. Summary of estimates for the frequency and amplitude of GW signatures associated with PBH discussed in points (i), (ii), (iii).A R is the amplitude of the primordial spectrum of curvature fluctuations, see Eq.(2.8