# When Primordial Black Holes from Sound Speed Resonance Meet a Stochastic Background of Gravitational Waves

###### Abstract

As potential candidates of dark matter, primordial black holes (PBHs) are within the core scopes of various astronomical observations. In light of the explosive development of gravitational wave (GW) and radio astronomy, we thoroughly analyze a stochastic background of cosmological GWs, induced by over large primordial density perturbations, with several spikes that was inspired by the sound speed resonance effect and can predict a particular pattern on the mass spectrum of PBHs. With a specific mechanicsm for PBHs formation, we for the first time perform the study of such induced GWs that originate from both the inflationary era and the radiation-dominated phase. We report that, besides the traditional process of generating GWs during the radiation-dominated phase, the contribution of the induced GWs in the sub-Hubble regime during inflation can become significant at critical frequency band because of a narrow resonance effect. All contributions sum together to yield a specific profile of the energy spectrum of GWs that can be of observable interest in forthcoming astronomical experiments. Our study shed light on the possible joint probe of PBHs via various observational windows of multi-messenger astronomy, including the search for electromagnetic effects with astronomical telescopes and the stochastic background of relic GWs with GW instruments.

###### pacs:

98.80.Cq, 11.25.Tq, 74.20.-z, 04.50.Gh## I Introduction

Primordial black holes (PBHs) are hypothetical objects predicted by many fundamental theories to form soon after the big bang, and hence the search for PBHs offer an inspiring possibility to probe physics in the early Universe Zel’dovich and Novikov (1966); Hawking (1971); Carr and Hawking (1974). Since they could be a possible candidate for dark matter (DM), the studies on cosmological implications of PBHs are crucial for cosmologists Ivanov et al. (1994); Carr et al. (2016); Gaggero et al. (2017); Inomata et al. (2017a); Georg and Watson (2017); Fuller et al. (2017); Kovetz (2017); Cai et al. (2018a); Nakama and Wang (2019); Ballesteros et al. (2018); Carr and Kuhnel (2018); Deng et al. (2018); Dalianis (2018). Recently, it was pointed out that the PBHs might be responsible for some gravitational wave (GW) events Bird et al. (2016); Clesse and García-Bellido (2017); Sasaki et al. (2016); Nakamura et al. (1997), which were detected by GW instruments such as the LIGO Abbott et al. (2016). This observational possibility motivates many theoretical mechanisms generating PBHs, which often require a power spectrum of primordial density perturbations to be suitably large on certain scales that are associated with a particularly tuned background dynamics of the quantum fields in the early Universe (e.g. see Garcia-Bellido et al. (1996); Garcia-Bellido and Ruiz Morales (2017); Domcke et al. (2017); Kannike et al. (2017); Carr et al. (2017); Ballesteros and Taoso (2018); Hertzberg and Yamada (2018); Franciolini et al. (2018); Kohri and Terada (2018a); Özsoy et al. (2018); Biagetti et al. (2018) for studies within inflation, see Chen et al. (2017); Quintin and Brandenberger (2016) for discussions within bounce, and see Khlopov (2010); Sasaki et al. (2018) for recent comprehensive reviews).

Primordial density perturbations, that seeded the large-scale structure (LSS) of the Universe, are usually thought to arise from quantum fluctuations during a dramatic phase of expansion at early times, as described by inflationary cosmology, from which a nearly scale-invariant power spectrum with a standard dispersion relation is obtained Mukhanov et al. (1992). This was confirmed by various cosmological measurements such as the cosmic microwave background (CMB) radiation Ade et al. (2015); Akrami et al. (2018) and LSS surveys at extremely high precision. It is interesting to note that, although density and tensor perturbations evolve independently through the early Universe at linear level, they couple with each other nonlinearly and hence can induce either the non-Gaussianities or stochastic background of relic GWs Ananda et al. (2007); Baumann et al. (2007). Once if these density perturbations could form PBHs, it becomes possible to constrain primordial non-Gaussianities with PBHs Bullock and Primack (1997); Ivanov (1998); Saito et al. (2008); Seery and Hidalgo (2006); Byrnes et al. (2012); Young and Byrnes (2013). Moreover, it is intriguing to search for PBHs via the measurements of the stochastic GW background induced by over large primordial density perturbations Saito and Yokoyama (2009, 2010); Bugaev and Klimai (2011); Garcia-Bellido et al. (2016); Inomata et al. (2017b); Garcia-Bellido et al. (2017); Kohri and Terada (2018b); Bartolo et al. (2018a); Cai et al. (2018b); Bartolo et al. (2018b); Wang and Zhang (2018); Unal (2018); Inomata and Nakama (2018); Clesse et al. (2018). The connection between PBHs and the induced GWs has been extensively studied in the literature, which mainly focused on the GW generation during the radiation-dominated phase when PBHs were formed. So far, the GW generation during the inflationary phase remain unclear, which is one of the main subjects of the present study.

Recently, a novel mechanism for PBHs formation, called sound speed resonance (SSR), was proposed in Cai et al. (2018a), where an oscillating sound speed square leads to the non-perturbative parametric amplification of certain perturbation modes during inflation. As a result, the power spectrum of the primordial density perturbations has a narrow major peak on small scales, while remains nearly scale-invariant on large scales as predicted by inflationary cosmology. Note that several minor peaks of the power spectrum of the primordial density perturbations on smaller scales are also predicted by this novel mechanism. In Cai et al. (2018a) it was found that the formation of PBHs caused by the resulting peaks in SSR can be very efficient, which could be testable in the future observational experiments.

Moreover, the enhanced primordial density perturbations are expected to induce large GWs signals according to the second-order cosmological perturbation theory. Motivated by the aforementioned reasons, in the present paper we turn to study the second-order GWs caused by the primordial density perturbations with peaks in the SSR scenario. Making use of this new mechanism, we perform a full analysis of the GWs signal, that evolves through the inflationary era until the present Universe.

First, we analyze the GWs induced by the enhanced primordial density perturbations when the relevant modes re-enter the Hubble horizon during the radiation-dominated era. Afterwards, we study GWs induced by the modes of primordial density perturbations at the super-Hubble scales during the inflationary era, which is often omitted in other works. Our calculation reveals that this part of contribution is usually suppressed by the slow-roll parameter when compared with that of the radiation-dominated era, but it might be important if the slow-roll condition could be violated for a short while during inflation. Finally, we compute the GWs induced by quantum fluctuations that remain inside the Hubble horizon during inflation, and we find that this sub-Hubble contribution can become very significant at critical frequency band due to the narrow resonance effect. We interestingly observe that, although the spectrum is damping out at small scales due to a blue tilt, there exists a narrow window where the induced GWs can be resonantly enhanced, and the resulting amplitude can be sizable when compared with the one derived in the radiation-dominated phase. Therefore, this work provides for the first time a comprehensive study on stochastic GWs background induced during both the inflationary and the radiation-dominated eras. And the resulting signal provides a new target for various future GW experiments, which also serves as an independent window of probing PBHs in the Universe.

The article is organized as follows. In Section II, we put forward a novel parametrization form of the power spectrum of primordial density perturbations that allows for several spikes and discuss the associated realization from the perspective of the SSR mechanism. In Section III, we work out the GWs induced by primordial density perturbations during radiation-dominated phase, and by the super-Hubble modes and the sub-Hubble modes during inflation, respectively. In Section IV, we derive the energy spectra of the induced stochastic GWs associated with the PBH formation and perform a comparison with the observational ability of the present and forthcoming GW experiments. Finally, we conclude with a discussion in Section V. The detailed calculation of induced GWs is presented in Appendices. Throughout the article, we adopt the natural units and the reduced Planck mass is defined as .

## Ii Sound Speed Resonance and Power Spectrum with Peaks

To generate PBHs within inflationary cosmology, the key point is to amplify the amplitude of primordial density perturbations for certain ranges of modes. For most of theoretical mechanisms studied in the literature, it requires a manifest enhancement around a unique comoving wavenumber, and accordingly, the mass spectrum of PBHs displays a single peak around a critical mass scale. However, it was recently pointed out in Carr and Kuhnel (2018) that the PBHs are likely to have an extended mass spectrum, in particular with multiple peaks, which has crucial implications for the interpretation of the observational constraints. This novel phenomenon happens to be realized also in Cai et al. (2018a) in terms of the SSR mechanism. We note that, for the sound speed, the first peak of the resonant power spectrum, which corresponds to the lowest comoving wavenumber, makes the larger contribution to the PBH mass spectrum, but the rest would affect the whole profile, in particular the tail of the mass spectrum.

The causal mechanism of generating power spectrum suggests that primordial density perturbations initially emerge inside the Hubble radius, and then exit in the primordial epoch, and eventually re-enter at late times. One often uses a gauge-invariant variable , the curvature perturbation in comoving gauge, to depict the primordial inhomogeneities. For convenience, one can introduce a canonical variable , where Armendariz-Picon et al. (1999); Garriga and Mukhanov (1999). Note that, is often regarded as the Hubble slow-roll parameter, is the Hubble parameter and is the sound speed parameter of the primordial Universe. The evolution of one Fourier mode for this variable satisfies: , where the prime represents for the derivative w.r.t. the conformal time .

To produce PBHs within inflationary cosmology, it requires a dramatic amplification of the primordial curvature perturbations for certain scales. In Cai et al. (2018a), a novel mechanism was proposed by introducing an oscillating correction to the sound speed parameter. In particular a parametric amplification of curvature perturbations caused by resonance with oscillations in the sound speed parameter, which provides an efficient way to enhance the primordial power spectrum around the astrophysical scales where PBHs, could account for DM in the current experimental bounds. Such an oscillation correction could arise when inflation models are embodied in UV-complete theories, such as D-brane dynamics in string theory Silverstein and Tong (2004); Alishahiha et al. (2004), or, by integrating out heavy modes from the effective field theory viewpoint Achucarro et al. (2011, 2012), or purely from a phenomenological construction Cai and Zhang (2009); Cai and Xia (2009); Cai and Xue (2009).

Specifically, the sound speed parameter can be parametrized as: with , where is a small dimensionless quantity that measures the oscillation amplitude and is the oscillation frequency. Note that is required such that is positively definite, and, the oscillation begins at , where needs to be deep inside the Hubble radius with . Moreover, we set before and assume that it can transit to oscillation smoothly for simplicity. In this mechanism, the perturbation equation can be re-expressed in form of a Mathieu equation Cai et al. (2018a), which then gives rise to a narrow parametric resonance, i.e., the perturbation modes in the neighborhood of the characteristic scale can be exponentially enhanced, while the power spectrum on large scales remains nearly scale invariant as predicted by the standard inflationary cosmology. However, it is interesting to observe that, this mechanism also predicts several secondary peaks on the scales with integer times of . Accordingly, by taking into account these peaks into the power spectrum, one can parametrize its form as follows,

(1) |

where is the amplitude of the power spectrum predicted by the conventional inflationary paradigm, is the Hubble slow-roll parameter as mentioned above, is the spectral index at the pivot scale Akrami et al. (2018). The resonant enhancements are characterized by the delta functions inside the square brackets in the second line of Eq. (II), where the amplitude of -th peak relative to the first one is quantified by a dimensionless parameter . Since we work in the perturbative regime, the height of the peaks in should be no more than unity. This condition constrains the upper bounds for the amplitudes of the induced GWs generating from the inflationary era and the radiation-dominated era.

To illustrate the technique of calculating GWs induced by the process of PBH production with multiple spikes, throughout the whole analyses we will only take into account the second and third peaks that are located at and respectively. This is because the amplitudes of other peaks at higher order are exponentially suppressed by a factor of about in the SSR mechanism and then soon be out of observable interest. According to the cosmological perturbation theory, different -modes of linear density perturbations can be mixed with each other nonlinearly and this mixing can play the role of generating tensor perturbations at second order Ananda et al. (2007); Baumann et al. (2007).

In the literature, there were extensive studies on the GWs induced by a single-peak pattern of the power spectrum of primordial density perturbations with the scalar modes re-enter the Hubble radius during the post-inflationary phase Saito and Yokoyama (2009, 2010); Bugaev and Klimai (2011); Garcia-Bellido et al. (2016); Inomata et al. (2017b); Garcia-Bellido et al. (2017); Kohri and Terada (2018b); Bartolo et al. (2018a); Cai et al. (2018b); Bartolo et al. (2018b); Wang and Zhang (2018); Unal (2018); Inomata and Nakama (2018); Clesse et al. (2018). These works have successfully demonstrated that a possible measurement of the stochastic background of GWs in future astronomical experiments could provide a powerful window to search for PBHs or to constrain the parameter space. However, it is crucial to notice that, this observational attempt requires a precise quantification of the profiles of energy spectra of the induced GWs, which also takes into account of the enhancement effects by the specific mechanisms of PBH generations in the primordial era.

To address the aforementioned issues, in the present article we perform a much more comprehensive analysis on how the GWs were induced by a resonantly enhanced power spectrum that began from the sub-Hubble regime of the inflationary era, then evolved to the super-Hubble regime during inflation, and finally re-entered the Hubble horizon during the radiation-dominated phase. The PBHs form in the radiation-dominated phase, but the GWs can be continuously produced throughout the whole evolution as shall be seen in the following section.

## Iii Gravitational waves generated in a primordial Universe

In this section, we present a complete analysis of the induced GWs throughout the whole evolution of a primordial Universe. In order to make the comparison with the pioneer works in the literature, we would like to show the results backward in time. We begin by discussing the induced GWs when the scalar modes re-enter the Hubble horizon during the radiation-dominated era.

### iii.1 Radiation-dominated era

In the Newtonian gauge, the line element in the perturbed metric is expressed as Ananda et al. (2007); Inomata et al. (2017b); Kohri and Terada (2018b)

(2) |

where is the conformal time, is the Bardeen potential, and is the tensor perturbation up to second order. The equation of motion (EoM) for the induced GWs in the Fourier space can be written as follows,

(3) |

where is the comoving Hubble parameter, denote two polarization modes of the GWs, and the source term during the radiation-dominated era is given by Kohri and Terada (2018b); Baumann et al. (2007); Ananda et al. (2007); Bartolo et al. (2018b)

(4) |

where the projections are defined as Saito and Yokoyama (2010). They take the form of for , and, for , where , is the angle between the wave vector (of the induced GWs) and (of the source), and is the azimuth angle of . The basic formulae above are presented in Appendix A in more detail. Accordingly, the EoMs for the induced GWs (3) are actually not equivalent for different polarizations, and hence we keep the indices of polarizations in the expressions. During the radiation-dominated era, the Bardeen potential in the source term (III.1) is related to the primordial comoving curvature perturbation via the relation , and the transfer function can be solved by the EoM for the Bardeen potential (see the Appendix B.1 for details) Ananda et al. (2007); Baumann et al. (2007).

Note that, while the evolution of the induced GWs during the radiation-dominated phase is relevant with the initial conditions that were inherited from inflationary phase, their generations mainly rely on the source term. Actually, the contribution induced from the source term would dominate over that from the cosmological evolution very soon in the radiation-dominated phase. The special solution of the induced GWs (3) can be determined by the Green function method. From the EoM (3) and the source term (III.1), one can see that the two-point correlation function of the induced GWs can be roughly estimated as the square of the two-point correlation function of Baumann et al. (2007); Ananda et al. (2007). After some lengthy calculations, one can derive the power spectrum of the induced GWs during the radiation-dominated phase as follows Kohri and Terada (2018b); Baumann et al. (2007); Ananda et al. (2007); Bartolo et al. (2018b),

(5) |

where we have introduced the variables and . Moreover, the function is given by

(6) |

where has been introduced. The expressions for the functions and are provided in the Appendix B.1 for detailed information. The formalism (III.1) is the general expression to calculate the power spectrum of the GWs induced by the primordial curvature perturbation when the associated Fourier modes of re-enter the Hubble horizon during the radiation-dominated era. The same modes of curvature perturbations can also induce over large density fluctuations that eventually could form PBHs after re-entering the Hubble horizon during the radiation-dominated phase.

We comment that, a significant property of the GWs induced by the multi-spike power spectrum of the primordial curvature perturbations is that the corresponding wave band is wider than that of the single-peak case Cai et al. (2019). This can lead to differences in the parametrization of the energy spectrum for the induced GWs, which shall be discussed in detail in Sec. IV. Note that, the average amplitude of the induced GWs is not altered by the multi-spike case, as this amplitude mainly rely on the first peak that in the SSR mechanism is always the dominant one.

### iii.2 Inflationary era

In this part, we investigate the induced GWs during the inflationary era, which is often omitted in previous works. In the SSR mechanism, the curvature perturbations are amplified deep inside the Hubble horizon at the beginning time while the sound speed start oscillating. And then, the curvature perturbations are exponentially enhanced due to the narrow resonance effect until the relevant modes exit the horizon, and they are frozen at the super-horizon scales. In principle, the GWs could be induced in both the sub-Hubble regime and the super-Hubble regime. Thus, the total contributions to the induced GWs during the inflationary era consist of these two regimes. First of all, we calculate the power spectrum for the GWs induced by the super-Hubble modes of the primordial curvature perturbations with the multiple peaks. Afterwards, we also calculate the power spectrum of GWs induced by the sub-Hubble modes of the primordial curvature perturbations which remain the quantum nature.

#### iii.2.1 Super-Hubble modes

During the inflationary era, the perturbed inflaton can provide the anisotropic stress which then sources the GWs nonlinearly. In this case, the source term in (3) becomes

(7) |

where is the perturbed inflaton in the Fourier space. We work in the spatially flat gauge, where the field fluctuation is related to the curvature perturbation via at the super-Hubble scales. After the lengthy calculations similar to (III.1), we obtain the power spectrum for the induced GWs as follows

(8) |

The function here is

(9) |

where , and is the characteristic frequency in the SSR mechanism. Comparing the formula (III.2.1) to (III.1), we can see that the difference between two power spectra of the induced GWs is from the functions and , since the source term during inflationary era (7) and radiation-dominated era (III.1) are different. We notice that the function involves the slow-roll parameter , and hence, the magnitude of the power spectrum (III.2.1) is sensitive to the value of . Consequently, the energy spectrum of the GWs induced by the super-Hubble modes is about suppressed comparing with the GWs induced during the radiation-dominated era, see Sec. IV. The explicit expression for induced by the single-peak pattern and multi-peak pattern of can be found in (B.2) and (55) of the Appendix B.2. According to our analysis, the power spectrum for the induced GWs in the super-Hubble regime is frozen at the end of inflation, and is shown in the Fig. 1, where is the conformal time at the end of inflation, i.e. . It is easy to see that the peak of the power spectrum (III.2.1) is located at , where the amplitude arrives at , and we use the values of parameters as , , , and which is in agreement with the latest CMB observations Akrami et al. (2018). However, the effects of multiple peaks are no longer manifest in the tail of the power spectrum due to the suppression of the higher peaks by a factor of about .

#### iii.2.2 Sub-Hubble modes

In this part, we continue to study the induced GWs from the quantum fluctuations of the perturbed inflaton inside Hubble horizon during the inflationary era. Due to the narrow resonance effect in the SSR mechanism, the sub-Hubble modes of in the neighborhood of the characteristic scale can be exponentially amplified, leading to the first major peak in the curvature power spectrum (II). Since the sub-Hubble modes are time-dependent, it is not easy to calculate the four-point correlator of the as the previous calculations. Moreover, the solution of the modes in the neighborhood of the characteristic scale is a complicated combination of Mathieu functions Cai et al. (2018a). Therefore, we use a semi-analytical method to compute the power spectrum for the GWs induced by the sub-Hubble modes, instead of the previous analytical approach used for the super-Hubble regime.

The power spectrum for the GWs induced by the sub-Hubble modes at time can be calculated by the formula Biagetti et al. (2013),

(10) | ||||

where is the Green’s function for the induced GWs, and is the mode function of the quantum fluctuation . The angle is spanned by the wave vector and . Equation (10) consists of two integrals, the phase space integral and the time integral. For the phase space integral, we use the thin-ring approximation (see the Appendix B.3), while the integral over interacting time is performed in a numerical way.

In the SSR mechanism, the full form of the mode function cannot be obtained analytically. Instead we use the sub-horizon approximation and neglect the Hubble friction term in the EoM for inflaton. After the approximation, the equation is of the standard Mathieu form and can be solved analytically in terms of Mathieu functions. The validity of the approximation depends on the physical wavelength of the considered mode. If its physical wavelength is small compared to the Hubble horizon, the approximation is expected to work well. Since here we are considering the sub-Hubble contributions, we only need to track down to its horizon-crossing time. Henceforth, this approximation should be valid for the most region of integration.

To better organize the calculation, we introduce a dummy variable . Thus the time integral can be rewritten as

(11) | ||||

where the function is

(12) |

As a shorthand notation, we denote the Matheiu sine and cosine funtions as , , and , .

Therefore the contribution from sub-horizon modes is written as

(13) |

where

(14) | ||||

Thus we arrive at the power spectrum induced by the first major peak in the power spectrum (II). The same treatments are applied to the analysis of the second and third resonance peaks in (II). Note that the formula (10) and the thin-ring approximation can also be applied to the calculations of the power spectrum for the GWs induced by the super-Hubble modes at the end of inflation, by changing the expressions for the mode functions and the interval of the integral to (, ). These two approaches give the same results, which confirm our previous calculations in (III.2.1) and (III.2.1).

Note that the power spectrum for GWs induced by the sub-Hubble modes is controlled by . This is the exponential of e-folding numbers between the triggering of resonance and the horizon-exit of the characteristic inflaton mode. Larger gives higher GW spectrum. However, this amplification cannot grow indefinitely. In order to keep the validity of our formalisms, we require and . The former bound is considered in Cai et al. (2018a) and gives

(15) |

Notice that and do not explicitly enter the expression for (13), but they do control the boundary of parameter region through (15).

In Fig. 1, we depict the power spectra in the super-Hubble regime and the sub-Hubble regime for a viable set of parameters indicated in the caption. We see that the average amplitude of the GWs induced by the sub-Hubble modes can be about larger than that of the super-Hubble modes. This gives us ample reasons to think that no significant corrections can arise from the evolution after (i.e. the super-Hubble regime). Physically speaking, it is also reasonable since the induced GWs possess a comparable wavelength as the characteristic mode. Therefore when the curvature perturbation of the characteristic mode is frozen outside the horizon, the induced GWs freeze as well. This conclusion is also confirmed by the numerical result, and thus the induced GWs from inflation is dominated by the contribution from the sub-Hubble regime, i.e. .

Furthermore, the induced GWs from the inflationary era are large enough to be of observational interest, and could be comparable to the induced GWs from the radiation-dominated era. We shall give more detailed discussion on this point in Sec. IV. Here we notice that the peak of the induced GWs by the sub-Hubble modes is located around , which in our specific example takes . Thus the location of the peak differs from the one of the induced GWs by the super-Hubble modes, as well as the one from the radiation-dominated era, both of which are located in the neighborhood of . This shift of the peak is due to the interplay of phase integral and time integral in (10). First, the mode functions take the complicated forms in terms of Mathieu functions at the sub-Hubble scales, the integral (11) of mode functions would have significant effects on the location of the peak of the power spectrum . Secondly, at small , the thin ring deforms to a shell due to the geometric effects, see the Appendix B.3 for details. Another different feature comparing with the power spectrum induced by the super-Hubble modes is that the tail of power spectrum in the sub-Hubble regime has sharp peaks.

We comment that the large GWs induced by the sub-Hubble modes arise from two effects. First, the expansion of the Universe from the beginning of the sound speed oscillation to the Hubble crossing of the mode: , where is the -folding number for this period of inflation. The GWs induced by the super-Hubble modes are actually contributed from the integral in (10), i.e., from the Hubble crossing of the mode to the end of inflation : . In the SSR mechanism, is quite larger than 1. Another effect is during the late stage of the resonance, the mode function oscillates in a trigonometric way, giving rise to a non-vanishing central value of which accumulates in the time integral in (10).

In Fig. 2, we choose three sets of independent parameters , and to study the dependence of inflationary power spectra on them. As expected, the power spectra are diminished by decreasing any of them while the shapes of the curves remain barely changed. This is reasonable since controls the duration of resonance and controls the rate of amplification. A decrease in either of them should cause to drop down. The slow-roll parameter appears in the overall normalization and has a simple quadratic dependence.

## Iv Results

In the previous sections, we provided theoretical analyses of cosmological GWs induced by primordial curvature perturbations with multiple spikes (II) from the inflationary era to the radiation-dominated era. We derived the power spectra for the induced GWs, which are listed in Eqs. (III.1), (III.2.1) and (10), respectively. Now we turn to forecast the observational implications on these induced GWs with the forthcoming GW experiments, e.g. the Laser Interferometer Space Antenna (LISA) Audley et al. (2017), Big-Bang Observer (BBO) Harry et al. (2006), Deci-hertz Interferometer Gravitational Wave Observer (DECIGO) Kawamura et al. (2011), and TianQin Luo et al. (2016). Moreover, the induced GWs at low frequency band within the range of may be accessible by the planned pulsar timing array experiments, e.g. the Square Kilometre Array (SKA) Santander-Vela et al. (2018), International Pulsar Timing Array (IPTA) Manchester (2013). To combine the observational windows of radio surveys and GW interferometers, the search for PBHs by virtue of probing the induced GWs is becoming promising in the era of multi-messenger astronomy.

When the induced GWs associated with the PBHs formation evolve into the present along with the expansion of our Universe, they become a stochastic GW background which can be characterized by their energy spectra Maggiore (2000); Allen and Romano (1999). Usually this is defined as energy density of the GWs per unit logarithmic frequency, and is the reduced dimensionless Hubble parameter at the present time . According to the definition of the effective energy for the GWs and the energy spectrum (see the Appendix A), we get the relation between the energy spectrum and the power spectrum for the GWs (see (35)) when the modes of GWs are well inside the Hubble horizon. For the induced GWs from the radiation-dominated era as shown in Eq. (III.1), the energy spectrum observed today is estimated as Saito and Yokoyama (2010); Bartolo et al. (2018a):

(16) |

where is some time near the end of the radiation-dominated era and is the present radiation energy density fraction. Since the energy density of GWs scales as radiation along with the cosmic expansion, the energy spectrum for the induced GWs do not dilute during the radiation-dominated era.

For the induced GWs during the inflationary era, which are given by Eqs. (III.2.1) and (10), the present energy spectrum can be approximately written as Boyle and Steinhardt (2008); Zhao and Zhang (2006)

(17) |

for the frequency Hz.

The present energy spectra of the induced GWs combined with the sensitivity curves of LISA, SKA and IPTA are shown in Fig. 3 and Fig. 4 respectively. Besides, the energy density of PBHs to DM ratio is another important phenomenological parameter. In each figure, the constraints from the corresponding electromagnatic window on spectra produced by the SSR mechanism are also shown. Here the relation between PBH mass and the frequency of the induced GWs follows an inverse-square law, namely . For the sensitivity of LISA, we have chosen the scale Mpc associated a PBH of mass Bartolo et al. (2018a), and the scale near the sensitivities of SKA and IPTA corresponds to the PBHs with mass about .

Both in Fig. 3 and Fig. 4, we can see that the present energy spectra from the inflationary era are dominated by the sub-Hubble contribution, which are comparable to the present energy spectrum from the radiation-dominated era when we choose the values of parameters as follows: , , and . In Fig. 3, the energy spectra of the induced GWs, both from the inflationary era and the radiation-dominated era, exceed the sensitivity of LISA, and their peaks are located in the sensitive region of LISA. Moreover, the frequency of the peak of energy spectrum contributed from the inflationary era is around lower than that of radiation-dominated era, so that it’s hopeful to distinguish these two signals in the future LISA experiments. In Fig. 4, the induced GWs from the inflationary era and the radiation-dominated era both exceed the sensitivities of SKA and IPTA, nevertheless the signal from the inflationary era is smaller than the signal from the radiation-dominated era. However, the above conclusions are based on the choice of the characteristic scale in sound speed oscillation. If we choose an appropriate scale , the signal of the induced GWs from the inflationary era becomes detectable by LISA, SKA and IPTA. Hence, through the GW observations we can extract constraints on parameters , and in the SSR mechanism.

The shapes of energy spectra of stochastic GW background are crucial for the GW observations.
As an approximation, we introduce the parametrization of the energy spectrum in the power-law form^{1}^{1}1The parametrization of stochastic GW background is found to be powerful in probing the cosmic history such as in Kuroyanagi et al. (2018).:

(18) |

where represents the amplitude of energy spectrum at the critical frequency , i.e. . and are the indices of the spectrum. For the induced GWs from the inflationary era, the parameters are given by

(19) |

For the induced GWs from the radiation-dominated era, we have

(20) |

Note that due to the term appears in power spectrum (III.1), the amplitudes of the parameterized energy spectra from the radiation-dominated era would change in different frequency ranges. As we have mentioned before, for the induced GWs from the sub-Hubble regime, the term does not explicitly enter the expression for (13), hence the amplitudes of the parameterized energy spectra are same in both frequency ranges of LISA and SKA (IPTA).

## V Conclusion

In this article we perform a comprehensive analysis of the stochastic background of GWs induced nonlinearly by over large primordial density perturbations that are accompanied with a process of PBH production. We report for the first time that the induced GWs can be resonantly enhanced within the sub-Hubble regime during inflation and hence make significant contribution to the energy spectra that are of observable interest in the forthcoming GW experiments. Our study also confirms that the contribution of induced GWs in the super-Hubble regime during inflation is secondary due to the nature of slow-roll suppression. Accordingly, by summing up all contributions the energy spectra of the induced GWs display a unique double-peak pattern that is innovative when compared with other works. To develop the technique of computation, we proposed a novel parametrization of the power spectrum of primordial density perturbations with several spikes as inspired by the SSR mechanism. It is acknowledged that the first two peaks of power spectrum would make the most important contribution to induce GWs nonlinearly. While the rest peaks could continue to contribute at a secondary level, there is a steep tail on the profile of energy spectra at the high frequency band that damps out soon. This is manifestly different from the single peak case where the energy spectra were almost cutoff at certain frequencies. In order to precisely characterize the energy spectra of the induced GWs from SSR mechanism, we put forward an envelope parametrization of their profiles which are expected to be measured in the future astronomical observations.

Note that, depending on the mass scales of the PBHs, the characteristic peaks of energy spectra of the induced GWs may be located within the frequency band of which is sensitive to the satellite projects of GW astronomy, or even as low as the frequency band of which then is most relevant with radio-astronomy. It is interesting to observe that, even for those PBHs with extremely light masses that have already been evaporated in the history of our Universe, they may leave a relic of induced GWs at high frequency band, namely within the LIGO range or even the kHz regime. This remarkable feature implies that, the probe of the energy spectra of induced GWs in multi-messenger astronomy has crucial implications for the search of PBHs at almost full frequency bands, even if those black holes may already disappear due to the mass loss via Hawking radiation.

Moreover, by comparing the observational abilities of GW astronomy with those of traditional telescopes upon the PBH mass spectrum, our results show explicitly that the observational window of GW instruments shall be more promising. To be specifically, one could constrain the parameter space of the SSR mechanism by testing the energy spectrum of induced GWs within the scope of LISA, while the dominant peak of the corresponding mass spectrum might be far from the sensitivities of traditional telescopes. This fact, from another perspective, has well illustrated that the development of multi-messenger astronomy shall become more and more promising in particular on the power of future GW detections.

In the end, we wish to highlight the implications of the present analysis on future studies from several perspectives. First of all, although our analysis is based on the SSR mechanism, from the perspective of methodology, the techniques developed in the present article can be easily extended into other scenarios with extended mass spectrum for PBHs. Moreover our results indicate that, for the study of induced GWs, the traditional approach may be incomplete, since it only focused on the radiation dominated stage, but there can also be significant contributions from the sub-Hubble regime during inflation. Therefore it is worth checking the GWs signals from other PBH generation mechanisms, and also taking into account the contribution from inflation stage. Phenomenologically, an important lesson from our study is that the probe of the information about the very early Universe is no longer limited by the traditional CMB and LSS surveys, but also include other astronomical telescopes at much smaller scales as well as a brand new window of GW experiments. On the other hand, in the era of multi-messenger astronomy, the search for PBHs are becoming more and more likely, making for a more and more serious DM candidate, which would inspire appropriate designs for the forthcoming experiments. In particular, a possible measurement of energy spectra of stochastic GW background with high precision could shed light on the nature of black holes existing in our Universe.

Acknowledgments.–
We are grateful to A. Achúcarro, R. Brandenberger, B. Carr, C. Germani, S. Pi, M. Sasaki, and E. Saridakis
for stimulating discussions and valuable comments.
This work is supported in part by the National Youth Thousand Talents Program of China, by the NSFC (Nos. 11722327, 11653002, 11421303, J1310021), by the CAST Young Elite Scientists Sponsorship (2016QNRC001), and by the Fundamental Research Funds for Central Universities.
DGW is supported by a de Sitter Fellowship of the Netherlands Organization for Scientific Research (NWO).
All numerics are operated on the computer clusters LINDA & JUDY in the particle cosmology group at USTC.

## Appendix A Basic Formulae for the Induced Gravitational Waves

We review the basic formulae for the induced GWs Maggiore (2007); Boyle and Steinhardt (2008); Baumann et al. (2007); Ananda et al. (2007). The EoM for the GWs induced by a source can be written as

(21) |

where is the comoving Hubble parameter, the prime denotes the derivative w.r.t. the conformal time and the projector operator is used to extract the transverse and traceless part of the source, i.e.

(22) |

Note that repeated Latin indices refer to summation, and denote the polarizations of the GWs. The polarization tensor is expressed in terms of a pair of polarization vectors and , both of which are orthogonal to the wave vector ,

(23) |

Then one can write the EoM (21) for the GWs in Fourier space as

(24) |

where

(25) |

Here we take the Fourier transformations as

(26) |

and

(27) |

The particular solution of the EoM (24) is given by the Green function method,

(28) |

where the Green’s function satisfies

(29) |

The canonical quantization of the GWs (27) can be written as , and the operator is Hermitian, i.e. . The annihilation and creation operators and satisfy the ordinary canonical commutation relations at the same time,

(30) |

The correlator for the GWs is defined as

(31) |

where is the dimensionless power spectrum for the GWs of each polarization.

Similarly, the power spectrum for the comoving curvature perturbation is given by

(32) |

A stochastic background of the GWs can be characterised by its energy density fraction Maggiore (2000); Allen and Romano (1999), which is defined as energy density of the GWs per unit logarithmic frequency

(33) |

where is the critical energy density at the conformal time , and is the Hubble parameter. The effective energy density of the GWs is usually defined as Boyle and Steinhardt (2008); Maggiore (2000)

(34) |

where the bracket means the time average over several periods of the GWs. When the relevant modes of the GWs are well inside the Hubble radius, one relates the and asBoyle and Steinhardt (2008); Maggiore (2000)

(35) |

where the overbar denotes the time average over several periods of the GWs.

## Appendix B Power Spectrum for the Induced Gravitational Waves

### b.1 Radiation-dominated era

There are many early works on the induced GWs during the radiation-dominated era, and the details of calculations are presented in, e.g. Bartolo et al. (2018b); Baumann et al. (2007); Ananda et al. (2007) and the references therein. Here we list the major formulae.

The relation between the 1st-order Bardeen potential and comoving curvature perturbations in the radiation-dominated era is,

(36) |

where the transfer function can be solved from the EoM for the Bardeen potential

(37) |

where .

We rewrite the source term (26) as

(38) |

where

(39) |

Note that, the prime still denotes the derivative w.r.t. the conformal time .

During the radiation-dominated era, it’s more convenient to use the canonical form for , and the correlator for is expressed as

(40) |

where the Green’s function satisfies the equation,

(41) |

and one can find the solution during the radiation-dominated era,

(42) |

From the source term (38), the two-point correlator can be expressed by the four-point correlator . By the Wick’s theorem, we obtain Baumann et al. (2007); Ananda et al. (2007)

(43) |

Note that since with , this connection does not contribute to the four-point correlator .

After having calculated the two-point correlator of the induced GWs (B.1), we get the power spectrum for the induced GWs

(44) |

where

(45) |

and the functions and are given by

(46) | ||||

(47) | ||||

respectively, where .

### b.2 Inflationary era

In the spatially flat gauge, one can write the source term in (21) contributed by the comoving curvature perturbation as

(48) |

where is the perturbed inflaton, and is the sound speed of . The relates to by . In the SSR mechanism, the modes of which are in the neighbourhood of the characteristic scale can be exponentially enhanced by the narrow resonance effect, and lead to the narrow peaks in the power spectrum (II). The source term in (24) during inflation reads off