Scalar Induced Gravitational Waves Review

Abstract

We provide a review on the state-of-the-art of gravitational waves induced by primordial fluctuations, so-called induced gravitational waves. We present the intuitive physics behind induced gravitational waves and we revisit and unify the general analytical formulation. We then present general formulas in a compact form, ready to be applied. This review places emphasis on the open possibility that the primordial universe experienced a different expansion history than the often assumed radiation dominated cosmology. We hope that anyone interested in the topic will become aware of current advances in the cosmology of induced gravitational waves, as well as becoming familiar with the calculations behind.

1. Introduction

Cosmology with Gravitational Waves (GWs) is a new and unique doorway to the primordial universe. The furthest in time we can observe with electromagnetic waves is around neutrino decoupling, roughly one second after the Big Bang and some time before Big Bang Nucleosynthesis (BBN). Prior to that, we have scarce evidence of the evolution of the primordial universe. See Ref. [1] for a review on possible expansion histories of the early universe. Nevertheless, we learned a lot from studies of the Cosmic Microwave Background1 (CMB). For instance, we know that there are Gaussian primordial fluctuations with an almost scale invariant spectrum [2,3], which by causality must have been generated before the Big Bang. This fact provides strong evidence that there was a period of accelerated expansion in the primordial universe, so-called inflation [4,5,6,7]. During inflation, quantum vacuum fluctuations are stretched out to the largest scales and become primordial fluctuations [8,9,10,11,12], which we later see as CMB anisotropies and galaxies. Thus, through the observation of primordial fluctuations we have access to inflation. However, the CMB only probes the largest scales and, therefore, only a small fraction of inflation. Any information on much smaller scales,2 and the last stages of inflation, is basically erased by complicated astrophysical processes. In contrast, GWs barely interact with intervening matter. This means that with GWs we can explore the universe before neutrino decoupling and even before the Big Bang towards the last stages of inflation. GWs provide a unique opportunity to complete our picture of the very early universe. A quantitative picture of the potential of cosmology with GWs can be found in Figure 1.

GWs generated by a cosmological process in the early universe will today appear as randomly distributed in all directions and with a very large number of unresolved sources.3 The superposition of all incoming unresolved GWs leads to the Stochastic GW Background (SGWB). Encouragingly, the typical frequency4 of such cosmic GWs falls right into the frequency range of future GWs detectors. For example, Pulsar Timing Arrays (PTA) [13,14,15,16,17,18], cover around ${10}^{-9}$∼${10}^{-7}$ $\mathrm{Hz}$, LIGO, VIRGO, KAGRA and ET [19], are sensitive roughly for 10∼10³ Hz. In between, we will have LISA [20,21], DECIGO [22,23,24], AION/MAGIS [25], Taiji [26] and Tainqin [27]. For an illustration see Figure 1.

Typical sources of cosmological GWs in the very early universe include (for more details see [29]): phase transitions, which may lead to collisions of bubbles or a universe filled with cosmic strings, resonances during reheating, quantum (gravity) fluctuations during inflation (so-called primordial GWs) and GWs induced from large primordial fluctuations. Such large primordial fluctuations may also collapse to form primordial black holes (PBHs). Out of these cosmological sources, GWs induced by primordial fluctuations, and to some extend PBHs, are our direct window to the latest stages of inflation. The information we could gain about inflation by the so-called induced GWs potentially covers scales around $k\sim {10}^7$–${10}^{18}{\mathrm{Mpc}}^{-1}$, otherwise inaccessible by any other probe. Even the absence of induced GWs will place new constraints on the primordial spectrum in unexplored regimes, potentially down to ${\mathcal{P}}_{\mathcal{R}}\sim {10}^{-4}$–${10}^{-5}$ [30]. On top of these very promising prospects, induced GWs are generated when primordial fluctuations re-enter the horizon sometime between the end of inflation and BBN. Since we have no evidence of the content and expansion history of the universe around that time, induced GWs not only provide access to the last stages of inflation but they also contain information on the content of the primordial universe. In this review, we will be opening to the possibility that the universe much before BBN was not dominated by radiation. In the future, information on the primordial spectrum obtained by induced GWs will complement those from other probes such as spectral distortions [31,32] in the multimessenger cosmology era [33].

1.1. Induced GWs History

The possibility of having GWs induced by density fluctuations was first noticed, as far as the author is aware, by K. Tomita in 1967 [34]. They were later rediscovered in the 1990s by Matarrese, Pantano and Saez [35,36] when studying second order cosmological perturbations in a dust dominated universe. Quite interestingly, in 1997 Matarrese, Mollerach and Bruni [37] noticed that induced GWs in a dust dominated universe suffer from gauge ambiguities. They proceeded to argue that only the oscillating part of the induced tensor modes were identified as true gravitational waves. We will discuss more about the gauge ambiguities in Section 7. However, Refs. [35,36,37] concluded that induced GWs were too small to be practically observed.

It was not until 20 years later that Ananda, Clarkson, and Wands [38] started to uncover the potential of induced GWs. In Ref. [38] they proposed to use induced GWs, generated in a radiation dominated universe, to constrain the spectral tilt of primordial fluctuations. A blue tilted primordial spectrum even with the CMB normalization might end up yielding large enough induced GWs. Some months later, Baumann, Steinhardt, Takahashi and Ichiki [39] looked into more detail at the transfer function of induced GWs taking into account the radiation-matter equality at around $z\sim 3400$ and the anisotropic stress due to neutrinos (see also [40] on the latter). They found that assuming the CMB normalisation of an almost scale invariant spectrum, modes with $k\sim k_{\mathrm{eq}}\sim 0.01{\mathrm{Mpc}}^{-1}$ were enhanced due to the matter dominated stage. Notably, the peak in the induced GW background can be larger than the primordial one if the tensor to scalar ratio is $r < 0.1$. However, the frequency of the peak of such GWs is extremely low, around $f_{\mathrm{eq}}\sim {10}^{-17}\mathrm{Hz}$. This is only observable perhaps by CMB B-mode polarization experiments or cosmic shear [41]. Nevertheless, this was an indication that an early epoch of dust domination might yield interesting results. In the same period, Martineau and Brandenberger [42] derived a lower bound to induced GWs for various inflationary models and alternatives, although they referred to it as lower bound from “backreaction” of scalar fluctuations. Another application of GWs induced by primordial fluctuations in the curvaton scenario, sourced between the end of inflation and the curvaton decay, was studied by Bartolo, Matarrese, Riotto and Väihkönnen [43]. An action formalism for a scalar field acting as a perfect fluid and its induced GWs was studied by Boubekeur, Creminelli, Noreña and Vernizzi [44].

An important realization was done by Saito and Yokoyama in 2008 [45,46]: if the primordial spectrum of fluctuations is large enough on small scales, it does not only induce GWs but might also collapse to form PBHs. This means that the induced GW spectrum can be used in the future to place bounds on the PBH abundance. Even more, if PBHs were found, one should also find the induced GW counterpart. These ideas were further pursued by Bugaev and Klimai [47,48,49]. Extending the work of Saito and Yokoyama a bit further, Assadullahi and Wands proposed induced GWs as probes of the primordial spectrum [50]. They also considered the early dust dominated case [51] which confirmed the large enhancement of induced GWs. The same group, now including Arroja, Assadullahi, Koyama and Wands [52], studied the matching conditions on superhorizon scales for such GWs induced at second order. Later, in 2012 the induced GWs from particular models of inflation were studied by Alabidi, Kohri, Sasaki and Sendouda in radiation [53] and dust domination [54]. Around the same time, a curvaton scenario with a blue tilted primordial spectrum yielding large induced GWs was proposed by Kawasaki, Kitajima and Yokoyama [55]. More investigations on the upper bounds of induced GWs by overproduction of PBHs were done by Nakama and Suyama [56,57]. Other works related to induced GWs include the decay of two curvatons by Suyama and Yokoyama [58], the impact of anisotropic stress due to free streaming particles by Saga, Ichiki and Sugiyama [59] and the B-modes due to induced GWs in the CMB by Fidler et al. [60].

In 2016, LIGO reported the detection of GWs from the merger of a binary black hole [61] and the amount of new works on induced GWs and their PBH counterpart exploded. The amount of works is so large that we will not attempt to go through all of them in chronological order but instead we will classify them into seven main directions below. We give a more detailed discussion in the corresponding sections.

–

General semi-analytical formulation: Since induced GWs are a second order effect one needs to integrate in time and momenta over the linear evolution of primordial fluctuations. The analytical transfer functions for radiation domination are derived in Refs. [62,63] and later generalized to constant equation of state parameter in Ref. [64];

–

Induced GWs for different expansion histories and different contents of the universe: Induced GWs may have been generated in a non-radiation dominated universe. This leaves characteristic signatures in the induced GW spectrum. Studies in early matter era can be found in Refs. [65,66,67]. The extension to an early PBH dominated epoch is investigated in Refs. [68,69,70,71]. More general thermal histories are studied in Refs. [64,72,73,74,75,76]. The impact of additional free streaming particles is studied in Ref. [77];

–

Induced GW spectral features: There are cases where the induced GW spectrum may be investigated semi-analytically. These are for example, the low frequency tail [74,78,79], the UV tail [80,81] and the log-normal peak in the primordial spectrum [82]. Furthermore, the primordial spectrum may also present oscillatory features which are captured into the induced GW spectrum [83,84,85,86,87]. On top of that, large primordial non-Gaussianities may have a non-trivial impact on the induced GW spectrum [81,88,89,90,91,92,93]. Other effects include: anisotropic non-gaussianities, which may be a source of superhorizon tensor modes [94], resonances that may occur during inflation, which enhance the induced GW spectrum [95,96,97,98], and non-Bunch Davies initial conditions in inflation [99], although the latter does not yield an observable signature;

–

Explanations of current observations: Induced GWs have been extensively used as counterpart of the PBH scenario as a totality or a fraction of dark matter. For example, the induced GWs from various inflationary models can be found in Refs. [100,101,102,103,104,105,106,107,108,109,110,111,112,113,114,115,116,117,118] and [62,116,119,120] in the context of Higgs inflation. In particular, the large induced GW counterpart to PBH as the totality of dark matter is studied in Refs. [88,110,121,122,123,124,125]. Other possibilities include an explanation to the LIGO observations [121,126,127,128] and the NANOGrav results [129,130,131,132,133,134,135];

–

Current and future GW constraints: It is important to place constraints on the current absence of induced GWs and also asses future capabilities to constraint/find different models. A study using current PTA and LIGO data on SGWBs can be found in Refs. [136,137,138,139] and an analysis of future GW prospects in Refs. [33,75,140,141,142,143,144,145,146];

–

SGWB anisotropies: In the far future we may not only observe the GW spectrum but also its anisotropies. The anisotropies due to induced GWs associated to PBH formation in the case of a monochromatic mass function are investigated in Refs. [122,123,147,148];

–

The gauge issue of induced GWs: Tensor modes are subject to gauge ambiguities at second order due to mode mixing. Since the work of Hwang, Jeong and Noh [149] in 2017 there has been an extensive discussion on the gauge issue of induced GWs [150,151,152,153,154,155,156,157,158,159,160,161,162]. The source of the problem and the applicability of predictions is by now well understood, although the gauge issue persists in the strictest sense.

1.2. Structure and Scope of the Review

Induced GWs are the most promising probe of the primordial universe. In this review, we will revisit and unify the current analytical formulation and main predictions for induced GWs. Note that the scope of this review is basically limited to the author’s works on induced GWs [64,70,71,74,81,132,150,161], although we discuss and properly cite the relevant literature. As interesting as it is, we will not enter in the details of PBH formation nor on the GW data analysis. Nevertheless, we briefly discuss these relevant topics and give appropriate credit, when needed.

For convenience, we proceed to briefly explain the organization of the review. In this way, the reader interested in certain aspects may jump directly to the relevant section. In Section 2 we present an intuitive picture of the induced GW generation. We estimate their typical frequency, the spectral shape and we briefly comment on their PBH counterpart for different expansion histories in the primordial universe. The details of the calculations and derivation of analytical approximations are reviewed in Section 3 and Section 4. The analytical approximations for the induced GW spectrum for typical primordial curvature spectra are presented in Section 5. The special case of the dust dominated universe is discussed in some detail in Section 6. We also include the case where PBH dominate the universe in Section 6.2. The gauge issue of induced GWs is explained in Section 7. We list other GW counterparts associated to PBHs in Section 8 and discuss the future observational prospects in Section 9. We also provide a summary of the main formulas in Section 10. These formulas are ready for calculating the induced GW spectrum for general primordial spectra. Thus, the reader interested solely in the final formulas of the induced GW may jump directly to Section 5 and/or Section 10. We end with a discussion on future directions in the conclusions Section 11. We present many details of the numerical factors, formulas and calculations used in the main text in the Appendices.

We work in natural units where $\hslash =c=1$. We keep the reduced Planck mass $M_{\mathrm{pl}}={\left(8\pi G\right)}^{-1/2}$ in the equations, except in Section 2.2 and Section 4 where for simplicity in the calculations we set it to $M_{\mathrm{pl}}=1$. At the beginning of each section we list the literature in which the section is based, except for Section 8 and Section 9 which are relatively short. We provide a list of useful related reviews below.

A list of other useful reviews. The context of induced GWs is very broad. It covers cosmological perturbation theory up to second order, stochastic GW backgrounds, primordial black holes and primordial non-Gaussianity. Unfortunately, all the details cannot be covered in this review. Instead, we list below existing reviews on topics directly related to induced GWs, which have been useful in writing this manuscript.

• Cosmological perturbation theory. The classical reviews on cosmological perturbation theory at linear order are the one by Kodama and Sasaki [12] and by Mukhanov, Feldman and Brandenberger [163]. A typical reference for second order cosmological perturbation theory is the review by Malik and Wands [164]. Since then there have been many other reviews on cosmological perturbation theory. Some that we found particularly up to date and useful are Refs. [165,166,167,168]. For a take on multi-field inflation we suggest Ref. [169]. For primordial features in the primordial spectrum we refer the reader to Ref. [170];
• Stochastic GW backgrounds. Induced GWs are not the only source of GW backgrounds. Reviews on the different cosmic and astrophysical sources can be found in the review by Caprini and Figueroa [29] and by Christensen [171]. A useful collection of cosmic GW spectra can be found in Kuroyanagi, Chiba and Takahashi [172], although it is technically not a review. A review focused on GWs from inflation is given by Guzzetti, Bartolo, Liguori and Matarrese [173];
• Primordial black holes. The literature on primordial black holes is very vast and currently under refinement. A review that has been used in particular is the one by Sasaki, Suyama, Tanaka and Yokoyama [174]. Other interesting reviews are Refs. [175,176,177,178]. A complementary review on induced GWs with more focus on the PBH counterpart is given by Yuan and Huang [179];
• Primordial non-Gaussianity. Although quantum fluctuations during inflation are drawn from a Gaussian distribution, they can develop small departures from such Gaussian distribution due to gravitational or general interactions. The reader interested in primordial non-Gaussianities may check the reviews in the context of inflation and CMB observations, e.g., [180,181,182,183];
• Alternative expansion histories. A recent review encompassing many of the new physics of a primordial universe which is not filled with radiation is given in Ref. [1].

2. Estimates and Intuitive Picture

In this section we present the intuitive physical picture behind induced GWs. We do not intend to do a rigorous derivation but to lay out the basic assumptions behind the calculations. Before we go into the details of the induced GWs themselves, we first review in Section 2.1 what is considered to be the “observable” quantity for cosmic GWs. This will also be useful when dealing with the gauge issue in Section 7. We later in Section 2.2 derive intuitively, based on heuristic calculations and order of magnitude estimations, the main features of the induced GW spectrum in general cosmological backgrounds. Lastly in Section 2.3, we give the relation between the typical frequencies of the induced GWs with the corresponding PBH masses also in general cosmological backgrounds.

Main references: Section 2.1 is mainly based on Misner, Thorne and Wheeler’s [184] and Maggiore’s [185] books, while the original work is essentially due to Isaacson [186,187]. We follow the notation of, e.g., Ref. [78,161]. Section 2.2 essentially follows from [74] although it builds up from previous works in the literature, e.g., [39,63,64,78]. Section 2.3 is mainly extracted from [64,132,174]. All relevant equations can be found in Appendix E.

2.1. The Spectral Density of GWs in Cosmology

A GW detector, such as an interferometer, measures the GW strain at a given frequency and time. In the case of the SGWB, one is measuring the time and/or ensemble average of the superposition of all incoming GWs. For convenience, cosmologists often use the energy density fraction of GWs as the “observable”. However, it is not an easy task to theoretically associate an energy density to GWs in general situations. This is because in general relativity there is no local notion of energy of the gravitational field. By the equivalence principle, one may always go to a locally flat frame and erase any sign of gravity. Of course, we know that GWs carry momentum and energy (after all, they have been detected). Thus, indeed one can define the energy of GWs in some limits of interest. For instance, when the frequency of the GW is much higher than that of a slowly varying background [185].

Relevant to the later discussion, let us briefly review the case of GWs propagating in a Ricci flat spacetime, i.e., a spacetime with no matter sources. We can then split the total metric as $g_{\mu \nu }={\overline{g}}_{\mu \nu }+h_{\mu \nu }$ where ${\overline{g}}_{\mu \nu }$ is the background metric and $h_{\mu \nu }\ll 1$ are the GWs, treated as a high frequency perturbation. Then, by the non-linear nature of gravity, these GW backreact onto the background metric. Mathematically speaking, if we expand the Einstein tensor up to second order we will find terms quadratic in $h_{\mu \nu }$ which can be thought of as a “matter” source. This means that, if we somehow “integrate out” the high frequency perturbation, we can write5

$${\overline{G}}_{\mu \nu }\left[{\overline{g}}_{\mu \nu }\right]=M_{\mathrm{pl}}^{-2}{t}_{\mu \nu }^{\mathrm{GW}}\left[h_{\mu \nu }\right],$$

(1)

where $G_{\mu \nu }$ is the Einstein tensor and $t_{\mu \nu }^{\mathrm{GW}}$ is the (pseudo)-energy momentum tensor of GWs. A way to “integrate out” the high frequency modes is by taking an average. For instance, we may focus on a small box and take the volume (and/or time) average. The box shall not be too small, so that several GW wavelengths fit in. In this case, we have that the (pseudo)-energy momentum tensor of GWs reads [186,187]

$$t_{\mu \nu }^{\mathrm{GW}}=\frac{M_{pl}^2}{4}〈{\partial }_{\mu }{h}^{\alpha \beta }{\partial }_{\nu }{h}_{\alpha \beta }-\frac{1}{2}{\overline{g}}_{\mu \nu }{\partial }_{\sigma }{h}^{\alpha \beta }{\partial }^{\sigma }{h}_{\alpha \beta }〉,$$

(2)

where all contractions of indices are done with the background metric.

In cosmology we face the extra difficulty of having a matter source to the Einstein equations which drives the expansion of the universe. Because of the expansion, we also have a cosmological horizon roughly given by $1/H$ where H is the Hubble parameter and quantifies the expansion rate. For these reasons, the definition of GWs in a cosmological background is more subtle than in Ricci flat space times. Nevertheless, with some tweaks we shall write in analogy

$$G_{\mu \nu }=M_{\mathrm{pl}}^{-2}\left(T_{\mu \nu }+t_{\mu \nu }^{\mathrm{GW}}\right),$$

(3)

where $T_{\mu \nu }$ is the matter energy-momentum tensor. This time, however, we must be more careful on the averaging procedure. Furthermore, we must introduce a slightly different, but important, notation. First, in the helicity decomposition of the metric, the transverse-traceless component of the spatial metric is referred to as tensor mode. A tensor mode with a co-moving wavenumber k stays constant if its physical size is larger than the cosmological horizon, i.e., $k/a\ll H$. In the opposite limit, a tensor mode behaves as a GW if its physical size is much smaller than the cosmological horizon, that is $k/a\gg H$. Thus, Equation (3) only makes sense for the high frequency tensor modes deep inside the horizon. This also means that the small boxes we used to take the average have to be much smaller than the horizon size but still larger than several GW wavelengths. In addition to that, in cosmology we often deal with a SGWB and, by ergodicity, we may take an ensemble average over all small boxes. Under these assumptions we may use Equation (2) to compute the energy density of GWs in cosmology. Be aware though that these assumptions are strictly speaking not enough as we shall argue in Section 7.

Proceeding with the cosmological SGWB, we consider GWs propagating in a Friedmann-Lemaître-Robertson-Walker (FLRW) metric, namely6

$$d{s}^2=-d{t}^2+a^2\left(t\right)({\delta }_{ij}+h_{ij})d{x}^{i}d{x}^j,$$

(4)

where $a\left(t\right)$ is the scale factor. The evolution of the scale factor is dictated by the Friedmann equations, which we provide in Appendix E. Using Equation (2) under the ensemble average procedure, we find that the energy density of GWs (the 00 component of $t_{\mu \nu }^{\mathrm{GW}}$) in Fourier space (see Appendix D for notations) reads

$${\rho }_{\mathrm{GW}}=M_{\mathrm{pl}}^2\int dlnk\frac{k^3}{16{\pi }^2}\left\{\sum _{\lambda }{\langle {\dot{h}}_{\mathit{k},\lambda }{\dot{h}}_{-\mathit{k},\lambda }\rangle }^{!}+\frac{k^2}{a^2}{\langle h_{\mathit{k},\lambda }{h}_{-\mathit{k},\lambda }\rangle }^{!}\right\},$$

(5)

where $\dot{}\equiv d/dt$, $\lambda$ is the GW polarization and we already used that the coincident ensemble average of two tensor modes in a homogeneous and isotropic universe is proportional to a Dirac delta,7 namely

$$\langle h_{\mathit{k},\lambda }{h}_{{\mathit{k}}^{\prime },\lambda }\rangle ={\langle h_{\mathit{k},\lambda }{h}_{-\mathit{k},\lambda }\rangle }^{!}{\left(2\pi \right)}^3{\delta }^{\left(3\right)}(\mathit{k}+{\mathit{k}}^{\prime }).$$

(6)

The notation “!” in Equations (5) and (6) is to indicate that the Dirac delta has been factored out, so that ${\langle h_{\mathit{k},\lambda }{h}_{-\mathit{k},\lambda }\rangle }^{!}$ is directly related to the power spectrum, which by isotropy is a function of $k=\left|\mathit{k}\right|$ only. Equation (5) is the total energy density of the GWs filling the universe. It is often more convenient to use the spectral density fraction, i.e., the GW power per logarithm of a given wavenumber k (or frequency) over the critical density, which is defined as

$${\Omega }_{\mathrm{GW}}\left(k\right)\equiv \frac{1}{3M_{\mathrm{pl}}^2{H}^2}\frac{d{\rho }_{\mathrm{GW}}}{dlnk}=\frac{k^2}{12a^2{H}^2}\sum _{\lambda }{\mathcal{P}}_{h,\lambda }\left(k\right),$$

(7)

where $H=\dot{a}/a$. This is often what is used to compare theoretical predictions with current constraints and future observational prospects. Note that in Equation (7) we used the definition of the dimensionless power spectrum, concretely

$${\langle h_{\mathit{k},\lambda }{h}_{-\mathit{k},\lambda }\rangle }^{!}=\frac{2{\pi }^2}{k^3}{\mathcal{P}}_{h,\lambda }\left(k\right).$$

(8)

We also used the approximation that for a freely propagating wave we have that

$$\begin{array}{c}\hfill {\dot{h}}_{\mathit{k},\lambda }\approx \frac{k}{a}{h}_{\mathit{k},\lambda }.\end{array}$$

(9)

With these ingredients we can use Equation (7) to estimate the amplitude and spectral shape of induced GWs.

2.2. GWs Induced by Primordial Fluctuations

Let us derive some estimates by considering a rough, yet instructive, toy model for induced GWs. Take a tensor mode $h_{ij}$ in a FLRW background (4) sourced by scalar field fluctuations $\delta \phi$. Its equations of motion, in terms of conformal time $d\tau =dt/a$, read

$$h_{ij}^{″}+2\mathcal{H}{h}_{ij}^{\prime }-{\partial }_k{\partial }^k{h}_{ij}={\mathcal{P}}^{ab}{}_{ij}\left\{T_{ab}\right\}\approx {\mathcal{P}}^{ab}{}_{ij}\left\{{\partial }_a\delta \phi {\partial }_b\delta \phi \right\},$$

(10)

where ${}^{\prime }\equiv d/d\tau$, $\mathcal{H}=H/a=a^{\prime }/a$, $T_{ij}$ is the spatial components of the energy-momentum tensor of a scalar field (see Equation (45)) and ${\mathcal{P}}^{ab}{}_{ij}$ is the transverse-traceless projector defined in Appendix E. Without understanding much about the evolution of these fluctuations, we can already see that

$${\mathcal{P}}_h\equiv \sum _{\lambda }{\mathcal{P}}_{h,\lambda }\propto {\mathcal{P}}_{\delta \phi }^2.$$

(11)

If the scalar field fluctuations were sourced during inflation, we can relate them to the curvature perturbation $\mathcal{R}$ so that ${\mathcal{P}}_{\delta \phi }\sim {\mathcal{P}}_{\mathcal{R}}$. Furthermore, the spectral density of GWs in a radiation dominated universe remains constant, as GWs behave as radiation themselves. Thus, we can boldly estimate that the amplitude of induced GW spectral density measured today is given by

$${\Omega }_{\mathrm{GW}}^{\mathrm{induced}}{h}^2\sim \frac{1}{12}{\Omega }_{r,0}{h}^2\times {\mathcal{P}}_{\mathcal{R}}^2,$$

(12)

where ${\Omega }_{r,0}\equiv {\rho }_{r,0}/\left(3H_0^2{M}_{\mathrm{pl}}^2\right)$ is the density fraction of radiation, $H_0$ is the Hubble rate today and $h\equiv H_0$/(100 km/Mpc/s). We have introduced ${\Omega }_{r,0}$ to take into account the dilution of the GWs (as radiation) as the universe expands and goes through the cold dark matter and dark energy dominated stages. We give more details on this factor in Section 5. h considers the uncertainty in the exact value of $H_0$. From the latest Planck 2018 results [2] we have that ${\Omega }_{r,0}{h}^2\approx 4\times {10}^{-5}$. With all the above simplifications, we arrive at

$${\Omega }_{\mathrm{GW}}^{\mathrm{induced}}{h}^2\sim {10}^{-6}{\mathcal{P}}_{\mathcal{R}}^2.$$

(13)

The estimate (13) can be now contrasted with future observational prospects. The sensitivity of future GW detectors after collecting data over several years might reach ${\Omega }_{\mathrm{GW}}^{\mathrm{induced}}{h}^2\sim {10}^{-16}$, if we use the power-law sensitivity curves of Thrane and Romano [28]. This optimistic value is for the best DECIGO sensitivity [23,188]. Then, in the absence of an induced GW detection we can use Equation (13) to place an optimistic upper bound to the primordial spectrum on the corresponding scales (roughly at ${10}^{16}$–${10}^{12}$ ${\mathrm{Mpc}}^{-1}$ ) of about ${\mathcal{P}}_{\mathcal{R}} < {10}^{-5}$. This is already three orders of magnitude better than current constraints from the absence of PBHs (around ${\mathcal{P}}_{\mathcal{R}} < {10}^{-2}$ [174]), although it is still four orders of magnitude above the CMB extrapolation which would be around ${\mathcal{P}}_{\mathcal{R}}\sim {10}^{-9}$ [3]. The most interesting point though is that PBHs which formed by the collapse of large primordial fluctuations, which requires ${\mathcal{P}}_{\mathcal{R}}\gtrsim {10}^{-2}$, must have an observable large induced GW counterpart. Thus, the absence of induced GWs completely rules out the PBH scenario from large primordial fluctuations [189,190]. Note that PBH could have formed by other mechanisms though. For example, PBH could form in first order phase transitions [191,192], tunneling during inflation [193], the collapse of Q-balls [194,195] and might also be the result of long range interactions stronger than gravity [196,197,198].

In the above estimations we have completely neglected the evolution of the scalar fluctuations since the end of inflation onwards, which is by no means irrelevant. As shown in Refs. [51,54,65,66], a period of dust domination might substantially enhance the amplitude of induced GWs. The amplification is such that GWs induced by tiny primordial fluctuations might observable [54,66]. This case exemplifies the importance of the expansion history of the primordial universe. At this point, we would like to emphasise that contrary to studies of the CMB we have scarce evidence of the content of the primordial universe. Therefore, when deriving theoretical predictions, and eventually contrasting them with the data, it would be a good idea to be open and acknowledge our ignorance by allowing a free equation of state parameter $w\equiv \rho /p$ and a free propagation speed of fluctuations $c_s$, which describe the unknown content of the primordial universe. A nice review on the implications of various expansion histories can be found in Ref. [1]. For these reasons, let us estimate the general impact of the equation of state w and sound speed $c_s$ onto the induced GW spectrum. To do that, we continue with the toy model of the scalar field. This time, however, we will consider the evolution of $\delta \phi$. Fluctuations of a K-essence massless scalar field roughly follow the Klein-Gordon equation

$$\delta {\phi }_k^{″}+2\mathcal{H}\delta {\phi }_k^{\prime }+c_s^2{k}^2\delta {\phi }_k=0,$$

(14)

where for simplicity we considered a constant w and $c_s$. The former implies constant expansion rate, i.e.,

$$ϵ\equiv -\frac{\dot{H}}{H^2}=\frac{3}{2}(1+w)=\mathrm{constant}.$$

(15)

Such a constant expansion rate implies a scale factor which goes as a power-law of conformal time, explicitly from the equations in Appendix E we have

$$a\propto {\tau }^{1+b}\mathrm{where}b=\frac{1-3w}{1+3w}.$$

(16)

We introduced the parameter b for later convenience. $b=0$ corresponds to a universe filled with radiation ($w=1/3$). Then, $b < 0$ and $b > 0$ respectively correspond to a stiffer and softer fluid. Note that if we allow $w\gg 1$ the range of b goes as $-1 < b < \infty$ for $\infty > w > -1/3$.

For the sake of simplicity, we solve the Klein-Gordon Equation (14) in the super-(sound)horizon, $c_{s}k\ll \mathcal{H}$, and sub-(sound)horizon, $c_{s}k\gg \mathcal{H}$, limits. This leads us to

$$\delta {\phi }_k\approx \delta {\phi }_i\left\{\begin{array}{cc}\mathrm{constant}\hfill & \hfill c_{s}k\tau \ll 1\\ a^{-1}{e}^{i{c}_{s}k\tau }\hfill & \hfill c_{s}k\tau \gg 1\end{array}\right.,$$

(17)

where $\delta {\phi }_i$ is an initial value. The above solutions are all we need to solve for the induced GWs. For analytical viability, let us take that there are fluctuations of the scalar field only in a single scale $k_p$. Namely the scalar fluctuations are Dirac-delta like in Fourier space, i.e.,

$$\delta {\phi }_i\sim \delta {\phi }_p\times \delta (ln(k/k_p)),$$

(18)

where $\delta {\phi }_p$ is the primordial amplitude, e.g., set by inflation. All the above simplifications lead to a simple differential equation for the tensor modes given by

$$h_k^{″}+2\mathcal{H}{h}_k^{\prime }+k^2{h}_k\approx k_p^2\delta {\phi }_p^2.$$

(19)

We readily see that, if we focus on $\delta {\phi }_p=\mathrm{constant}$ and we evaluate $h_k$ at horizon crossing, that is when the mode k is of the size of the comoving horizon, i.e., $k=aH$, then the amplitude of induced GWs is proportional to $h_k\sim \delta {\phi }_p^2$ for $k\sim k_p$, as used in our previous estimate (12). Now, to derive the spectral features we can focus on two interesting limits: $\left(i\right)$ the infra-red (IR) tail for modes with $k\ll k_p$ and $\left(ii\right)$ the near peak behaviour for modes with $k\sim k_p$.

2.2.1. The IR Tail

If we focus on modes with $k\ll k_p$ there will be no competition between the $k^2$ and $k_p^2$ terms nor any resonance. Thus, we can peacefully solve the tensor modes in the superhorizon ($k\tau \ll 1$) and subhorizon ($k\tau \gg 1$) regimes. First, in the superhorizon regime we find that

$$h(k\tau \ll 1)\approx \left\{\begin{array}{cc}{\left(k_p\tau \right)}^2\hfill & \hfill c_s{k}_p\tau \ll 1\\ \mathrm{constant}+{\left(k_p\tau \right)}^{-2b}\hfill & \hfill c_s{k}_p\tau \gg 1\end{array}\right..$$

(20)

What Equation (20) tells us is that tensor modes initially grow as ${\tau }^2$ due to a constant source. Then, after the scalar source enters the sub(sound)-horizon regime, i.e., at $c_s{k}_p{\tau }_p=1$, the tensors developed a constant amplitude. Now, if $b > 0$ ($w<1/3$) the scalar source decays much faster than the cosmological background, that is $\delta {\phi }_k^2/{\mathcal{H}}^2\propto {\left(k\tau \right)}^{-b}\to 0$, and the tensor growth stops. However, for $b<0$ ($w > 1/3$) the scalar source decays slower than the background expansion, which yields a second superhorizon growth of the tensor modes [74]. To follow tensor mode evolution, the second line of Equation (20) has to be matched at horizon crossing for the tensor mode ($k\tau =1$) with the subhorizon solution. By doing so, we arrive at

$$h(k\tau \gg 1)\approx {\left(k\tau \right)}^{-1-b}{e}^{ik\tau }\left(\mathrm{constant}+{\left(k_p/k\right)}^{-2b}\right).$$

(21)

The final ingredient is that somehow the universe transitions to a radiation dominated universe to recover the standard cosmology. For simplicity, we take an instantaneous reheating at $\tau ={\tau }_{\mathrm{rh}}$ with the corresponding reheating scale $k_{\mathrm{rh}}={\mathcal{H}}_{\mathrm{rh}}$. From that moment on the spectral density of induced GWs stays constant and we can derive our predictions. Inserting Equation (21) evaluated at $\tau ={\tau }_{\mathrm{rh}}$ into Equation (5) we find that the IR slope of the induced GW spectrum is in general given by

$$\frac{dln{\Omega }_{\mathrm{GW}}^{\mathrm{induced}}}{dlnk}(k_{\mathrm{rh}}\ll k\ll k_p)\approx 3-2\left|b\right|.$$

(22)

We later find that for $\beta \in \mathbb{Z}$ there is a logarithmic correction [74]. Additionally, if the scalar spectrum of fluctuations is really a Dirac delta, the spectrum slope is given by $2-2\left|b\right|$ due to the unphysical divergence of $\delta {\phi }_p$ at $k=k_p$. Note that for tensor modes which entered during the radiation domination era we have to match the first line of (20) to the subhorizon solution. This yields a spectral index given by

$$\frac{dln{\Omega }_{\mathrm{GW}}^{\mathrm{induced}}}{dlnk}(k\ll k_{\mathrm{rh}})\approx 3,$$

(23)

which is in agreement with the general results of Ref. [78]. These naive derivations will be later recovered from the analytical formulas in Section 4. It is interesting to see that for $b > 3/2$ ($w<-1/15$) the induced GW spectrum peaks at $k\sim k_{\mathrm{rh}}$ which is clearly different than the peak in the primordial spectrum at $k\sim k_p$.

The intuitive physical picture behind (22) is the following. In a universe with constant b, modes of a massless field (e.g., the tensor modes $h_{ij}$ or the scalar field $\delta \phi$) which enter the horizon early get more diluted than modes which enter later. For example, consider that right after horizon crossing, i.e., when $k={\mathcal{H}}_k$, we have that the energy density of the massless field starts to dilute as $\rho \propto {(a_k/a)}^4$. Since we have that $a_k\propto {\tau }_k^{1+b}\propto k^{-1-b}$ then $\rho \propto a_k^4\propto k^{-4b}$ and $\rho /H^2\propto a^{-2b/(1+b)}\propto k^{-2b}$. Thus, if the IR tail of induced GWs in radiation domination goes as $k^3$, we expect ${\Omega }_{\mathrm{GWs}}\propto k^{3-2b}$ due to dilution [78]. However, as we have shown, this is incomplete for induced GWs from a peaked spectrum. For $b<0$, the density fraction of the scalar field ${\rho }_{\delta \phi }$ grows as $a^{-2b/(1+b)}$. This means that induced tensor modes which are generated later have a larger source and a larger amplitude than those generated earlier. In other words, for $b<0$ small k have a larger amplitude than large k by a factor $a_k^{-4}\propto k^{4b}$. Then, for $b<0$ the density fraction of induced GWs goes as ${\Omega }_{\mathrm{GWs}}\propto k^{3-2b+4b}\propto k^{3+2b}$ [74]. This gives the general formula ${\Omega }_{\mathrm{GWs}}\propto k^{3-2\left|b\right|}$ that we derived in (22). Whether this absolute value of b is unique to induced GWs is something yet to explore.

Before moving on to the near peak limit, let us comment on the implications of Equation (22). While the far most IR tail of the induced GW spectrum indeed must go as $k^3$, we may have an intermediate regime where ${\Omega }_{\mathrm{GWs}}\propto k^{3-2\left|b\right|}$. This has interesting consequences when comparing with GWs from other sources. Let us model a general GW spectrum by a broken power-law around a characteristic scale $k_c$, as done in Ref. [172]. Then, for example, the IR spectrum of induced GWs for $b=2$ ($w=-1/9$) resembles that of GWs generated by a strong first order phase transition which have ${\Omega }_{\mathrm{GWs}}\propto k^{2.8}$ for $k<k_c$ and ${\Omega }_{\mathrm{GWs}}\propto k^{-1}$ for $k > k_c$. This gives another very good reason to be open about a primordial universe with a different equation of state than radiation.

2.2.2. The near Peak Regime

Looking closely at Equation (19), we might suspect that there could be a resonance if the wavenumber k of a tensor mode approaches that of twice the scalar mode $2k_p$. This resonance may be understood in two ways. First, we know that a harmonic oscillator with an external force has a resonance when the period of the force matches that of the harmonic oscillator. If that occurs, the amplitude of the oscillations diverges in time. In short, the external force is kicking the oscillator at the right times. The second way to see it is that two scalar modes with physical momentum $c_s{k}_p$ have a preferred window to produce a tensor mode with physical momentum equal to twice the scalar momenta, i.e., $k=2c_s{k}_p$, if allowed by momentum conservation. To see this mathematically, let us redefine the tensor modes in Equation (19) to remove the Hubble friction by

$$h_k=v_k/a.$$

(24)

With this new variable $v_k$ Equation (19) now reads

$$v_k^{″}+\left(k^2-\frac{b(1+b)}{{\tau }^2}\right)v_k=\frac{k_p^2}{a}\delta {\phi }_p^2{e}^{2i{c}_s{k}_p\tau },$$

(25)

where we only focused on the oscillatory behaviour of the scalar field for $c_s{k}_p\tau \gg 1$. The particular solution to $v_k$ can be found by the Green’s function method (see Appendix B for the details). Neglecting the term proportional to $1/{\tau }^2$ in Equation (25), the homogeneous solutions to $v_k$ are $e^{\pm ik\tau }$. With such homogeneous solutions, the particular solution by the Green’s method reads

$$v_k\approx {\int }_{{\tau }_i}^{\tau }d\tilde{\tau }sin\left(k(\tau -\tilde{\tau })\right){\tilde{\tau }}^{-1-b}{e}^{2i{c}_s{k}_p\tilde{\tau }}.$$

(26)

We clearly see from the integrand in Equation (26) that there is a resonance when the frequency of the sine equals that of the exponential. Then the oscillations interfere and leave the integral of a power-law. Since we are only interested in how the amplitude diverges, let us only focus on the divergent part of (26) which is given by

$$v_k\propto \int d\tilde{\tau }{\tilde{\tau }}^{-1-b}{e}^{-i\tilde{\tau }(k-2c_s{k}_p)}\propto \Gamma [-b]{(k-2c_s{k}_p)}^b.$$

(27)

We see that when $b\le 0$ the resonance is effective and the amplitude diverges. Note that for $b=0$ the divergence is logarithmic. When $b > 0$ the resonance is not effective enough to cause a divergence as the amplitude of the integrand decays faster than ${\tau }^{-1}$. However, since we so far used a Taylor expansion around $k\sim 2c_s{k}_p$ we expect that for $1 > b > 0$ Equation (27) still yields the leading contribution after the constant term. Thus, for $1 > b > 0$ we also expect a peak at around $k\sim 2c_s{k}_p$, albeit not divergent. For $b > 1$ the situation is less clear and a detailed analysis is necessary. Nevertheless for $b<1$, we can roughly write that the tensor modes in the near peak regime have an amplitude given by

$$h_k(k\sim 2c_s{k}_p,b<1)\propto \mathrm{constant}-{(k-2c_s{k}_p)}^b.$$

(28)

This also implies that the induced GW spectrum has a peak near $k\sim 2c_s{k}_p$ with a growth rate proportional to

$$\begin{array}{c}\hfill {\Omega }_{\mathrm{GW}}^{\mathrm{induced}}(k\sim 2c_s{k}_p)\propto \left\{\begin{array}{cc}{\left|k-2c_s{k}_p\right|}^{-2\left|b\right|}\hfill & \hfill b<0\\ \mathrm{constant}-{\left|k-2c_s{k}_p\right|}^b\hfill & \hfill 1 > b > 0\end{array}\right..\end{array}$$

(29)

These estimates will be later checked by the general analytical solutions in Section 4. Thus, the resonant peak in the induced GW spectrum has information on the sound horizon $c_s{k}_p$ by the position of the peak and on the expansion history b by the growth rate around the peak. This type of resonant structure is characteristic of induced GWs. Thus, the observation of such resonant peak could give strong evidence for induced GWs. Then, the sound speed of fluctuations $c_s$ and the peak $k_p$ could be disentangled by the simultaneous observation of the PBH counterpart, which we proceed to discuss. Note that in (29) the radiation domination case $b=0$ is “special” because it is exactly when the amplitude of the source term is proportional to the background density.

2.3. Primordial Black Hole Counterpart

In this section we briefly introduce the PBH counterpart to the induced GWs. It is by no means a review on PBH formation and it does not intend to be so. A detailed review on PBHs can be found in Ref. [174]. Other interesting reviews are listed at the Section 1.2. PBHs from large primordial fluctuations will form in those Hubble patches where the density contrast $\delta \rho /\rho$ is above a certain threshold8 ${\delta }_c$ [190,203,204,205,206,207,208], which depends on the EoS parameter w at horizon re-entry [209,210,211,212,213,214,215,216]. Since fluctuations generated during inflation closely follow a Gaussian distribution centred around zero, only the tail of the distribution has a large enough amplitude to collapse to PBH. This means that PBH formation is an extremely rare event. The fraction of Hubble patches where PBH forms is exponentially suppressed and sensitive to the root mean squared of the distribution, which is related to the power spectrum. This is why we generally need ${\mathcal{P}}_{\mathcal{R}}\sim {10}^{-2}$ to form a substantial fraction of PBHs. Precise estimations of the fraction of PBH is an active and ongoing research field, e.g., see Refs. [202,214,215,216,217,218,219,220,221] and references therein. For our purposes, it is enough to know that only large induced GWs have a large associated fraction of PBH and, inversely, the absence of large induced GWs rules out PBHs from large primordial fluctuations.

For simplicity, we continue with the assumptions of the previous section and take a primordial spectrum with a Dirac delta peak at $k_p$. This essentially leads to an almost monochromatic9 PBH mass function. When the peak in the scalar spectrum enters the horizon, those fluctuations which are above ${\delta }_c$ collapse into a PBH. By the end of the collapse, only a fraction $\gamma$ of the total mass enclosed inside the horizon at horizon crossing ($H_{⊛}{a}_{⊛}=k_p$) will end up as a PBH. We use the symbol ⊛ to refer to evaluation at PBH formation time. In mathematical terms, we have that PBHs mainly form with a mass equal to the mass enclosed inside the Hubble sphere, that is

$$M_{\mathrm{PBH},⊛}=\gamma \frac{4\pi }{3}\frac{{\rho }_{⊛}}{H_{⊛}^3}=\frac{4\pi \gamma M_{\mathrm{pl}}^2}{H_{⊛}}.$$

(30)

To grasp the magnitude of the actual PBH mass, let us write it in terms of a solar mass $M_{\odot }$ by referring all quantities in (30) in terms of the size of the horizon at matter-radiation equality, which is very well measured by Planck [2]. The numerical factors and formulas used to derive this and the following estimates can be found Appendix A. After some numerical manipulations we arrive at

$$\begin{array}{cc}\hfill M_{\mathrm{PBH},⊛}\approx 3.97\times {10}^{-12}{M}_{\odot }{\left(\frac{k_p}{k_{\mathrm{rh}}}\right)}^{-2-b}& {\left(\frac{k_{\mathrm{rh}}}{{10}^{12}{\mathrm{Mpc}}^{-1}}\right)}^{-2}\hfill \\ & \times \left(\frac{\gamma }{0.2}\right){\left(\frac{g_{*}\left(T_{\mathrm{rh}}\right)}{106.75}\right)}^{1/2}{\left(\frac{g_{*s}\left(T_{\mathrm{rh}}\right)}{106.75}\right)}^{-2/3},\hfill \end{array}$$

(31)

where we kept our ignorance on the primordial universe by leaving a free b (or w) (16). We took $\gamma \sim 0.2$ in a radiation dominated universe [206] but it may differ for different expansion histories [209,210,211,212]. $g_{*}\left(T_{\mathrm{rh}}\right)$ and $g_{*s}\left(T_{\mathrm{rh}}\right)$ are the effective degrees of freedom in the energy density and entropy evaluated at reheating. It should be noted that we also assumed an instantaneous reheating. In that case, the horizon scale at reheating $k_{\mathrm{rh}}$ can be related to the reheating temperature by

$$\begin{array}{c}\hfill k_{\mathrm{rh}}=1.2\times {10}^{12}{\mathrm{Mpc}}^{-1}\left(\frac{T_{\mathrm{rh}}}{5\times {10}^4\mathrm{GeV}}\right){\left(\frac{g_{*}\left(T_{\mathrm{rh}}\right)}{106.75}\right)}^{1/2}{\left(\frac{g_{*s}\left(T_{\mathrm{rh}}\right)}{106.75}\right)}^{-1/3}.\end{array}$$

(32)

As we have seen in Section 2.2, the peak in the primordial curvature power spectrum induces GWs. These GWs have a typical frequency $f_p$ associated to the wavenumber $k_p$. Since the typical mass and the typical frequency are related by the same typical scale $k_p$, we can write one in terms of the other. In this way, the typical frequency evaluated today reads

$$f_p=\frac{k_p}{2\pi a_0}\approx 1.54\times {10}^{-3}\mathrm{Hz}{\left(\frac{M_{\mathrm{PBH},⊛}}{3.97\times {10}^{-12}{M}_{\odot }}\right)}^{\frac{-1}{2+b}}{\left(\frac{f_{\mathrm{rh}}}{1.54\times {10}^{-3}\mathrm{Hz}}\right)}^{\frac{b}{2+b}}.$$

(33)

For a radiation dominated universe ($b=0$) we have a one-to-one correspondence between the frequency and the mass. If we fix the reheating scale, the relation (33) roughly tells us that PBH with a given mass have an induced GW counterpart with a peak frequency at around $f_p$. The exact relation between $f_p$ and $M_{\mathrm{PBH},⊛}$ depends on the value of b. For instance, if $b < 0$ ($w > 1/3$) the expansion rate was faster and for a fixed $k_p$ corresponds to smaller PBH mass at formation. The opposite holds for $b > 0$ ($w < 1/3$). It should also be noted that for $b > 3/2$ ($w < -1/15$) there is another peak in the induced GWs at $k\sim k_{\mathrm{rh}}$. This may break the correspondence between the PBH mass and the peak of the induced GWs if not all of the GW spectrum is seen.

3. General Formalism

In this section we present the general formulation of induced GWs. This includes the derivation of the equations of motion in Section 3.1 and the general form of the solutions in Section 3.2. Then in Section 3.3 we include the leading terms due to local-type non-Gaussianity. In Section 3.1, we derive the equations of motion using the action formalism and so we do not follow the conventional derivation, e.g., of Refs. [35,38]. The final result is obviously the same, but the derivation is rather quick and clear. Throughout this section and on, we shall neglect the effect of vector perturbations as they typically decay. Anyone interested in the approximations and applications may jump directly to Section 4.

Main references: In Section 3.1 we closely follow Ref. [150] although we work at the level of the Lagrangian. For alternative derivations starting from the Einstein equations see for example Refs. [35,38,39]. Then, Section 3.2 is extracted from [64] which is built upon and uses the notation of Ref. [63] except for the numerical factor 2 involving $h_{ij}$. Lastly, Section 3.3 on the primordial non-Gaussianity is mainly based on Ref. [81] updated and corrected by Ref. [92].

3.1. Derivation from the Action

We need to arrive at the equations of motion for the transverse-traceless component of the spatial metric at the second order in cosmological perturbation theory. Although it is a second order calculation, we shall derive it quickly from the action formalism; if we work in a particular gauge and we only focus on the interaction of tensor with scalars. It is particularly convenient to work in a gauge similar to the Newtonian gauge10 with the exponential notation typical of Misner, Throne and Wheeler [184]. With this choice and assuming a perturbed FLRW universe, the metric reads

$$d{s}^2=g_{\mu \nu }d{x}^{\mu }d{x}^{\nu }=-{\mathrm{e}}^{2\Psi }d{t}^2+a^2\left(t\right){\mathrm{e}}^{2\Phi }{\gamma }_{ij}d{x}^{i}d{x}^j,$$

(34)

where $g_{\mu \nu }$ is the metric of our space-time, $i=\{1,2,3\}$ are the spatial components, $a\left(t\right)$ is the scale factor and we have further used the conformal decomposition of the spatial metric such that

$$\frac{\partial }{\partial t}det\gamma ={\gamma }^{ij}\frac{\partial }{\partial t}{\gamma }_{ij}=0.$$

(35)

Note that for a flat FLRW background in Cartesian coordinates we have ${\gamma }_{ij}={\delta }_{ij}$. The conformal decomposition of the spatial metric is very convenient and simplifies calculations. It is also used, e.g., in Numerical Relativity [222]. In short, the conformal decomposition completely splits the trace (volume changing) from the trace-free (volume preserving) degrees of freedom. In the Newtonian gauge this decomposition directly splits the scalar and tensor modes of the spatial metric. This means that ${\gamma }_{ij}$ only contains the transverse-traceless degrees of freedom. As we shall see, it also leaves a clean splitting between ${\gamma }_{ij}$ and $\Phi$ at the level of the action.

To consider a sensible cosmological set up, we include a canonical scalar field $\phi$ in the perturbed FLRW metric. We can later generalize it to a perfect fluid quite straightforwardly. The action, without any decomposition, is given by

$$S=\int d^{4}x\sqrt{-g}\left\{\frac{1}{2}R-\frac{1}{2}{g}^{\mu \nu }{\partial }_{\mu }\phi {\partial }_{\nu }\phi -V\left(\phi \right)\right\},$$

(36)

where g is the determinant of $g_{\mu \nu }$, R is the 4D Ricci scalar, ${\partial }_{\mu }\equiv \partial /\partial x^{\mu }$ and $V\left(\phi \right)$ is the potential of $\phi$. In the (3 + 1) conformal decomposition (34), after some algebra and integration by parts, the action becomes

$$\begin{array}{cc}\hfill S=\int d^{3}xdt& \left\{a{\mathrm{e}}^{\Psi +\Phi }\right(\frac{1}{2}{R}^{\left(3\right)}\left[{\gamma }_{ij}\right]-2{\gamma }^{ij}{D}_i{D}_j\Phi -{\gamma }^{ij}{D}_i\Phi D_j\Phi -\frac{1}{2}{\gamma }^{ij}{D}_i\phi D_j\phi \hfill \\ & +a^3{\mathrm{e}}^{3\Phi -\Psi }\left(\frac{1}{8}{\gamma }^{ij}{\gamma }^{kl}{\dot{\gamma }}_{ik}{\dot{\gamma }}_{jl}-3{\left(H+\dot{\Phi }\right)}^2+\frac{1}{2}{\dot{\phi }}^2\right)-a^3{\mathrm{e}}^{3\Phi +\Psi }V\left(\phi \right)\},\hfill \end{array}$$

(37)

where $R^{\left(3\right)}\left[{\gamma }_{ij}\right]$ and $D_i$ are respectively the 3D Ricci scalar and the covariant derivative associated to ${\gamma }_{ij}$. Since we work in Cartesian coordinates, we already used that $det\gamma =1$. In going from (36) to (37) we took a big leap in the (3 + 1) and conformal decomposition of the 4D Ricci scalar. Some steps can be found in Appendix C. For more details, the interested reader is referred to E. Poisson’s book [223] for the (3 + 1) or ADM decomposition [224], and to the appendix of R. Wald’s book [225] for the conformal transformation rules.

Now, we could take the variation of (37) with respect to ${\gamma }_{ij}$, having in mind that the result should still be transverse and traceless, and obtain the transverse-traceless spatial component of Einstein Equations. However, this would not be very illuminating. Before taking the variation, let us use cosmological perturbation theory and expand the action. A smart way to decompose the spatial metric and to keep the requirement (35) is to take the exponential matrix, i.e.,

$$\begin{array}{c}\hfill {\gamma }_{ij}={\left(e^h\right)}_{ij}={\delta }_{ij}+h_{ij}+\frac{1}{2}{\delta }^{kl}{h}_{ik}{h}_{jl}+O\left(h^3\right),\end{array}$$

(38)

where $h_{ij}\ll 1$ are the transverse-traceless (tensor) modes and satisfy ${\delta }^{ij}{h}_{ij}={\partial }^i{h}_{ij}=0$. This is used, e.g., in Maldacena’s work on non-gaussianities [226]. For the details of the expansion in general gauges see [150]. From now on, spatial indices are contracted with the background spatial metric ${\delta }_{ij}$, for example $h^{ij}{h}_{ij}\equiv {\delta }^{ij}{\delta }^{kl}{h}_{ik}{h}_{jl}$. The inverse metric is then ${\gamma }^{ij}={\delta }^{ik}{\delta }^{jl}{\left(e^{-h}\right)}_{kl}$. With this expansion, the 3D Ricci scalar up to second order is given by

$$R^{\left(3\right)}\left[e^h\right]=-\frac{1}{4}{\partial }_i{h}_{kl}{\partial }^i{h}^{kl}+O\left(h^3\right).$$

(39)

We stopped at second order since we are only interested in scalar-scalar-tensor interactions. We split the other variables as

$$\begin{array}{c}\hfill \Psi =\Psi (t,\mathbf{x}),\Phi =\Phi (t,\mathbf{x}),\phi =\overline{\phi }\left(t\right)+\delta \phi (t,\mathbf{x}),\end{array}$$

(40)

where $\overline{\phi }\left(t\right)$ is the background solution and $\Phi$, $\Psi$ and $\delta \phi$ are the perturbations. Using the perturbative expansions (38) and (40) and only selecting the terms with two scalars and one tensor, we arrive at

$$\begin{array}{cc}\hfill S=\int d^{3}xdt& \left\{\frac{a^3}{8}{\dot{h}}^{ij}{\dot{h}}_{ij}-\frac{a}{8}{\partial }_i{h}_{kl}{\partial }^i{h}^{kl}-2a{h}^{ij}{\partial }_i(\Phi +\Psi ){\partial }_j\Phi +a{h}^{ij}{\partial }_i\Phi {\partial }_j\Phi +\frac{a}{2}{h}^{ij}{\partial }_i\delta \phi {\partial }_j\delta \phi \right\},\hfill \end{array}$$

(41)

where the first two terms correspond to the second order Lagrangian of the tensor modes. By taking the variation with respect to $h_{ij}$ we obtain the equations of motion of induced GWs, namely

$${\ddot{h}}_{ij}+3H{\dot{h}}_{ij}-a^{-2}\Delta h_{ij}=a^{-2}{\mathcal{P}}^{ab}{}_{ij}\left\{-8{\partial }_a(\Phi +\Psi ){\partial }_b\Phi +4{\partial }_a\Phi {\partial }_b\Phi +2{\partial }_a\delta \phi {\partial }_b\delta \phi \right\},$$

(42)

where $H=\dot{a}/a$ is the Hubble parameter and ${\mathcal{P}}^{ab}{}_{ij}$ is the transverse-traceless projector, so that the equation is consistent with the transverse-traceless degrees of freedom. For the moment, it is enough to know that it satisfies the requirements of a transverse-traceless object for both pairs of indices. Also note that in the current case of study, we have no source of anisotropic stress. Thus, for the sake of simplicity and consistency, we momentarily cheat and use the traceless part of the spatial component of the linear Einstein equations in Appendix E, which in the absence of anisotropic stress yields

$$\Psi +\Phi =0.$$

(43)

Note that this is not strictly true in the early universe where neutrinos have a large mean free path and give a tiny contribution to the anisotropic stress [39,40].

Before going to the general solutions of the induced GWs, let us translate the current calculation to the perfect fluid picture. A perfect fluid with energy density $\rho$ and pressure P is described by the following energy-momentum tensor:

$$T_{\mu \nu }=\left(\rho +P\right)u_{\mu }{u}_{\nu }+P{g}_{\mu \nu },$$

(44)

where $u^{\mu }$ is the fluid 4-velocity. The perfect fluid is specified once the equation of state $w=P/\rho$ is given. In the perturbative expansion, one takes $\rho =\overline{\rho }+\delta \rho$ and $u_i={\partial }_{i}v$ where v is the velocity perturbation and we neglected vector modes. The energy momentum tensor (44) has to be compared with the one for the scalar field, which reads

$$T_{\mu \nu }^{\phi }={\partial }_{\mu }\phi {\partial }_{\nu }\phi -g_{\mu \nu }\left(\frac{1}{2}{\partial }_{\alpha }\phi {\partial }^{\alpha }\phi +V\left(\phi \right)\right).$$

(45)

The perfect fluid description of the scalar field follows from the identification

$$u_{\mu }=\frac{{\partial }_{\mu }\phi }{\sqrt{-{\partial }_{\alpha }\phi {\partial }^{\alpha }\phi }}.$$

(46)

If we simply focus on the spatial component $u_i$ we find that

$$\delta \phi \leftrightarrow v\sqrt{\rho +P}.$$

(47)

Thus, in terms of the perfect fluid description we have that

$${\ddot{h}}_{ij}+3H{\dot{h}}_{ij}-a^{-2}\Delta h_{ij}=a^{-2}{\mathcal{P}}^{ab}{}_{ij}\left\{4{\partial }_a\Phi {\partial }_b\Phi +2(\rho +P){\partial }_{a}v{\partial }_{b}v\right\}.$$

(48)

This result coincides with the equations derived in Ref. [38,39]. Note that Equation (48) will slightly change in the case of a multi-perfect fluid system. We expect contributions from the relative velocities of the fluids. For a radiation-matter dominated universe, the source term to induced GWs can be found, e.g., in Ref. [70].

3.2. General Solutions

To solve for the induced GWs, we must solve beforehand the first order equations of motion. Continuing with the scalar field fluctuations $\delta \phi$, we find that they are related by the momentum constraint to $\Phi$ by [183]

$$\delta \phi =-\sqrt{\frac{2}{ϵ}}\left(\Phi +\frac{{\Phi }^{\prime }}{\mathcal{H}}\right).$$

(49)

Thus, we must only solve for the only scalar degree of freedom $\Phi$, which we will do briefly. For the moment, we take a general approach and split $\Phi$ into an initial value ${\Phi }_{\mathbf{k}}$ and a transfer function $T_{\Phi }$ as

$$\Phi (\mathbf{k},\tau )=T_{\Phi }(k,\tau ){\Phi }_{\mathbf{k}}.$$

(50)

If we are considering primordial fluctuations, the initial spectrum for ${\Phi }_{\mathbf{k}}$ is set on superhorizon scales ($k\tau \ll 1$) by quantum fluctuations during inflation. However, in general, any source of fluctuations to the curvature perturbation $\Phi$ will lead to induced GWs. We will discuss this possibility later in Section 6. With Equations (49) and (50) we can formally solve the equations of motion for induced GWs [38,39], which in Fourier space read

$$h_{\mathbf{k},\lambda }^{″}+2\mathcal{H}{h}_{\mathbf{k},\lambda }^{\prime }+k^2{h}_{\mathbf{k},\lambda }={\mathcal{S}}_{\mathbf{k},\lambda },$$

(51)

with

$${\mathcal{S}}_{\mathbf{k},\lambda }=4\int \frac{d^{3}q}{{\left(2\pi \right)}^3}{e}_{\lambda }^{ij}\left(k\right)q_i{q}_j{\Phi }_{\mathbf{q}}{\Phi }_{|\mathbf{k}-\mathbf{q}|}f(\tau ,q,|\mathbf{k}-\mathbf{q}\left|\right),$$

(52)

and

$$\begin{array}{cc}\hfill f(\tau ,q,|\mathbf{k}-\mathbf{q}\left|\right)=& T_{\Phi }\left(q\tau \right)T_{\Phi }\left(\right|\mathbf{k}-\mathbf{q}\left|\tau \right)\hfill \\ & +\frac{1+b}{2+b}\left(T_{\Phi }\left(q\tau \right)+\frac{T_{\Phi }^{\prime }\left(q\tau \right)}{\mathcal{H}}\right)\left(T_{\Phi }\left(\right|\mathbf{k}-\mathbf{q}\left|\tau \right)+\frac{T_{\Phi }^{\prime }\left(\right|\mathbf{k}-\mathbf{q}\left|\tau \right)}{\mathcal{H}}\right),\hfill \end{array}$$

(53)

where we wrote $ϵ$ in terms of b using Equations (15) and (16). To derive Equation (51) from Equation (48), we used the Fourier transform for $h_{ij}$ in terms of the polarization tensors $e_{ij}\left(k\right)$ which are defined in Appendix D. Then, we used that the projection of $e_{ij}$ with the Fourier transform of the transverse-traceless projector ${\tilde{\mathcal{P}}}^{ab}{}_{ij}$ is trivial, i.e., $e^{ij}{\tilde{\mathcal{P}}}^{ab}{}_{ij}=e^{ab}$. Now, with the Green’s method we can find a formal solution to (51) given by

$$h_{k,\lambda }\left(\tau \right)={\int }_{{\tau }_i}^{\tau }d\tilde{\tau }{G}_h(\tau ,\tilde{\tau }){\mathcal{S}}_{\mathbf{k},\lambda }\left(\tilde{\tau }\right),$$

(54)

where $G_h(\tau ,\tilde{\tau })$ is the Green’s function of the homogeneous solutions to (51) given in Appendix B. The formal solution (54) is defined such that at the initial time ${\tau }_i$ we have $h_{\mathbf{k},\lambda }\left({\tau }_i\right)=h_{\mathbf{k},\lambda }^{\prime }\left({\tau }_i\right)=0$. Any GW with primordial origin, i.e., generated during inflation, may be simply added to the solution (54).

Since we are interested in the main observable of the SGWB, that is the GW power, let us compute the 2-point correlation function of the induced GWs which is given by

$$\begin{array}{c}\hfill \langle h_{\lambda }(k,\tau )h_{\lambda }(k^{\prime },\tau )\rangle ={\int }_0^{\tau }d{\tau }_1{\int }_0^{\tau }d{\tau }_{2}G(\tau ,{\tau }_1)G(\tau ,{\tau }_2)\langle S_{\lambda }(k,{\tau }_1)S_{\lambda }(k^{\prime },{\tau }_2)\rangle .\end{array}$$

(55)

Here we have neglected any primordial signal but we could add one without further issue. In fact, this was considered in Ref. [88] where they may have a non-trivial cross-correlated spectrum by tensor-scalar-scalar couplings during inflation. We will not pursue this possibility further in this review. Continuing with Equation (55), we can compute the two-point function of the source term by

$$\begin{array}{cc}\hfill \langle S_{\lambda }(k,{\tau }_1)S_{\lambda }(k^{\prime },{\tau }_2)\rangle =& 4^2\int \frac{d^{3}q}{{\left(2\pi \right)}^3}\int \frac{d^3{q}^{\prime }}{{\left(2\pi \right)}^3}{e}_{\lambda }^{ij}\left(k\right)q_i{q}_j{e}_{\lambda }^{ij}\left(k^{\prime }\right)q_i^{\prime }{q}_j^{\prime }\hfill \\ & \times f({\tau }_1,q,|\mathbf{k}-\mathbf{q}|\left)f\right({\tau }_2,q^{\prime },|k^{\prime }-{\mathbf{q}}^{\prime }\left|\right)\langle {\Phi }_q{\Phi }_{|k-\mathbf{q}|}{\Phi }_{q^{\prime }}{\Phi }_{|k^{\prime }-{\mathbf{q}}^{\prime }|}\rangle .\hfill \end{array}$$

(56)

We see that the induced GW spectrum depends on the 4-point function of the scalar fluctuations. Such a 4-point function may be decomposed in general into a disconnected, i.e., the product of 2-point functions, and a connected piece [89]. If we consider that the stochastic fluctuations are drawn from a Gaussian distribution, the connected 4-point function vanishes and we are led to11

$$\begin{array}{cc}\hfill \langle {\Phi }_q{\Phi }_{|\mathbf{k}-\mathbf{q}|}{\Phi }_{q^{\prime }}{\Phi }_{|{\mathbf{k}}^{\prime }-{\mathbf{q}}^{\prime }|}\rangle & =\frac{2{\pi }^2}{q^3}{\mathcal{P}}_{\Phi }\left(q\right)\frac{2{\pi }^2}{{|\mathbf{k}-\mathbf{q}|}^3}{\mathcal{P}}_{\Phi }\left(\right|\mathbf{k}-\mathbf{q}\left|\right)\hfill \\ & \times {\left(2\pi \right)}^6{\delta }^{\left(3\right)}(\mathbf{q}+{\mathbf{q}}^{\prime }){\delta }^{\left(3\right)}(\mathbf{k}+{\mathbf{k}}^{\prime }-\mathbf{q}-{\mathbf{q}}^{\prime })+(\mathbf{q}\leftrightarrow \mathbf{k}-\mathbf{q}).\hfill \end{array}$$

(57)

We discuss the case of mildly non-Gaussian fluctuations in Section 3.3. With further manipulations of Equations (55) and (8), integrating one internal momentum using a Dirac delta and writing the remaining integral in spherical (momentum) coordinates, we arrive at

$$\overline{{\mathcal{P}}_h}=8{\int }_0^{\infty }dv{\int }_{|1-v|}^{1+v}du{\left(\frac{4v^2-{(1-u^2+v^2)}^2}{4uv}\right)}^2\overline{I^2(\tau ,k,v,u)}{\mathcal{P}}_{\Phi }\left(ku\right){\mathcal{P}}_{\Phi }\left(kv\right),$$

(58)

where we followed the notation of Kohri and Terada [63]. In Equation (58) we already summed over polarizations (see Appendix D for the details) and we have introduced for convenience

$$v\equiv \frac{q}{k},u\equiv \frac{|\mathbf{k}-\mathbf{q}|}{k}.$$

(59)

We also took the oscillation average, i.e., we integrated over half period and divided by $\pi$, to take into account that observations of the SGWB actually measure an average over many wavelengths. Furthermore, we have inserted all the time dependence into a “kernel” or “transfer function” defined by

$$I(\tau ,k,u,v)\equiv {\int }_{{\tau }_i}^{\tau }d\tilde{\tau }G(\tau ,\tilde{\tau })f(\tilde{\tau },k,u,v).$$

(60)

The (averaged) power spectrum (58) is the main quantity needed for the calculations of induced GWs. Knowing (58) we can compute the spectral density of GWs by Equation (7). Note that the induced tensor power spectrum (58) is valid for any curvature perturbation regardless of its origin as long as it is Gaussian. In the next section, we discuss in more detail the case when the fluctuations $\Phi$ are generated during inflation.

3.3. Inclusion of Primordial Non-Gaussianity

Primordial fluctuations generated out of quantum fluctuations during inflation are very close to being Gaussian, although not exactly. Small departures of a Gaussian distribution are expected due to gravitational interactions [183,226]. These perturbative departures from the Gaussian distribution are often referred to as Non-Gaussianity (NG).12 Depending on the formation mechanism, e.g., due to large interactions among fields, we may obtain some significant level of NG. The possibility of observing NG of primordial fluctuations is very exciting, as it might provide information on the particle content during inflation. For this reason, the study of NG is sometimes referred to as cosmological collider physics [227].

For practical convenience, the predictions of primordial fluctuations generated during inflation are often given in terms of the curvature perturbation $\mathcal{R}$. This is because $\mathcal{R}$ is a non-linearly conserved quantity on superhorizon scales [228]. This makes it easy to set well-defined initial conditions on superhorizon scales for whatever scalar variable one might want to consider, independently of what happened previously (e.g., how inflation ended). For instance, since we are working in the Newtonian gauge, we have to relate $\mathcal{R}$ with $\Phi$. It can be shown that on superhorizon scales the linear relation is given by [229]

$$\mathcal{R}=\frac{5+3w}{3(1+w)}\Phi =\frac{2b+3}{b+2}\Phi .$$

(61)

In the case at hand, we are interested on NGs generated during inflation, i.e., primordial NG. These primordial NGs are often parametrized by the amplitude and shape of 3-point function or the bispectrum [180,181]. For simplicity, we usually consider only the local type NG which can be expressed as a local perturbative expansion around the Gaussian curvature perturbation ${\mathcal{R}}^g$ by

$$\mathcal{R}\left(x\right)={\mathcal{R}}^g\left(x\right)+\frac{3}{5}{f}_{NL}{\left({\mathcal{R}}^g\left(x\right)\right)}^2,$$

(62)

where the factor $3/5$ is by convention.13 Note, however, that this is just a particular shape of the NG. Currently though, there is no study on the impact of other shapes of NG on the induced GW spectrum. Now, if we go to Fourier space the squared term in (62) becomes a convolution. Then, we can write that for the ${\Phi }_q$ modes the perturbative expansion is given by

$${\Phi }_q={\Phi }_q^g+F_{NL}\int \frac{d^{3}l}{{\left(2\pi \right)}^3}{\Phi }_l^g{\Phi }_{|\mathbf{q}-\mathbf{l}|}^g,$$

(63)

where to avoid unnecessary numerical factors we have introduced

$$F_{NL}\equiv \frac{3}{5}\frac{2b+3}{b+2}{f}_{NL}.$$

(64)

We can then expand perturbatively the 4-point function in (56) as

$$\begin{array}{cc}\hfill & \langle {\Phi }_q{\Phi }_{|\mathbf{k}-\mathbf{q}|}{\Phi }_{q^{\prime }}{\Phi }_{|{\mathbf{k}}^{\prime }-{\mathbf{q}}^{\prime }|}\rangle =\langle {\Phi }_q^g{\Phi }_{|\mathbf{k}-\mathbf{q}|}^g{\Phi }_{q^{\prime }}^g{\Phi }_{|{\mathbf{k}}^{\prime }-{\mathbf{q}}^{\prime }|}^g\rangle \hfill \\ & +F_{NL}^2\int \frac{d^{3}l}{{\left(2\pi \right)}^3}\int \frac{d^3{l}^{\prime }}{{\left(2\pi \right)}^3}\langle {\Phi }_q^g{\Phi }_l^g{\Phi }_{|\mathbf{k}-\mathbf{q}-\mathbf{l}|}^g{\Phi }_{q^{\prime }}^g{\Phi }_{l^{\prime }}^g{\Phi }_{|{\mathbf{k}}^{\prime }-{\mathbf{q}}^{\prime }-{\mathbf{l}}^{\prime }|}^g\rangle \hfill \\ & +\left(\right|\mathbf{k}-\mathbf{q}|\leftrightarrow q)+(|{\mathbf{k}}^{\prime }-{\mathbf{q}}^{\prime }|\leftrightarrow q^{\prime })+(q\leftrightarrow q^{\prime };|\mathbf{k}-\mathbf{q}|\leftrightarrow |{\mathbf{k}}^{\prime }-{\mathbf{q}}^{\prime }\left|\right)\hfill \\ & +F_{NL}^2\int \frac{d^{3}l}{{\left(2\pi \right)}^3}\int \frac{d^3{l}^{\prime }}{{\left(2\pi \right)}^3}\langle {\Phi }_l^g{\Phi }_{|\mathbf{q}-\mathbf{l}|}^g{\Phi }_{l^{\prime }}^g{\Phi }_{|\mathbf{k}-\mathbf{q}-{\mathbf{l}}^{\prime }|}^g{\Phi }_{q^{\prime }}^g{\Phi }_{|{\mathbf{k}}^{\prime }-{\mathbf{q}}^{\prime }|}^g\rangle \hfill \\ & +(q\leftrightarrow q^{\prime };|\mathbf{k}-\mathbf{q}|\leftrightarrow |{\mathbf{k}}^{\prime }-{\mathbf{q}}^{\prime }\left|\right)+O\left(F_{NL}^4\right).\hfill \end{array}$$

(65)

The first term on the right hand side of (65) is already computed in Equation (57). If we focus on the second line in the NG contribution of (65), we see that there are 6 possible non-vanishing Wick contractions in the first 6 point function of ${\Phi }^g$ which, together with the 4 remaining possibilities of the third line, makes 24 terms for the first NG contribution. In the fourth line of (65) there are 8 possible non-vanishing contractions, which make 16 possibilities counting the fifth line. The total is then 40 non-vanishing contractions.14 The NG contribution can be further classified into three different terms following the notation of Ref. [92]: the “H”, the “C”, the “Z”. Using Equation (65) in (8) and (55), these NG contributions to the induced tensor spectrum can be respectively written as

$$\begin{array}{cc}\hfill \overline{{\mathcal{P}}_h^{\mathrm{H}}}=2^5{F}_{NL}^2{k}^3\sum _{\lambda }\int \frac{d^{3}q}{2\pi }{\left(e_{\lambda }^{ij}\left(k\right)q_i{q}_j\right)}^2& \overline{I^2(\tau ,q,|\mathbf{k}-\mathbf{q}\left|\right)}\hfill \\ & \times \frac{{\mathcal{P}}_{\Phi }\left(q\right)}{q^3}\int \frac{d^{3}l}{2\pi }\frac{{\mathcal{P}}_{\Phi }\left(l\right)}{l^3}\frac{{\mathcal{P}}_{\Phi }\left(\right|\mathbf{k}-\mathbf{q}-\mathbf{l}\left|\right)}{{|\mathbf{k}-\mathbf{q}-\mathbf{l}|}^3},\hfill \end{array}$$

(66)

$$\begin{array}{cc}\hfill \overline{{\mathcal{P}}_h^{\mathrm{C}}}=2^6{F}_{NL}^2{k}^3\sum _{\lambda }\int \frac{d^{3}q}{2\pi }\int \frac{d^{3}l}{2\pi }{e}_{\lambda }^{ij}\left(k\right)q_i{q}_j& e_{\lambda }^{ij}\left(k\right)l_i{l}_j\overline{I(\tau ,q,|\mathbf{k}-\mathbf{q}\left|\right)I(\tau ,l,|\mathbf{k}-\mathbf{l}\left|\right)}\hfill \\ & \times \frac{{\mathcal{P}}_{\Phi }\left(\right|\mathbf{k}-\mathbf{l}\left|\right)}{{|\mathbf{k}-\mathbf{l}|}^3}\frac{{\mathcal{P}}_{\Phi }\left(l\right)}{l^3}\frac{{\mathcal{P}}_{\Phi }\left(\right|\mathbf{q}-\mathbf{l}\left|\right)}{{|\mathbf{q}-\mathbf{l}|}^3},\hfill \end{array}$$

(67)

and

$$\begin{array}{cc}\hfill \overline{{\mathcal{P}}_h^{\mathrm{Z}}}=2^6{F}_{NL}^2{k}^3\sum _{\lambda }\int \frac{d^{3}q}{2\pi }\int \frac{d^{3}l}{2\pi }{e}_{\lambda }^{ij}\left(k\right)q_i{q}_j& e_{\lambda }^{ij}\left(k\right)l_i{l}_j\overline{I(\tau ,q,|\mathbf{k}-\mathbf{q}\left|\right)I(\tau ,l,|\mathbf{k}-\mathbf{l}\left|\right)}\hfill \\ & \times \frac{{\mathcal{P}}_{\Phi }\left(q\right)}{q^3}\frac{{\mathcal{P}}_{\Phi }\left(l\right)}{l^3}\frac{{\mathcal{P}}_{\Phi }\left(\right|\mathbf{k}-\mathbf{q}-\mathbf{l}\left|\right)}{{|\mathbf{k}-\mathbf{q}-\mathbf{l}|}^3}.\hfill \end{array}$$

(68)

These derived formulas agree with those of Ref. [92]. We will assume that we are in the perturbative regime, so that the NG contributions (66) and (67) should be thought of as a correction to the Gaussian contribution (58). The first contribution (66) is also referred to as “hybrid” [90]. Additionally, note that the second line in (66) can be reabsorbed into a redefinition of the primordial power spectrum to include non-Gaussianities. For instance, we could write

$$\begin{array}{c}\hfill {\mathcal{P}}_{\mathcal{R}}^{NL}\left(q\right)={\mathcal{P}}_{\mathcal{R}}\left(q\right)+F_{NL}^2{q}^3\int \frac{d^{3}l}{2\pi }\frac{{\mathcal{P}}_{\mathcal{R}}\left(l\right)}{l^3}\frac{{\mathcal{P}}_{\mathcal{R}}\left(\right|\mathbf{q}-\mathbf{l}\left|\right)}{{|\mathbf{q}-\mathbf{l}|}^3},\end{array}$$

(69)

which follows from the computation of $\langle \Phi \left(k\right)\Phi \left(k^{\prime }\right)\rangle$ using Equation (63). This is the term considered by Cai, Pi and Sasaki [89] and acts as a local redefinition of the primordial spectrum. With the form of Equation (69) it is much easier to compute the induced GW spectrum up to $O\left(F_{NL}^4\right)$ since it just accounts for a replacement of ${\mathcal{P}}_{\mathcal{R}}$ in (58) for ${\mathcal{P}}_{\mathcal{R}}^{NL}$.The second contribution (67) was first considered by Unal [90] and the third (68) by Ref. [81]. Unal also considered the terms of $O\left(F_{NL}^4\right)$ [90]. Most recently, Adshead, Lozanov and Weiner [92] presented a detailed and complete analysis of all the contributions up to $O\left(F_{NL}^4\right)$. The interested reader is referred to Ref. [92] for the detailed formulas. We discuss the main features of the NG corrections in Section 5. It should be noted that the effect of primordial NG might compete with the terms due to non-linearities, e.g., from studied third order contributions to the source term of induced GWs. The reader is referred to Refs. [143,179] for the details on the third order expansion. Lastly, it should also be noted that for simplicity we assumed a scale invariant $f_{NL}$ in (62). However, in general situations one expects that $f_{NL}$ is scale dependent. For a recent study on the impact of the scale dependence see Ref. [93].

4. Analytical Transfer Functions

In general situations the kernel (60) and the power spectrum (58) have to be computed numerically. The problem is that the time integral in (60) is extremely demanding for $k\tau \gg 1$, as the integrand is the product of three oscillating functions with frequencies that depend on the tensor and scalar mode wavenumbers. In fact, $k\tau \gg 1$ is the range of scales we are interested in. If we extend the integral until today, say at $\tau ={\tau }_0\propto H_0^{-1}$, and we look at scales much smaller than the horizon, we trivially have $k{\tau }_0\gg 1$. However, there is more we can learn. For instance, scales accessible to future GW detectors range from $k\sim {10}^7$ to ${10}^{18}{\mathrm{Mpc}}^{-1}$. These very small scales entered the horizon much before BBN, which roughly correspond to the time when modes with $k\sim {10}^3{\mathrm{Mpc}}^{-1}$ entered the horizon.15 Luckily, this implies that we can evaluate the kernel (60) some time prior to BBN when induced GWs are really inside the horizon and in a radiation dominated universe. Then, from that moment on, we may simply treat these induced GWs as a radiation fluid with $w=1/3$. Note that we will also consider the case that the modes of interest enter the horizon in an epoch where $w\ne 1/3$. In that case, we must follow induced GWs until the universe transitions to the standard radiation dominated era. We will provide more details later.

Surprisingly, there are relevant regimes where the integrals in (60) can be done analytically. This is the case of a constant equation of state parameter w and constant propagation speed of fluctuations $c_s$. For example, this includes a perfect fluid whether it is made of a scalar field16 or an adiabatic perfect fluid, which are commonly used in cosmology. Even when there are transitions between one perfect fluid domination to another, there will be periods in which w and $c_s$ are almost constant. Thus, the constant w and $c_s$ limit is a good approximation in a cosmological set up if the modes of interest enter the horizon outside the transition periods. With these assumptions, we present the solutions to first order cosmological perturbations in the Newtonian gauge in Section 4.1 and derive a very simple form of the source term (53) for the induced GWs. In the general situation where $w\ne 1/3$, the integral (60) can be analytically done for modes which either enter the horizon before or after the transition. In this way, in Section 4.2 and Section 4.3 we respectively derive the kernel for modes which are subhorizon and superhorizon before the transition. When then follow our solutions to the radiation dominated epoch.

Main references: The whole section is a generalisation of Refs. [64,74] to account for not only general constant w but general constant $c_s$. The notation used here improves and unifies those of [64,74]. The main base of Refs. [64,74] was started in [39,63] in the radiation dominated universe.

4.1. First Order Solutions

At the moment, the formal solution (54) for the induced tensor modes does not tell us much. To extract some meaningful information, we have to solve the equations of motion for first order perturbations, which are given in detail in Appendix E. After some simplifications, the equation for the Newtonian potential for general w and $c_s$ and in the absence of isocurvature17 fluctuations reads18

$${\Phi }^{″}+\left(2ϵ-\eta \right)\mathcal{H}{\Phi }^{\prime }-\left(\eta +2s\left(1+ϵ-\eta -2s+\frac{\dot{s}}{Hs}\right)\right){\mathcal{H}}^2\Phi +c_s^2{k}^2\Phi =0,$$

(70)

where we have also introduced $ϵ$ as in Equation (15) as well as

$$\eta \equiv \frac{\dot{ϵ}}{Hϵ}\mathrm{and}s\equiv \frac{{\dot{c}}_s}{H{c}_s}.$$

(71)

The notation $ϵ$, $\eta$ and s is typical of inflationary models and simplifies the form of the equations. One may of course write Equation (70) in a more conventional notation, e.g., that in Mukhanov’s book [229] which for constant $c_s$ reads

$${\Phi }^{″}+3\mathcal{H}\left(1+c_w^2\right){\Phi }^{\prime }+\left(2{\mathcal{H}}^{\prime }+(1+3c_w^2){\mathcal{H}}^2\right)\Phi +c_s^2{k}^2\Phi =0,$$

(72)

where $c_w^2=\dot{P}/\dot{\rho }$ and $c_s^2=\delta P/\delta \rho$. In the case of an adiabatic perfect fluid one has $c_s^2=c_w^2=w$. In the case of a canonical scalar field we have $c_w^2=w$ and $c_s^2=1$. As we already mentioned, for analytical viability we are mostly interested in the case of constant w and $c_s$. This requirement in turn implies $\eta =s=0$, so that the Newtonian potential in Equation (70) behaves as a massless field. The solution to (70) for general $c_s$ and w is given in terms of Bessel functions of the first and second kind, explicitly

$$\Phi \left(k\tau \right)={\left(c_{s}k\tau \right)}^{-b-3/2}\left(C_1{J}_{b+3/2}\left(c_{s}k\tau \right)+C_2{\gamma }_{b+3/2}\left(c_{s}k\tau \right)\right),$$

(73)

where b in terms of w is defined in Equation (16). Now, we must give some kind of initial conditions to $\Phi$. If the primordial spectrum was generated by quantum fluctuations during inflation, the initial conditions are well-defined and constant on super(sound)horizon scales, that is when $c_{s}k\tau \ll 1$. By picking up the asymptotically constant term in (73), we find that $C_2=0$ and

$$\Phi \left(k\tau \right)={\Phi }_k{2}^{b+3/2}\Gamma [b+5/2]{\left(c_{s}k\tau \right)}^{-b-3/2}{J}_{b+3/2}\left(c_{s}k\tau \right),$$

(74)

where ${\Phi }_k$ is the primordial (stochastic) value set by inflation. Note that the general solution (73) is not valid for $c_s=0$ where the gradient term in (70) is absent. Nevertheless, the solution for constant w is just $\Phi ={\Phi }_{\mathbf{k}}$ on all scales. The case of $c_s=0$ is quite particular as the fact that $\Phi$ does not decay implies a constant source to induced GWs. This makes the detailed dynamics of the transition to radiation domination very important for the final GW spectrum [65,66], when the source term experiences drastic changes. We dedicate Section 6 to this particular case.

With an exact analytical formula for the first order perturbations, we can attempt to compute the resulting induced GWs. For that, we first need to know the two independent homogeneous solutions to the tensor mode Equation (51), say $h_1$ and $h_2$. For constant w they are given as well in terms of Bessel functions with one order less than in $\Phi$, concretely

$$h_1\left(k\tau \right)={\left(k\tau \right)}^{-b-1/2}{J}_{b+1/2}\left(k\tau \right)\mathrm{and}{h}_2\left(k\tau \right)={\left(k\tau \right)}^{-b-1/2}{\gamma }_{b+1/2}\left(k\tau \right).$$

(75)

In Equation (75), $h_1$ and $h_2$ respectively are the growing and decaying modes. Note that there is no $c_s$ since in general relativity tensor modes propagate at the speed of light. This yields a Green’s function (see Appendix B) given by

$$G(\tau ,\tilde{\tau })=\frac{\pi }{2k}\frac{{\left(k\tilde{\tau }\right)}^{b+3/2}}{{\left(k\tau \right)}^{b+1/2}}\left(J_{b+1/2}\left(k\tilde{\tau }\right){\gamma }_{b+1/2}\left(k\tau \right)-J_{b+1/2}\left(k\tau \right){\gamma }_{b+1/2}\left(k\tilde{\tau }\right)\right).$$

(76)

This finishes the list of necessary ingredients for our induced GWs calculations. In passing, since there always appears the combination $k\tau$ or $k\tilde{\tau }$ we introduce the useful notation

$$x\equiv k\tau ,$$

(77)

which we will use from now on. Additionally, a tilde on x implies a tilde on $\tau$, i.e., $\tilde{x}=k\tilde{\tau }$.

Thus far it seems that we may have to deal with several integrals containing at least three Bessel functions, which is challenging to do analytically in general. However, the source term to the induced GWs ends up acquiring a simple enough form after several simplifications. For instance, using the Bessel functions properties (see Appendix F), we find that Equation (53) reduces to

$$\begin{array}{cc}\hfill f(x,u,v)=& \frac{2^{2b+3}{\Gamma }^2[b+5/2]}{(2b+3)(b+2)}{\left(c_{s}x\right)}^{-2b-1}{\left(uv\right)}^{-b-1/2}\hfill \\ & \times \left(J_{b+1/2}\left(c_{s}vx\right)J_{b+1/2}\left(c_{s}ux\right)+\frac{b+2}{b+1}{J}_{b+5/2}\left(c_{s}vx\right)J_{b+5/2}\left(c_{s}ux\right)\right),\hfill \end{array}$$

(78)

where we used x from (77) and u and v from (59). To arrive at (78) one needs to write the Bessel function $J_{b+3/2}$ that appears in $\Phi$ in a combination of $J_{b+1/2}$ and $J_{b+5/2}$ which can be found in Appendix F. So far, it is not entirely clear whether there is any physical reason behind the simplification (78) which turns out to render the source term as a sum of the squared of Bessel functions of the same order. This simple form of the source term is crucial to obtain analytical formulas for general b and $c_s$ that we derive below.

Replacing (77) and (78) into Equation (60) we find that the kernel can be written as

$$\begin{array}{cc}\hfill I(x,u,v)=& \pi 4^b{\Gamma }^2[b+3/2]\frac{2b+3}{b+2}{\left(c_s^{2}uvx\right)}^{-b-1/2}\left(J_{b+1/2}\left(x\right){\mathcal{I}}_{\gamma }-{\gamma }_{b+1/2}\left(x\right){\mathcal{I}}_J\right),\hfill \end{array}$$

(79)

where we defined for compactness

$$\begin{array}{cc}\hfill {\mathcal{I}}_{J/\gamma }\equiv {\int }_0^{x}d\tilde{x}& {\tilde{x}}^{1/2-b}\left\{\begin{array}{c}{J}_{b+1/2}\left(\tilde{x}\right)\\ {\gamma }_{b+1/2}\left(\tilde{x}\right)\end{array}\right\}\hfill \\ & \times \left(J_{b+1/2}\left(c_{s}v\tilde{x}\right)J_{b+1/2}\left(c_{s}u\tilde{x}\right)+\frac{b+2}{b+1}{J}_{b+5/2}\left(c_{s}v\tilde{x}\right)J_{b+5/2}\left(c_{s}u\tilde{x}\right)\right).\hfill \end{array}$$

(80)

Unfortunately, we are not aware of a general analytical formula for the time integrals (80) unless $b\in \mathbb{Z}$, in which case the Bessel functions reduce to spherical Bessel functions and can be written in terms of transcendental functions. This is for example the case of radiation domination $b=0$ ($w=1/3$), where the general kernel can be found in Ref. [63]. Nevertheless, we can integrate (80) in two relevant limits, $x\gg 1$ (subhorizon) and $x\ll 1$ (superhorizon).

Let us explain the assumptions needed for the approximations before jumping into the actual calculations. We are considering that right after inflation there is a period when the universe has a free b (or w) and free $c_s$. However, to recover the standard hot big bang cosmology, such an epoch must end, and the universe must enter a radiation domination stage. Thus, whatever solution we find in this section for the induced GWs must be followed until we reach the radiation dominated stage. For simplicity, we will assume an instantaneous transition. Then, we will refer to the time of transition as “reheating” time ${\tau }_{\mathrm{rh}}$ and the comoving size of the horizon at that time as “reheating” scale $k_{\mathrm{rh}}$. The value of $k_{\mathrm{rh}}$ in terms of the “reheating” temperature $T_{\mathrm{rh}}$ can be found in Equation (32). The assumption of an instantaneous transition might be important for those tensor modes with wavenumber $k\sim k_{\mathrm{rh}}$. Namely, those modes which entered the horizon around the transition. If the transition is gradual, then our approximations will not be very good for modes which entered during the last stages of the transition. In any case, as we shall see, the instantaneous approximation gives a very good idea of the shape of the induced GW spectrum even for $k\sim k_{\mathrm{rh}}$. On top of that, we will assume that the scale corresponding to the peak of the primordial spectrum enters the horizon well before reheating so that during and after reheating there is no significant source of induced GWs. In the case of a top-hat primordial spectrum, this means that the scale corresponding to the IR cut-off enters the horizon much before reheating. Under these two assumptions, i.e., sudden transition and “localized” primordial spectrum, we can proceed to derive our analytical formulas. Note that this does not apply to the case when $c_s=0$ which we discuss in Section 6.

4.2. General Subhorizon Kernel

Let us start with the simplest situation. This corresponds to the case when all modes of interest entered the horizon much before reheating, that is $k\gg k_{\mathrm{rh}}$. Then, we can safely send the upper integration limit in (80) to infinity, i.e., $x=k\tau \gg 1$. This is because when we eventually match our solution at reheating to a free GW propagating in a radiation dominated universe, the upper limit would be $x_{\mathrm{rh}}=k{\tau }_{\mathrm{rh}}\sim k/k_{\mathrm{rh}}\gg 1$. The integrals (80) for $x\to \infty$ are derived by Gervois and Navelet [233] and given in Appendix G. Here we just give the form of (79) after integration, which reads

$$\begin{array}{cc}\hfill I(x\gg 1,u,v)=& x^{-b-1}{4}^b{\Gamma }^2[b+3/2]\frac{2b+3}{b+2}\frac{|1-y^2{|}^{b/2}}{c_s^{2}uv}\hfill \\ & \times \{cos\left(x-\frac{b\pi }{2}\right)\left({\mathsf{P}}_b^{-b}\left(y\right)+\frac{b+2}{b+1}{\mathsf{P}}_{b+2}^{-b}\left(y\right)\right)\Theta [c_s(u+v)-1]\hfill \\ & +\frac{2}{\pi }sin\left(x-\frac{b\pi }{2}\right)\left({\mathsf{Q}}_b^{-b}\left(y\right)+\frac{b+2}{b+1}{\mathsf{Q}}_{b+2}^{-b}\left(y\right)\right)\Theta [c_s(u+v)-1]\hfill \\ & -\frac{2}{\pi }sin\left(x-\frac{b\pi }{2}\right)\left({\mathcal{Q}}_b^{-b}(-y)+2\frac{b+2}{b+1}{\mathcal{Q}}_{b+2}^{-b}(-y)\right)\Theta [1-c_s(u+v)]\},\hfill \end{array}$$

(81)

where $\Theta \left[x\right]$ is the Heaviside function, ${\mathsf{P}}_b^{-b}\left(y\right)$ and ${\mathsf{Q}}_b^{-b}\left(y\right)$ are the Ferrer’s function of the first and second kind, valid for $\left|y\right| < 1$, and ${\mathcal{Q}}_b^{-b}\left(y\right)$ is the associated Legendre function of the second kind, valid for $\left|y\right| > 1$. Their explicit expressions in terms of Hypergeometric functions can be found in Appendix H. We also defined for simplicity

$$\gamma =1-\frac{1-c_s^2{(u-v)}^2}{2c_s^{2}uv}=-1-\frac{1-c_s^2{(u+v)}^2}{2c_s^{2}uv},$$

(82)

and we Taylor expanded the Bessel functions for large argument. Note that y is related to the area of the triangle made by the three momenta, $c_s|\mathbf{k}-\mathbf{q}|$, $c_s\left|\mathbf{q}\right|$ and k. The presence of the Heaviside function in (81) signals the possible resonance we anticipated in Section 2.2 when the wavenumber of a tensor modes equals the sum of two typical scalar modes $2c_s{k}_p\sim k$. In more practical terms, when $c_s(u+v)\sim 1$ the three Bessel functions in the definite integral (80) for $x\to \infty$ conspire (by rendering their product as coherent oscillations) and yield a divergent integral. Additionally, note that Equation (81) has the same time dependence as a subhorizon tensor mode propagating in a constant b and $c_s$ universe, e.g., see Equation (75). In other words, it decays as $a^{-2}$ and oscillates with a frequency proportional to x. This was expected as the source term for induced GWs decays once inside the horizon and, therefore, induced GWs behave as a freely propagating GW. This does not apply in the case of a dust dominated universe though.

In the current form Equation (81) is general enough so that it can also be used in the C (67) and Z (68) NG terms. However, in this review we are more interested with the Gaussian contribution (58). Thus, we are more concerned about the averaged square kernel. After some simple algebra, that is squaring (81) and calculating the oscillation average, we arrive at

$$\begin{array}{cc}\hfill \overline{I^2(x,u,v)}=& x^{-2(b+1)}{4}^{2b}{\Gamma }^4[b+3/2]{\left(\frac{2b+3}{b+2}\right)}^2\frac{|1-y^2{|}^b}{2c_s^4{u}^2{v}^2}\hfill \\ & \times \{{\left({\mathsf{P}}_b^{-b}\left(y\right)+\frac{b+2}{b+1}{\mathsf{P}}_{b+2}^{-b}\left(y\right)\right)}^2\Theta [c_s(u+v)-1]\hfill \\ & +\frac{4}{{\pi }^2}{\left({\mathsf{Q}}_b^{-b}\left(y\right)+\frac{b+2}{b+1}{\mathsf{Q}}_{b+2}^{-b}\left(y\right)\right)}^2\Theta [c_s(u+v)-1]\hfill \\ & +\frac{4}{{\pi }^2}{\left({\mathcal{Q}}_b^{-b}(-y)+2\frac{b+2}{b+1}{\mathcal{Q}}_{b+2}^{-b}(-y)\right)}^2\Theta [1-c_s(u+v)]\}.\hfill \end{array}$$

(83)

With the averaged kernel (83) we are ready to compute the induced GWs in quite general situations. We are only left with the integral over momenta which can be easily computed numerically. Analytical formulas for the final GW spectrum are discussed in Section 5. Let us note that although (83) might look technically complicated to implement, it is by far more useful than attempting a numerical integration of three Bessel functions. Furthermore, the superhorizon approximation is very good for modes with $k\gg k_{\mathrm{rh}}$ regardless of the shape of the primordial spectrum. This is because induced GWs with $k\gg k_{\mathrm{rh}}$ are essentially already free GWs. Thus, soon after horizon re-entry and thereon they will behave as relativistic particles in an expanding universe. This also means that we can evaluate (7) and (83) at the instantaneous reheating time, i.e., $\tau ={\tau }_{\mathrm{rh}}$, to estimate the observable spectral density of induced GWs. We will do this shortly. If the transition is gradual, the kernel (83) is a good approximation for modes which enter during a period with almost constant b.

In the next subsections, we will have a closer inspection to the kernel (83). After matching to radiation domination, we will look at the behaviour of the kernel near the resonant point and in the IR tail.

4.2.1. Matching to Radiation Domination

Let us discuss the matching to radiation domination in more technical terms. After matching the background quantities at $\tau ={\tau }_{\mathrm{rh}}$, that is the scale factor a and the Hubble parameter $\mathcal{H}$ we find that a and $\mathcal{H}$ in the radiation dominated period follows the typical solution (found in (16) with $w=1/3$) with a shifted conformal time, which reads

$$\begin{array}{c}\hfill \overline{\tau }\equiv \tau -\frac{b}{1+b}{\tau }_{\mathrm{rh}},\end{array}$$

(84)

where b is related to the equation of state of the previous epoch by (16). Now, for the matching of tensor modes it is important to realize that the continuity of $h_{ij}$ and its first derivative imply the continuity of the kernel (60) and its first derivative. Thus, from now on we simply deal with the matching of the kernel (60). Since we are dealing with subhorizon scales we have that the kernel before the transition is

$$\begin{array}{cc}\hfill I(x\gg 1,u,v)\approx & x^{-b-1}\left(A_{1,b}sin\left(x-\frac{b\pi }{2}\right)+A_{2,b}cos\left(x-\frac{b\pi }{2}\right)\right),\hfill \end{array}$$

(85)

where the coefficients, not important at the moment, can be extracted from Equation (81). The kernel after the transition is also for subhorizon modes, which must have the form of a freely propagating GW in a radiation dominated universe, namely

$$\begin{array}{cc}\hfill I_{RD}(x\gg 1,u,v)\approx & A_{1,RD}\frac{sin\overline{x}}{\overline{x}}+A_{2,RD}\frac{cos\overline{x}}{\overline{x}}.\hfill \end{array}$$

(86)

It is easy to convince oneself that after matching (85) and (86) at $\tau ={\tau }_{\mathrm{rh}}$ one finds that

$$\begin{array}{c}\hfill A_{1,RD}^2+A_{2,RD}^2={\left(\frac{k{\tau }_{\mathrm{rh}}}{1+b}\right)}^2\left(A_{1,b}^2+A_{2,b}^2\right).\end{array}$$

(87)

This implies that although we should be calculating

$${\Omega }_{\mathrm{GW}}(k\gg k_{\mathrm{rh}},\tau \gg {\tau }_{\mathrm{rh}})=\frac{k^2}{12a^2{H}^2}\overline{{\mathcal{P}}_h^{RD}},$$

(88)

with

$$\overline{{\mathcal{P}}_h^{RD}}=8{\int }_0^{\infty }dv{\int }_{|1-v|}^{1+v}du{\left(\frac{4v^2-{(1-u^2+v^2)}^2}{4uv}\right)}^2\overline{I_{RD}^2(x\gg 1,v,u)}{\mathcal{P}}_{\Phi }\left(ku\right){\mathcal{P}}_{\Phi }\left(kv\right),$$

(89)

it is equivalent to evaluate

$${\Omega }_{\mathrm{GW}}(k\gg k_{\mathrm{rh}},\tau \gg {\tau }_{\mathrm{rh}})=\frac{k^2}{12a^2{H}^2}\overline{{\mathcal{P}}_h}{|}_{\tau ={\tau }_{\mathrm{rh}}},$$

(90)

with

$$\overline{{\mathcal{P}}_h}=8{\int }_0^{\infty }dv{\int }_{|1-v|}^{1+v}du{\left(\frac{4v^2-{(1-u^2+v^2)}^2}{4uv}\right)}^2\overline{I^2(x\gg 1,v,u)}{\mathcal{P}}_{\Phi }\left(ku\right){\mathcal{P}}_{\Phi }\left(kv\right),$$

(91)

where the kernel is that given in (83).

It is important to note that our subhorizon approximation is valid for modes with $k{\tau }_{\mathrm{rh}}\gg 1$. Later, we may try to join the subhorizon approximation with the superhorizon approximation valid for $k{\tau }_{\mathrm{rh}}\ll 1$ by extrapolating up to $k{\tau }_{\mathrm{rh}}\sim 1$. Let us call such scale $k_m$. However, there is an important subtlety in joining the two approximations. If we look at the scale that last crossed the horizon at reheating, namely

$$\begin{array}{c}\hfill k_{\mathrm{rh}}={\mathcal{H}}_{\mathrm{rh}}=\frac{1+b}{{\tau }_{\mathrm{rh}}},\end{array}$$

(92)

we see that if we compare with $k_m$ we have that

$$\begin{array}{c}\hfill \frac{k_m}{k_{\mathrm{rh}}}=\frac{1}{1+b}.\end{array}$$

(93)

For $b > 0$ we have that $k_m<k_{\mathrm{rh}}$. This means that our subhorizon approximation is even valid for $k\sim k_{\mathrm{rh}}$ and, therefore, we expect a good matching with the superhorizon approximation at $k\sim k_{\mathrm{rh}}$ (since for $k<k_{\mathrm{rh}}$ the modes are superhorizon before the transition). However, for $b<0$ we find $k_m > k_{\mathrm{rh}}$. This time the subhorizon approximation breaks down for $k > k_{\mathrm{rh}}$ and so we expect the matching with the superhorizon approximation to be at $k\sim k_m$. This is what we eventually find in Section 4. In addition to the above, we have that the subhorizon approximation has corrections of $O\left(x^{1+b}\right)$ [74]. So for $b<0$ the corrections quickly become important for $x\sim 1$. A detailed studied for $b<0$, perhaps with numerical methods, is needed to clarify the exact shape of the GW spectrum. Nevertheless, our results should be a good order of magnitude estimate.

Let us remind the reader that when implementing (83) inside the tensor spectrum (58), we must bear in mind that if we are considering primordial curvature fluctuations from inflation, we have to rescale the spectrum of $\Phi$ accordingly, namely

$${\mathcal{P}}_{\Phi }={\left(\frac{b+2}{2b+3}\right)}^2{\mathcal{P}}_{\mathcal{R}}.$$

(94)

4.2.2. Resonances

Here we focus on the general behaviour near the resonant point for the averaged kernel (83) that we discussed in Section 2.2. We start by noting that we have three vectors, the one from the tensor mode $\mathbf{k}$ and the two scalar modes, $\mathbf{q}$ and $\mathbf{k}-\mathbf{q}$, which of course satisfy momentum conservation. Recall that from Equation (59) q corresponds to v and $|\mathbf{k}-\mathbf{q}|$ to u. This means that our three variables satisfy a triangle inequality, that is

$$\begin{array}{c}\hfill |\mathbf{k}-\mathbf{q}|<|\mathbf{k}-\mathbf{q}|<k+q\mathrm{or}\left|\right|\mathbf{k}-\mathbf{q}|-\mathbf{q}|<k<|\mathbf{k}-\mathbf{q}|+q,\end{array}$$

(95)

which translates into

$$\begin{array}{c}\hfill |1-v|<u<1+v\mathrm{or}|u-v|<1<u+v.\end{array}$$

(96)

The area A of the triangle made by k, q and $|\mathbf{k}-\mathbf{q}|$ is related to the angle between $\mathbf{k}$ and $\mathbf{q}$ by

$$\begin{array}{c}\hfill sin{\theta }_k=\frac{2A}{kq}.\end{array}$$

(97)

This means that the projection with the polarization tensor (see Appendix D) is proportional to the area, namely

$$\begin{array}{c}\hfill e_{\lambda }^{ij}{q}_i{q}_j\propto {sin}^2{\theta }_k\propto A^2.\end{array}$$

(98)

This implies that the saturation of the triangle inequality at $u+v=1$ corresponds to zero area and, therefore, does not contribute to the integral (58). Now, from the arguments in Section 2.2 we expect a resonance when

$$\begin{array}{c}\hfill c_s(u+v)=1.\end{array}$$

(99)

Thus, if we have $c_s^2=1$ the resonant point coincides with the saturation of the triangle inequality and, thus, has vanishing contribution. In this way we see that, by momentum conservation, there is no resonance for $c_s^2=1$ in the induced GW spectrum. Nevertheless, for $c_s^2<1$ we do have a resonance inside the integration region (see Figure 2).

We can confirm these expectations by expanding the kernel (83) around $c_s(u+v)\sim 1$ which corresponds to $\gamma \sim -1$, where the associated Legendre functions might have a divergence. The asymptotic behaviours around $\gamma \sim -1$ are given in Appendix H. Focusing first on the case $b<0$, we find that the kernel diverges as

$$\begin{array}{cc}\hfill \overline{I_{\mathrm{res}}^2(x,u,v,b<0)}\approx & x^{-2(b+1)}{4}^{2b}{\Gamma }^4[b+5/2]\frac{8}{c_s^4{u}^2{v}^2}\frac{{csc}^2\left(b\pi \right)}{{\Gamma }^2[3+b]}{|1-|\gamma \left|\right|}^{2b}.\hfill \end{array}$$

(100)

This recovers Equation (29) which showed that a very peaked primordial spectrum, so that $u\sim v\sim k_p/k$, leads to

$$\begin{array}{c}\hfill {\Omega }_{\mathrm{GWs}}^{\mathrm{res}}(b<0)\propto {|1-|\gamma \left|\right|}^{2b}\propto {|k-2c_s{k}_p|}^{2b}.\end{array}$$

(101)

The case $b=0$ is explicitly given in Equation (209). In the case when $b > 0$, we arrive at

$$\begin{array}{cc}\hfill \overline{I_{\mathrm{res}}^2(x,u,v,b > 0)}\approx & x^{-2(b+1)}{4}^{3b}{\Gamma }^4[b+5/2]\frac{32}{{\pi }^2{c}_s^4{u}^2{v}^2}{\left(\frac{(1+b+b^2)}{(b+1)(b+2)}\frac{\Gamma \left[b\right]}{\Gamma [2b+3]}\right)}^2,\hfill \end{array}$$

(102)

which is independent of y. Thus, as expected there is no resonance for $b > 0$. Nevertheless, as shown in [64] for $1 > b > 0$ there is still a peak at $k=2c_s{k}_p$ in the induced GWs spectrum for a peaked primordial spectrum. We can roughly understand the peak by noting that since ${\mathcal{Q}}_b^{-b}\left(y\right)\sim {(1-y)}^{-b}$, the next leading order must be ${\mathcal{Q}}_b^{-b}\left(y\right)\sim {(1-y)}^{-b}+\mathrm{constant}$ since $b<1$. Then, inserting this into the kernel (83) we find that

$$\begin{array}{c}\hfill {\Omega }_{\mathrm{GWs}}^{\mathrm{res}}(1<b<0)\propto {\mathrm{constant}-|1-|y\left|\right|}^b\propto \mathrm{constant}-{|k-2c_s{k}_p|}^b.\end{array}$$

(103)

For $b > 1$ and $c_s\ne 1$ there is no detailed study of the kernel (83). Nevertheless, by looking at the explicit expression for $b=2$ (213) we see that there is nothing special at $c_s(u+v)=1$. Before ending this section, note that the kernel is symmetric around the resonant point $c_s(u+v)=1$. It displays the same behaviour when approaching from both directions, $c_s(u+v)\to 1^{\pm }$.

4.2.3. Infrared Regime

Here we show the behaviour of the kernel when $u\sim v\gg 1$ which is the relevant regime for the IR tail of the GW spectrum as $v\sim k_p/k$ and $k\ll k_p$. This limit means that, if there is any peak in the spectrum of scalar fluctuations, we are looking at scales much larger than the scale of the peak. Thus, such limit of the kernel (83) is relevant for wavenumbers $k_{\mathrm{rh}}\ll k\ll k_p$, where we included $k_{\mathrm{rh}}$ since it is the limit of the subhorizon approximation. The cases $b<0$ and $b > 0$ exhibit different behaviours as anticipated in Section 2.2. In the limit $y\to 1$, the average kernel (83) is respectively given by

$$\begin{array}{cc}\hfill \overline{I_{\mathrm{IR}}^2(x\gg 1,u,v,b<0)}\approx & x^{-2(b+1)}{4}^{2b}{\Gamma }^4[b+5/2]\frac{8}{c_s^4{u}^2{v}^2}\frac{{csc}^2\left(b\pi \right)}{{\Gamma }^2[3+b]}{|1-y|}^{2b},\hfill \end{array}$$

(104)

for $b<0$ and

$$\begin{array}{cc}\hfill \overline{I_{\mathrm{IR}}^2(x\gg 1,u,v,b > 0)}\approx & x^{-2(b+1)}{4}^{3b}{\Gamma }^4[b+5/2]\frac{32}{{\pi }^2{c}_s^4{u}^2{v}^2}{\left(\frac{(1+b+b^2)}{(b+1)(b+2)}\frac{\Gamma \left[b\right]}{\Gamma [2b+3]}\right)}^2,\hfill \end{array}$$

(105)

for $b > 0$. They are actually very similar to (100) and (102), except for the fact that we are now looking at $u\sim v\gg 1$, so that $1-y\sim v^{-2}$. Evaluating (104) and (105) at the time of reheating and taking into account the factor $k^2$ in (58), we see that whatever the IR scaling in the radiation dominated universe (which for localized sources is $k^3$ [78]) we have to correct the exponent by a factor $-2\left|b\right|$, i.e., to $k^{3-2\left|b\right|}$ for a localized sourced. This is the same result as in Section 2.2. The absolute value is due to an extended superhorizon growth of tensor modes for $b<0$, as explained in Section 2.2.

4.3. Superhorizon Kernel Approximation

In the previous subsections we have studied in detail the subhorizon kernel. However, there are modes which are still superhorizon at the time of reheating. Thus, for such modes we cannot take the upper limit in (80) to $x\to \infty$ as $x\ll 1$ at reheating. In addition to the approximation $x\ll 1$, we will use the assumption that the peak in the scalar spectrum ${\mathcal{P}}_{\Phi }$ is on scales $k_p\gg k_{\mathrm{rh}}$. This means that we are only concerned with the region of integration of (80) around the peak in ${\mathcal{P}}_{\Phi }$, where we have that $v\sim k_p/k\gg 1$ and therefore $u\sim v\gg 1$. Using this fact we can expand the first Bessel function for a small argument in (80) and integrate in the limit of $u\sim v$. By doing so we arrive at Equation (20) but for the kernel, namely

$$\begin{array}{cc}\hfill I(x\ll 1,u,v)\approx & B_{1,b}+B_{2,b}{x}^{-2b},\hfill \end{array}$$

(106)

where the exact coefficients are given by

$$\begin{array}{c}\hfill B_{1,b}=-\frac{{(3+2b)}^2(1+b+b^2)}{4b{(1+b)}^2(2+b)}{\left(c_{s}v\right)}^{-2},B_{2,b}=\frac{4^{1+b}{\Gamma }^2[b+5/2]}{b(1+b)(2+b)\pi }{\left(c_{s}v\right)}^{-2(1+b)}.\end{array}$$

(107)

The details on the integrations can be found in Appendix G. Thus, we have found the kernel for modes which are superhorizon before reheating. However, these induced tensor modes are not yet GWs so we must follow them after they enter the horizon in the radiation domination era.

Matching to Radiation Domination

The kernel (106) is valid for superhorizon modes during an arbitrary $b=\mathrm{constant}$ period. After this period ends with a sudden reheating, we shall assume that the source term is shut off and tensor modes evolve as freely propagating massless tensor mode. As explained in Section 4.3 the continuity of $h_{ij}$ implies continuity of the kernel. This means that the kernel goes from (106) to the kernel of a superhorizon tensor mode in radiation domination, which reads

$$\begin{array}{cc}\hfill I_{RD}(\overline{x}\ll 1,u,v)\approx & B_{1,RD}+B_{2,RD}{\left(k\overline{\tau }\right)}^{-1}\mathrm{with}\overline{\tau }\equiv \tau -\frac{b}{1+b}{\tau }_{\mathrm{rh}},\hfill \end{array}$$

(108)

where we introduced $\overline{\tau }$ from the continuity of the background FLRW metric, i.e., we imposed continuity of a and $\mathcal{H}$ at $\tau ={\tau }_{\mathrm{rh}}$. Matching (108) and its first derivatives with (106), we find that

$$\begin{array}{c}\hfill B_{1,RD}=B_{1,b}+\frac{1-b}{1+b}{B}_{2,b}{\left(k{\tau }_{\mathrm{rh}}\right)}^{-2b},B_{2,RD}=\frac{2b}{{(1+b)}^2}{B}_{2,b}{\left(k{\tau }_{\mathrm{rh}}\right)}^{1-2b}.\end{array}$$

(109)

Thus, we have effectively continued the superhorizon kernel from the $b=\mathrm{constant}$ era to the radiation era. We can now follow the kernel down to subhorizon scales, which must behave as a tensor mode in a radiation dominated universe as there is no source, that is Equation (75) with $b=0$. Then, taking the averaged square kernel yields

$$\begin{array}{cc}\hfill \overline{I_{RD}^2(k\ll k_{\mathrm{rh}},\tau \gg {\tau }_{\mathrm{rh}})}\approx & \frac{1}{2{\overline{x}}^2}\left(B_{1,RD}^2+B_{2,RD}^2\right).\hfill \end{array}$$

(110)

After picking up the most relevant contributions, that is the leading terms for $v\gg 1$, we arrive at

$$\begin{array}{cc}\hfill & \overline{I_{RD}^2(k\ll k_{\mathrm{rh}},\tau \gg {\tau }_{\mathrm{rh}})}\approx \frac{1}{2{\overline{x}}^2}{\left(\frac{{(3+2b)}^2}{4b{(1+b)}^2(2+b)}\right)}^2\hfill \\ & \times {\left((1+b+b^2){\left(c_{s}v\right)}^{-2}-4^{1+b}{\Gamma }^2[b+3/2]\frac{(1-b)}{\pi }{\left(c_{s}v\right)}^{-2(1+b)}{\left((1+b)\frac{k}{k_{\mathrm{rh}}}\right)}^{-2b}\right)}^2.\hfill \end{array}$$

(111)

This is the kernel squared that should be used in the induced GW spectrum (58) for scales $k\ll k_{\mathrm{rh}}$. In this way, the kernel (111) yields the right IR scaling for localized sources in a radiation dominated universe, i.e., ${\Omega }_{\mathrm{GWs}}^{\mathrm{IR}}\propto k^3$ [78]. The matching between the superhozion (111) and subhorizon (83) approximations at $k\sim k_{\mathrm{rh}}$ has shown to be very good in the instantaneous reheating case [74]. Note that for a gradual transition, the superhorizon approximation should still be good enough except for those modes which enter during the transition. We can then calculate the spectral density of induced GWs by

$${\Omega }_{\mathrm{GW}}(k\gg k_{\mathrm{rh}},\tau \gg {\tau }_{\mathrm{rh}})=\frac{k^2}{12a^2{H}^2}\overline{{\mathcal{P}}_h^{RD}},$$

(112)

with

$$\overline{{\mathcal{P}}_h^{RD}}=8{\int }_0^{\infty }dv{\int }_{|1-v|}^{1+v}du{\left(\frac{4v^2-{(1-u^2+v^2)}^2}{4uv}\right)}^2\overline{I_{RD}^2(x\gg 1,v,u)}{\mathcal{P}}_{\Phi }\left(ku\right){\mathcal{P}}_{\Phi }\left(kv\right),$$

(113)

and using Equation (111). Nevertheless, it is worth noting that for a Dirac delta or a fairly sharp peak and $b > 0$ a good approximation to the full induced GW spectrum is to stop the subhorizon spectral density (90) at $k\sim 3/4k_{\mathrm{rh}}$ and match with a power-law ${\Omega }_{\mathrm{GW}}={\Omega }_{\mathrm{GW},\mathrm{rh}}{(4k/3k_{\mathrm{rh}})}^2$ [74]. The amplitude of ${\Omega }_{\mathrm{GW},\mathrm{rh}}$ can be found by matching the two.

It should be noted that Equation (111) does not recover the logarithmic correction typical of induced GWs in a radiation dominated universe [79]. The main reason is that we assumed the source term to stop at reheating and, therefore, that superhorizon tensor modes at that time behave as free tensor modes. This is a good approximation taken for simplicity. However, the precise matching would require to first match the Newtonian potential $\Phi$ from a $b=\mathrm{constant}$ universe to radiation domination and then continue the kernel (60) after reheating with the matched $\Phi$. This would give a more accurate calculation and would most likely recover the logarithmic correction in the IR tail.

5. Typical Induced GW Spectra

In Section 4 we derived analytical approximations for the kernel (or transfer functions) of induced GWs. However, we still have the integral over the scalar momenta. While the kernel was not so sensitive to the shape of the primordial spectrum, i.e., the approximation is valid as long as all the relevant primordial spectrum enters the horizon before reheating, the integral over momenta crucially depends on the shape of the primordial spectrum. In this section, we will present the final form of the induced GW spectrum and discuss their main features. We will not be so concerned in the series of approximations to achieve analytical formulas of the GW spectrum, which can be found in the papers. Before doing so, we will carefully derive in Section 5.1 what is the spectral density of induced GWs measured today. Then we will focus on the simplest types of motivated primordial spectra. These are a very sharp peak Section 5.2, a log-normal peak Section 5.3, a broken power-law Section 5.4 and a scale invariant spectrum Section 5.4.2. We will also comment on possible oscillatory features Section 5.5 on top of the basic shapes and the impact of primordial NG on the final GW spectrum Section 5.6.

Main references: Although there are many studies on the induced GW spectrum for various primordial spectra, we focus on the cases where analytical formulas have been derived. In this way, Section 5.1 can be derived using the formulas in Appendix A and can also be found in [64,74]. Section 5.2 is mainly based on [74]. Section 5.3 follows [74,82]. Then, Section 5.4 is based on [81]. Section 5.5 briefly describes the results of Refs. [84,87]. Lastly, Section 5.6 reports the main findings from Refs. [89,90,92].

5.1. The GW Spectral Density Today

To find the power of GWs that we would observe we have to evaluate (7) today, namely

$${\Omega }_{\mathrm{GW},0}=\frac{1}{3M_{\mathrm{pl}}^2{H}_0^2}\frac{d{\rho }_{\mathrm{GW},0}}{dlnk}.$$

(114)

However, our computations in Section 4 were done up to and including the transition to a radiation dominated universe. Since deep inside the horizon we have that ${\rho }_{\mathrm{GW}}\propto {\rho }_r\propto a^{-4}$ then after reheating we have that ${\Omega }_{\mathrm{GW}}={\Omega }_{\mathrm{GW},\mathrm{rh}}=\mathrm{constant}$. To relate ${\Omega }_{\mathrm{GW},\mathrm{rh}}$ with ${\Omega }_{\mathrm{GW},0}$ we can make use of the fact that GWs behave as radiation and write

$$\begin{array}{cc}\hfill {\Omega }_{\mathrm{GW},0}{h}^2& ={\Omega }_{r,0}{h}^2\frac{1}{{\rho }_{r,0}}\frac{d{\rho }_{\mathrm{GW},0}}{dlnk}\hfill \\ & =1.62\times {10}^{-5}\left(\frac{{\Omega }_{r,0}{h}^2}{4.18\times {10}^{-5}}\right)\left(\frac{g_{*}\left(T_{\mathrm{rh}}\right)}{106.75}\right){\left(\frac{g_{*s}\left(T_{\mathrm{rh}}\right)}{106.75}\right)}^{-4/3}{\Omega }_{\mathrm{GW},\mathrm{rh}}.\hfill \end{array}$$

(115)

where ${\Omega }_{r,0}{h}^2\approx 4.18\times {10}^{-5}$ is the density fraction of radiation today given by Planck [2]. We also traced ${\rho }_{r,0}$ and ${\rho }_{\mathrm{GW},0}$ back to the reheating time taking into account that the effective degrees of freedom change with temperature. In the case that the induced GWs are generated during a radiation domination era we should replace the subscript “rh” for the evaluation at a time when the induced tensor modes start to behave as a GW. This is denoted with a subscript “c” by Inomata et al. [126].

5.2. Dirac Delta Peak

The simplest case is that of a very sharp peak in the primordial spectrum. While this type of spectrum is not possible to generate during single field inflation [103,109,234], it may be achieved in multi-field models of inflation [95,235,236,237,238,239,240,241,243]. If the peak is extremely sharp, it is often modelled by a Dirac delta as

$$\begin{array}{c}\hfill {\mathcal{P}}_{\mathcal{R}}\left(k\right)={\mathcal{A}}_{\mathcal{R}}\delta \left(ln(k/k_p)\right).\end{array}$$

(116)

Although this might be an unrealistic situation, it captures the essence of sharp peaks and gives insight on the kernel. We discuss a more realistic situation in Section 5.3. Now, for the power spectrum inside the integral (58) we have

$$\begin{array}{c}\hfill {\mathcal{P}}_{\mathcal{R}}\left(vk\right)={\mathcal{A}}_{\mathcal{R}}\delta \left(ln(vk/k_p)\right)={\mathcal{A}}_{\mathcal{R}}\frac{k_p}{k}\delta \left(v-\frac{k_p}{k}\right),\end{array}$$

(117)

and similarly for u. Recall, however, that the integral over u is bounded by momentum conservation to $|1-v| < u < 1+v$. So, when we integrate over v and u using the Dirac delta, the range is now $|k-k_p|<k_p < k+k_p$. This translates to a range in k given by $0 < k < 2k_p$. After taking this finite range of k into consideration, we may simply evaluate at $v=u=k_p/k$ the integrand (89) (and (113)) and insert it into (88) (and (112)) to arrive at

$$\begin{array}{c}\hfill {\Omega }_{\mathrm{GW},\mathrm{rh}}\left(k\right)=\frac{2}{3}{\left(\frac{k_p}{k_{\mathrm{rh}}}\right)}^2{\left(1-\frac{k^2}{4k_p^2}\right)}^2\overline{I_{RD}^2\left(k/k_{\mathrm{rh}},k/k_p\right)}\Theta (2k_p-k).\end{array}$$

(118)

Note that we are being general by allowing a $b=\mathrm{constant}$ phase before the standard radiation domination. Since we are assuming an instantaneous transition at $\tau ={\tau }_{\mathrm{rh}}$ from the $b=\mathrm{constant}$ to radiation domination, we have to include the scale $k_{\mathrm{rh}}$ corresponding to the last scale that entered the horizon at the transition. In more detail, $k_{\mathrm{rh}}$ is related to ${\tau }_{\mathrm{rh}}$ by

$$\begin{array}{c}\hfill k_{\mathrm{rh}}={\mathcal{H}}_{\mathrm{rh}}=\frac{1+b}{{\tau }_{\mathrm{rh}}}.\end{array}$$

(119)

Thus, we have that $x_{\mathrm{rh}}=k{\tau }_{\mathrm{rh}}=(1+b)k/k_{\mathrm{rh}}$. If we are considering the induced GWs in a radiation dominated universe we may simply replace $k_{\mathrm{rh}}\to k_p$ and recover the standard formula [63]. In Figure 3 and Figure 4 we see the GW spectrum for different values of b and respectively for $c_s^2=w$ and $c_s^2=1$, so that the resonance is evident.

Despite the fact that (118) is already a completely analytical formula, the functional complication does not provide much insight yet. Nevertheless, we can check one interesting limit, which is the IR tail of the spectrum. Expanding (118) for $k\ll k_p$ we arrive at a general analytical approximation for the low frequency tail of the GW spectrum. First, for $b < 0$ we have that

$$\begin{array}{cc}\hfill {\Omega }_{\mathrm{GW},\mathrm{rh}}(b<0,k\ll k_p)=& \frac{{\mathcal{A}}_{\mathcal{R}}^2}{12}{\left(\frac{2^{1+b}(2+b){\Gamma }^2[3/2+b]}{\pi c_s^{2(1+b)}{\left(1+b\right)}^{1+b}}\right)}^2{\left(\frac{k_{\mathrm{rh}}}{k_p}\right)}^{2+4b}\hfill \\ & \times \left\{\begin{array}{ccc}\hfill & \frac{2^{3+2b}}{\pi {\left(1+b\right)}^{2b}}{\left(\frac{k}{k_{\mathrm{rh}}}\right)}^2\hfill & \hfill (k\lesssim k_{\mathrm{rh}}/(1+b))\\ \hfill & {\left(\frac{\pi }{sin\left(b\pi \right)\Gamma [2+b]}\right)}^2{\left(\frac{k}{k_{\mathrm{rh}}}\right)}^{2+2b}\hfill & \hfill (k\gtrsim k_{\mathrm{rh}}/(1+b))\end{array}\right.,\hfill \end{array}$$

(120)

where we set the broken power-law knee at $k\sim k_{\mathrm{rh}}/(1+b)$ since it is when the subhorizon approximation breaks down for $b < 0$. A more detailed discussion can be found in Section 4.2. Second, for $b > 0$ we find

$$\begin{array}{cc}\hfill {\Omega }_{\mathrm{GW},\mathrm{rh}}(b > 0,k\ll k_p)=& \frac{{\mathcal{A}}_{\mathcal{R}}^2}{24\pi }{\left(\frac{(2+b)(1+b+b^2)}{c_s^{2}b{\left(1+b\right)}^2}\right)}^2{\left(\frac{k_{\mathrm{rh}}}{k_p}\right)}^2\hfill \\ & \times \left\{\begin{array}{ccc}\hfill & {\left(\frac{k}{k_{\mathrm{rh}}}\right)}^2\hfill & \hfill (k\lesssim k_{\mathrm{rh}})\\ \hfill & \frac{1}{2}{\left(\frac{2^{1+b}\Gamma [b+3/2]}{{\left(1+b\right)}^{1+b}}\right)}^2{\left(\frac{k}{k_{\mathrm{rh}}}\right)}^{2-2b}\hfill & \hfill (k\gtrsim k_{\mathrm{rh}})\end{array}\right.,\hfill \end{array}$$

(121)

where this time the matching is roughly at $k\sim k_{\mathrm{rh}}$. The above formulas agree with those in Ref. [74]. The difference between $b > 0$ and $b < 0$ is the continued superhorizon growth in the case of $b < 0$ which is explained in Section 2. The most important aspect of Equations (120) and (121) is that it clearly shows how the IR spectral tilt of the induced GW spectrum is directly related to the equation of state parameter at the time of induced GW generation. Thus, the observation of the IR tail gives insight on the expansion history of the primordial universe. Interestingly, for $b > 1$ ($w < 0$) the spectrum has a peak at $k\sim k_{\mathrm{rh}}$ which is larger than that at $k\sim k_p$. This means that induced GWs for $b < 1$ peaks at a different scale than the peak of the primordial spectrum. This has been used, e.g., in [132] to fit the NANOGrav results [18] and at the same time predict PBH with a mass smaller than the solar mass. If instead one assumes induced GWs in radiation domination to fit the NANOGrav frequency range, the associated PBH counterpart must be of the order of solar mass. We discuss more about these possibilities in Section 9. Another important point is that the amplitude of the spectrum at $k\sim k_p$ strongly depends on the ratio $k_{\mathrm{rh}}/k_p$ and often mildly on b. We thus see that the smaller the ratio $k_{\mathrm{rh}}/k_p$, the more suppressed is the spectrum at $k\sim k_{\mathrm{rh}}$. This is because the main induced GW generation occurs when the scalar mode $k_p$ enters the horizon. Then, since the modes with $k\sim k_{\mathrm{rh}}$ have been superhorizon until the transition they can only grow by a factor ${(k_{\mathrm{rh}}/k_p)}^2$. For $b < 0$ the argument is slightly different because of the second superhorizon growth and so the suppression is minor compared to $b > 0$ with a power index that depends on b.

5.3. Log-Normal Peak

We now turn to a more realistic situation where the peak in the primordial spectrum has a finite width. Inspired by several models of multi-field inflation the peak of the primordial spectrum may be parametrized by a log-normal [82], namely

$$\begin{array}{c}\hfill {\mathcal{P}}_{\mathcal{R}}\left(k\right)=\frac{{\mathcal{A}}_{\mathcal{R}}}{\sqrt{2\pi }\Delta }exp\left[-\frac{{ln}^2(k/k_p)}{2{\Delta }^2}\right],\end{array}$$

(122)

where $\Delta$ is the dimensionless width of the peak. Note that the power spectrum (122) is normalized, i.e., ${\int }_{-\infty }^{\infty }dlnk{\mathcal{P}}_{\mathcal{R}}\left(k\right)={\mathcal{A}}_{\mathcal{R}}$. The log-normal parametrization includes sharp peaks where $\Delta \ll 1$ [95,235,236,237,238,239,240,241,242,243] and broad peaks where $\Delta \gtrsim 1$ [62,85,102,128,238,244,245,246,247,248,249,250,251,252,253,254,255]. For $\Delta \to 0$ we recover the Dirac delta case.

In this review we are mainly interested in the case of a narrow peak ($\Delta \ll 1$). The reasons are its simplicity and that we can use our results for general b of Section 5.2. The interested reader for analytical approximation in the broad peak case ($\Delta \gtrsim 1$) in radiation domination is referred to the work of Pi and Sasaki [82]. Here we may directly borrow their result for $\Delta \ll 1$ which reads

$$\begin{array}{c}\hfill {\Omega }_{\mathrm{GW},\Delta }\left(k\right)=\mathrm{erf}\left(\frac{1}{\Delta }{sinh}^{-1}\frac{k}{2k_p}\right){\Omega }_{\mathrm{GW},\delta }\left(k\right),\end{array}$$

(123)

where $\mathrm{erf}\left(x\right)$ is the error function and ${\Omega }_{\mathrm{GW},\delta }\left(k\right)$ is the GW spectrum induced by a $\delta$-function peak given in Section 5.2. Then, Equation (123) provides the correction to the Dirac delta spectrum of Section 5.2 due to the small but finite width of the spectrum which leads to the $k^3$ IR scaling in radiation domination [78]. The most important effect then appears when the ratio of scales $k/k_p$ is smaller than the width $\Delta$. Indeed, when $k\ll k_p$ we have that

$$\begin{array}{c}\hfill {\Omega }_{\mathrm{GW},\Delta }\left(k\right)=\mathrm{erf}\left(\frac{k}{2k_p\Delta }\right){\Omega }_{\mathrm{GW},\delta }\left(k\right),\end{array}$$

(124)

which gives a correction of $\mathrm{erf}\left(k/\left(2k_p\Delta \right)\right)\sim k/\left(2k_p\Delta \right)$ when $\Delta \gg k/\left(2k_p\right)$. Thus, due to the presence of $\Delta$ our IR tail of Section 5.2 changes respectively for $2k_p\Delta < k_{\mathrm{rh}}$ and $2k_p\Delta > k_{\mathrm{rh}}$ to

$$\begin{array}{c}\hfill {\Omega }_{\mathrm{GWs}}(k\ll k_p,2k_p\Delta <k_{\mathrm{rh}})\propto {\mathcal{A}}_{\mathcal{R}}^2\left\{\begin{array}{cc}{\left(\frac{k}{k_p}\right)}^3\hfill & \hfill k\ll 2k_p\Delta \ll k_{\mathrm{rh}}\\ {\left(\frac{k}{k_p}\right)}^2\hfill & \hfill 2k_p\Delta \ll k\ll k_{\mathrm{rh}}\\ {\left(\frac{k}{k_p}\right)}^{2-2\left|b\right|}\hfill & \hfill k\ll 2k_p\Delta \ll k_p\end{array}\right.,\end{array}$$

(125)

or

$$\begin{array}{c}\hfill {\Omega }_{\mathrm{GWs}}(k\ll k_p,2k_p\Delta > k_{\mathrm{rh}})\propto {\mathcal{A}}_{\mathcal{R}}^2\left\{\begin{array}{cc}{\left(\frac{k}{k_p}\right)}^3\hfill & \hfill k\ll k_{\mathrm{rh}}\\ {\left(\frac{k}{k_p}\right)}^{3-2\left|b\right|}\hfill & \hfill k_{\mathrm{rh}}\ll k\ll 2k_p\Delta \\ {\left(\frac{k}{k_p}\right)}^{2-2\left|b\right|}\hfill & \hfill 2k_p\Delta \ll k\ll k_p\end{array}\right..\end{array}$$

(126)

Note that since $k_{\mathrm{rh}}/k_p\ll 1$ the case $2k_p\Delta > k_{\mathrm{rh}}$ still allows for a very sharp peak with $k_{\mathrm{rh}}/\left(2k_p\right) < \Delta \ll 1$. See how the inclusion of the finite width increases the richness of the induced GWs spectrum. We illustrate this in Figure 5. Additionally, note that if $b=0$ then the induced GW spectrum only has two different slopes and $k_{\mathrm{rh}}$ becomes meaningless.

Before ending this subsection, we briefly discuss the case when $\Delta > 1$. It is shown in Ref. [82] that in the induced GW spectrum in radiation domination for $\Delta > 1$ and near $k\sim k_p$ is well approximated by

$$\begin{array}{c}\hfill {\Omega }_{\mathrm{GWs}}^{\mathrm{peak}}(\Delta > 1)\approx 0.125\frac{{\mathcal{A}}_{\mathcal{R}}^2}{{\Delta }^2}{e}^{-\frac{{ln}^2(k/k_p)}{{\Delta }^2}}.\end{array}$$

(127)

This means that the GW spectrum induced by a log-normal peak in the primordial spectrum also follows a log-normal peak but with a shallower width given by $\Delta /\sqrt{2}$. The factor 2 is a reflection of the secondary nature of induced GWs [82].

5.4. Broken Power-Law

In Section 5.3 we have studied the log-normal sharp peak which is inspired in multi-field models of inflation. However, if the curvature perturbation is enhanced during single field inflation it often takes a broken power-law shape [103,109,234,256]. The induced GWs counterpart for a broken power-law primordial spectrum have been studied in Refs. [80,81,105,141,220]. Here we present the results of Ref. [81], which contain analytical approximations for the amplitude and spectral tilt of both the IR and UV tails of the induced GW spectrum. A detailed discussion on the approximations used can be found in [81]. At the end of this subsection, we also briefly discuss the case of a scale invariant primordial spectrum.

Let us parametrize the primordial curvature power spectrum as a strict broken power-law given by

$$\begin{array}{c}\hfill {\mathcal{P}}_{\mathcal{R}}={\mathcal{A}}_{\mathcal{R}}\left\{\begin{array}{cc}{\left(\frac{k}{k_p}\right)}^{n_{\mathrm{IR}}}\hfill & \hfill k\le k_p\\ {\left(\frac{k}{k_p}\right)}^{-n_{\mathrm{UV}}}\hfill & \hfill k\ge k_p\end{array}\right.,\end{array}$$

(128)

where $n_{\mathrm{IR}},n_{\mathrm{UV}} > 0$ respectively are the IR and UV spectral tilts of the spectrum which peaks at $k=k_p$. In single field models one often finds that $n_{\mathrm{IR}}\sim 4$ [103]. If the enhancement of the curvature fluctuations is due to a bump in the inflationary potential, then $n_{\mathrm{UV}}$ can be related to the second derivative of the potential at the bump [80,256]. In such models $n_{\mathrm{UV}}$ is also related to the local non-gaussianity parameter $f_{NL}$ [256]. In what follows, we will only focus on induced GWs in the radiation dominated universe and later discuss its generalisations to other values of the equation of state parameter. Currently, only analytical approximations for radiation domination are available. Now, using (91) and (128) we can compute the associated induced GW spectrum. Doing so for the IR tail, i.e., for modes with $k\ll k_p$, we find that

$$\begin{array}{c}\hfill {\Omega }_{\mathrm{GW},\mathrm{rh}}(k\ll k_p)\approx 12{\mathcal{A}}_{\mathcal{R}}^2\left(\frac{1}{2n_{\mathrm{IR}}-3}+\frac{1}{2n_{\mathrm{UV}}+3}\right){\left(\frac{k}{k_p}\right)}^3{ln}^2\left(\frac{k}{k_p}\right).\end{array}$$

(129)

Equation (129) is valid for $n_{\mathrm{IR}} > 3/2$ and $n_{\mathrm{UV}} > 0$ and gives a good approximation of the IR tail. The condition $n_{\mathrm{IR}} > 3/2$ comes from requiring convergence of the integral at large internal momenta ($u\sim v\gg 1$).

The UV tail of the induced GW spectrum presents two different behaviours. If $0 < n_{\mathrm{UV}} < 4$ the integral converges even in the strict limit $v\to 0$. This means that we can pull out of the integral (91) all the k dependence in ${\mathcal{P}}_{\mathcal{R}}$ and have that ${\Omega }_{\mathrm{GW},\mathrm{rh}}\sim {\mathcal{P}}_{\mathcal{R}}^2$. The remaining piece is just a numerical factor. In this way, we find that the induced GW spectrum in the UV reads

$$\begin{array}{c}\hfill {\Omega }_{\mathrm{GW},\mathrm{rh}}(k\gg k_p,n_{\mathrm{UV}}<4)\approx \frac{1}{12}{\mathcal{A}}_{\mathcal{R}}^{2}F\left(n_{\mathrm{UV}}\right){\left(\frac{k}{k_p}\right)}^{-2n_{\mathrm{UV}}},\end{array}$$

(130)

where $F\left(n_{\mathrm{UV}}\right)$, for $0 < n_{\mathrm{UV}} < 4$, is given by

$$\begin{array}{c}\hfill F\left(n_{\mathrm{UV}}\right)=8{\int }_0^{\infty }dv{\int }_{|1-v|}^{1+v}du{\left(\frac{4v^2-{(1-u^2+v^2)}^2}{4uv}\right)}^2{\left(uv\right)}^{n_{\mathrm{UV}}}\overline{I^2(n_{\mathrm{UV}},k,v,u)},\end{array}$$

(131)

and converges everywhere. The range of $n_{\mathrm{UV}}$ for which this approximation is valid depends on the value of b. We discuss this after presenting the results for radiation domination ($b=0$). Since (131) has to be computed numerically, we present here a numerical fit for $n_{\mathrm{UV}} < 4$, which reads

$$\begin{array}{c}\hfill F\left(n_{\mathrm{UV}}\right)\approx 41+\frac{16n_{\mathrm{UV}}^2}{\sqrt{16-n_{\mathrm{UV}}^2}}.\end{array}$$

(132)

The form of (132) has been found noting that the integral diverges in the exact limit $n_{\mathrm{UV}}\to 4$. After several trial and error attempts, a reasonable fit was found to be (132). A detailed study of (131) might yield better analytical insights than (132). For the moment Equation (132) gives a good order of magnitude estimate. Note that when $n_{\mathrm{UV}}=4$ the induced GW spectrum present a logarithmic divergence [80,81]. Lastly, for $n_{\mathrm{UV}} > 4$ the primordial spectrum has a sharp decay and, in analogy with the sharp peak, we do not expect that ${\Omega }_{\mathrm{GW},\mathrm{rh}}\sim {\mathcal{P}}_{\mathcal{R}}^2$. Instead, there should be a fast fall off at around $k\sim 2k_p$ by momentum conservation. Indeed when $n_{\mathrm{UV}} > 4$ we find that

$$\begin{array}{c}\hfill {\Omega }_{\mathrm{GW},\mathrm{rh}}(k\gg k_p,n_{\mathrm{UV}} > 4)\approx \frac{16}{3}{\mathcal{A}}_{\mathcal{R}}^2\left(\frac{1}{n_{\mathrm{UV}}-4}+\frac{1}{n_{\mathrm{IR}}+4}\right){\left(\frac{k}{k_p}\right)}^{-4-n_{\mathrm{UV}}}.\end{array}$$

(133)

We show the analytical estimates against the numerical results in Figure 6. See how even if the approximations are extrapolated to $k\sim k_p$ they still yield a sensible order of magnitude estimate of the induced GW spectrum.

With the above analytical estimates for the amplitude and the spectral tilt in the IR and UV limits we conclude that the induced GW spectrum has the following broken power-law shape

$$\begin{array}{c}\hfill {\Omega }_{\mathrm{GWs}}\left(k\right)\propto {\mathcal{A}}_{\mathcal{R}}^2\left\{\begin{array}{cc}{\left(\frac{k}{k_p}\right)}^3\hfill & \hfill k\ll k_p\\ {\left(\frac{k}{k_p}\right)}^{-\Delta }\hfill & \hfill k\gg k_p\end{array}\right.,\mathrm{with}\Delta =\left\{\begin{array}{ccc}\hfill & 2n_{\mathrm{UV}}\hfill & \hfill 0<n_{\mathrm{UV}}<4\\ \hfill & 4+n_{\mathrm{UV}}\hfill & \hfill n_{\mathrm{UV}} > 4\end{array}\right..\end{array}$$

(134)

Equation (134) tells us that the UV spectral tilt of the induced GW has information the UV tail of the primordial power spectrum. This is to be contrasted with the IR spectral tilt, which does not tell us much about the inflationary model unless we measure the amplitude very well. Lastly, if we focus on the UV tail of the induced GW spectrum (134), we see that for example it may be degenerate with the GW spectrum from the sound waves in a first order phase transition [172].

5.4.1. Alternative Expansion Histories

Here we briefly present the spectral tilts of the induced GW spectrum if instead of radiation the universe was dominated by another perfect fluid. In this case, we find that (134) changes to [81]

$$\begin{array}{c}\hfill {\Omega }_{\mathrm{GW}}\left(k\right)\propto {\mathcal{A}}_{\mathcal{R}}^2\left\{\begin{array}{cc}{\left(\frac{k}{k_p}\right)}^3\hfill & \hfill k\ll k_{\mathrm{rh}}\\ {\left(\frac{k}{k_p}\right)}^{3-2\left|b\right|}\hfill & \hfill k_{\mathrm{rh}}\ll k\ll k_p\\ {\left(\frac{k}{k_p}\right)}^{-\Delta -2b}\hfill & \hfill k\gg k_p\end{array}\right.,\end{array}$$

(135)

where b is defined in Equation (16) and this time

$$\begin{array}{c}\hfill \Delta =\left\{\begin{array}{cc}2n_{\mathrm{UV}}\hfill & \hfill 0<n_{\mathrm{UV}}+b<4\\ 4+n_{\mathrm{UV}}\hfill & \hfill n_{\mathrm{UV}}+b > 4\end{array}\right..\end{array}$$

(136)

Previously, in Section 2.2 and Section 4, we derived the IR broken power-law behaviour of the induced GW spectrum for general equation of state parameter. Now, in this section we have a power-law UV tail in the induced GW coming from the UV tail of the primordial spectrum. The additional factor $-2b$ with respect to (134) in the UV tail can be understood as follows. All modes with $k\gg k_p$ entered the horizon deep inside the $b=\mathrm{constant}$ era and before the scalar peak at $k\sim k_p$. After that the induced tensor modes effectively behave as free GWs irrespective of whether $b<0$ or $b > 0$. Thus, the free GWs experience the typical relative redshift with respect to the background, which is the factor $k^{-2b}$. Note that if $b<0$ the UV tail might have a blue tilt.

5.4.2. Scale Invariant Spectrum

A limiting case of the broken power-law studied (128), setting $n_{\mathrm{UV}}=n_{\mathrm{IR}}=0$, is a scale invariant spectrum typical of slow-roll inflation. In this case, it is easy to convince ourselves that the induced GW spectrum for general expansion histories is given by

$${\Omega }_{\mathrm{GW},\mathrm{rh}}={\mathcal{A}}_{\mathcal{R}}^2\mathcal{F}(b,c_s){\left(\frac{k}{k_{\mathrm{rh}}}\right)}^{-2b}.$$

(137)

where

$$\begin{array}{c}\hfill \mathcal{F}(b,c_s)=\frac{2}{3}{\int }_0^{\infty }dv{\int }_{|1-v|}^{1+v}du{\left(\frac{4v^2-{(1-u^2+v^2)}^2}{4uv}\right)}^2\overline{I^2(b,c_s,k,v,u)},\end{array}$$

(138)

and $\overline{I^2}$ is given by (83). Another explicit formula can be found at the summary Section 10. The integral (138) converges everywhere for any value of b and $c_s$. However, note that (138) with (83) is only valid for $k\gg k_{\mathrm{rh}}$. Nevertheless, since deep inside the radiation domination we expect a scale invariant induced GW spectrum it is reasonable to assume that for scales $k\sim k_{\mathrm{rh}}$ we can match the $k^{-2b}$ behaviour to the constant spectral density. A more detailed, perhaps numerical, analysis is needed to determine the spectrum for $k\sim k_{\mathrm{rh}}$ as induced GWs are being sourced before and after the transition.

5.5. Oscillatory Features

In the subsections above we have discussed the GWs induced by the simplest shapes of the primordial spectrum, that is a log-normal peak and a broken power-law. However, in some situations, it is expected that the primordial power spectrum presents a series of oscillations usually referred to as “primordial features”. For recent reviews on such features see, e.g., [257,258,259]. These oscillations can be classified according to the k dependence in their frequency. If the frequency is linear in k, the oscillations resulted from sharp dynamics during inflation. For example, a step in the potential in single field inflation [260,261,262,263,264,265,266,267] or a sharp turn in field space [84,86,254,255,268,269,270]. This case is also known as a sharp feature. If the frequency is logarithmic in k, the oscillations are called resonant features. They are due to resonances during inflation on subhorizon scales. For instance, resonances occur if the potential in single field inflation has wiggles or a massive scalar field is oscillating at the bottom of the potential [87,170,271,272,273,274,275,276,277]. In more general cases, such as general interactions [278,279] or alternative scenarios to inflation [280], the frequency of the oscillations in the primordial spectrum could be a power-law of k.

Remarkably, despite the induced GWs being a second order effect, oscillations are not completely washed out. The main difficulty is to derive analytical approximations to the induced GW spectrum. In this respect, Fumagalli, Renaux-Petel and Witkowski [84,87] and also Braglia, Chen and Hazra [86] did a detailed numerical analysis of the induced GW spectrum and derived motivated semi-analytical templates. Perhaps the most interesting case is when the oscillations in the primordial spectrum are $O\left(1\right)$. For example, this is the case when a sharp feature is responsible for the peak in the primordial spectrum [84,255]. Now, if the oscillations are large and separated linearly in k, we can gain intuition by looking at the sum of Dirac delta peaks [83,84]. In this case if we call ${\omega }_{\mathrm{lin}}$ to the frequency of the Dirac delta peaks in the primordial spectrum, then the frequency in the induced GW spectrum is ${\omega }_{\mathrm{lin}}^{\mathrm{GWs}}={\omega }_{\mathrm{lin}}/c_s$ [83,84] (in radiation domination $c_s=1/\sqrt{3}$). Interestingly, in the more realistic situation studied [84] the induced GW spectrum follows an envelope as if there were no oscillations and, on top of that, oscillations appear near the resonant peak at $k\sim 2c_s{k}_p$ and reach at most an amplitude of about $20\%$ in radiation domination [84]. We illustrate this possibility in Figure 7. Another example are $O\left(1\right)$ oscillations due to a resonant feature in the primordial spectrum. In this case, if we call ${\omega }_{\log }$ the frequency of the resonant feature, the resulting induced GW spectrum exhibits two oscillatory components with frequencies ${\omega }_{\log }$ and $2{\omega }_{\log }$. The relative amplitude of such components depend on both the frequency and the envelope of the primordial spectrum [87]. For the analytical templates and a more detailed discussion, see Refs. [84,87].

5.6. Impact of Non-Gaussianities

Here we briefly discuss the effects that primordial non-gaussianities may have on the induced GW spectrum. We will not dwell on the details, which are well explained in Refs. [89,90,92]. The leading formulas for local-type non-gaussianity can be found in Section 3.3. We describe the main effects of local-type non-gaussianity below.

First, let us focus on the very sharp peaks of Section 5.2 and Section 5.3. In this case, the most relevant effect is that the induced GW spectrum does not stop at $k\sim 2k_p$ but continues until $k\sim 3k_p$. This is roughly because when including primordial non-gaussianity, the source term to induced GWs has three scalars. Thus, we have a cut-off at $k\sim 3k_p$ by momentum conservation. Note that this effect is characteristic of the sharp peak in the primordial spectrum. For large non-gaussianity, that is $F_{NL}\gg 1$ in Equation (63), the spectrum does not have the clear resonant peak of Section 5.2 but presents a few peaks with similar amplitude instead [89,92]. Second, in the case where the peak is broad or a broken power-law, the non-gaussian contribution typically has the same IR and UV spectral tilts [81,92]. For large enough $F_{NL}\gtrsim 1$, the induced GW spectrum also shows an additional peak/bump at around $k\sim 3k_p$. If the non-gaussianity is large enough, i.e., ${\mathcal{A}}_{\mathcal{R}}{F}_{NL}^2 > O\left(1\right)$, the non-gaussian contribution starts to dominate over the gaussian one [89,92]. However, it is not so clear whether one can achieve ${\mathcal{A}}_{\mathcal{R}}{F}_{NL}^2 > 1$ in inflationary models as in this case the non-gaussian expression (63) cannot be regarded as a perturbative expansion [81]. A more detailed and general analysis is necessary. It should be noted that if we include third order perturbation theory terms, the induced GW spectrum for a Dirac delta peak has three resonant peaks instead of one [143,179].

The primordial non-gaussian contribution to the induced GW spectrum, as we have just presented it, might not yield much information on inflation. For example, one may design a gaussian primordial spectrum such that the resulting induced GW spectrum is exactly as that of a primordial log-normal peak with large non-gaussianities. Although this may be a very complicated task, it is in principle possible. In this sense, there is no telling primordial non-gaussianity apart just from the induced GW spectrum [179]. Nevertheless, the story becomes more interesting when considering the PBH counterpart. Since PBH formation is highly dependent on the tail of the distribution of primordial fluctuations, any tiny deviation from a normal distribution could have an exponential impact on the fraction of PBHs. For the details see the recent discussions in Refs. [217,219,220,221,256,281,282,283,284]. Then, ideally one may use the two observations, the PBH mass function and the induced GW counterpart, to perhaps disentangle the primordial non-gaussianity from the primordial spectrum.

6. The Dust Dominated Universe

In previous sections we have derived intuitive and analytical formulas for the spectrum of induced GWs under the assumption that the universe was not dominated by a pressure-less perfect fluid with vanishing propagation speed of fluctuations, i.e., $w=c_s^2=0$, so-called dust. The particularity of a dust dominated universe is that the Newtonian potential is constant on all scales. As can be seen from Equation (70) with $c_s^2=0$, $\Phi =\mathrm{constant}$ is a solution. This also implies that density fluctuations of the dust fluid grow. In fact, they grow as $\delta \rho /\rho \propto a$. This is the typical behaviour of a CDM dominated universe. This also means that there is a scale below which density fluctuations grow so much that they become non-linear, i.e., larger than $O\left(1\right)$. This occurs for modes with [51]

$$\begin{array}{c}\hfill k > k_{NL}\sim \sqrt{\frac{3}{2}}{\mathcal{P}}_{\Phi }^{-1/4}\left({\tau }_{\mathrm{rh}}\right)\mathcal{H}\left({\tau }_{\mathrm{rh}}\right).\end{array}$$

(139)

Note that it does not mean that if there are fluctuations on $k > k_{NL}$ there will be no production of induced GWs. It means that we cannot completely trust our calculations since non-linear dynamics might be important, e.g., due to GWs from mergers. Nevertheless, we quite often have $\Phi <1$ and thus perturbation theory for $\Phi$ is sensible.

In Section 6.1 we show and discuss the subtleties of the dust dominated universe and derive analytical estimates for the amplitude of induced GWs in the case of an instantaneous transition. We apply our results to PBH dominated epoch in Section 6.2 and use overproduction of induced GWs to place stringent constraints on the scenario.

Main references: Section 6.1 is based on Inomata et al. [65,66] with the notation of Refs. [70,71]. These works are built upon previous papers by Assadullahi and Wands [51] and Kohri and Terada [63]. Then, Section 6.2 follows the pioneering work by Papanikolau, Vennin and Langlois [69] further extended by [70] to account for the transition to radiation domination. Lastly, Ref. [68] contains useful details on the transition of the PBH dominated universe to radiation domination.

6.1. General Dust Domination

Let us now derive the kernels and some estimates for the induced GWs. First we note that the fact that $\Phi =\mathrm{constant}$ in an early Matter Dominated (eMD) era seemingly makes naive calculations of induced GWs easier. Just take a transfer function equal to unity in (53), which yields

$$\begin{array}{c}\hfill f(\tau ,q,|\mathbf{k}-\mathbf{q}\left|\right)=\frac{5}{3},\end{array}$$

(140)

where we took $b=1$ ($w=0$). Now we have to simply integrate (60) with a constant source. In the limit where $x\gg 1$, we have that the averaged kernel squared is just a constant [51,63]. In our notation we have

$$\begin{array}{c}\hfill \overline{I_{eMD}^2(x\gg 1,q,|\mathbf{k}-\mathbf{q}\left|\right)}\approx {\left(\frac{20}{3}\right)}^2.\end{array}$$

(141)

We can then use (58) to give a naive estimation of the induced GW during the early matter era. However, this is far from being a complete picture. As first realised by Inomata, Kohri, Nakama and Terada, induced GWs from a dust dominated are more subtle and depend on the nature of the transition to radiation domination [65,66]. If the transition is gradual, the induced GW spectrum is actually suppressed [65]. The suppression can be understood from the fact that the initial enhancement in the dust dominated universe found in Ref. [51] is due to the constant $\Phi$. If $\Phi$ decays, the enhancement is not as efficient. Instead, if the transition from dust to radiation is instantaneous, the induced GW spectrum gets amplified very much [66]. We will shortly explain the source of the amplification. Let us first note that the reason why we did not find such behaviour in Section 4 comes down to the particular behaviour of $\Phi$ in dust domination. In Section 4, due to the fact that we assumed $c_s^2\ne 0$, the general solution to $\Phi$ in a $w=\mathrm{constant}$ universe is a decaying and oscillating function. After the transition to radiation domination, $\Phi$ is also decaying and oscillating, albeit with different frequency. Thus, the matching between the two periods is well described within the WKB approximation. In this case, only those scales close to $k\sim k_{\mathrm{rh}}$ might receive at most a $O\left(1\right)$ correction.

What is special about the instantaneous reheating in the dust dominated universe is that we are suddenly going from ${\Phi }^{\prime }=0$ to ${\Phi }^{\prime }\ne 0$. Since $\Phi$ did not have time to decay and the frequency of the oscillations is proportional to the ratio of the largest wavenumbers k with the reheating scale $k_{\mathrm{rh}}$, the amplitude of ${\Phi }^{\prime }$ right after the sudden reheating will be huge. A more physical picture was presented in Ref. [66] and it goes as follows. Density fluctuations grow as $\delta \rho /\rho \propto a$ and attain a “large” amplitude by the end of the dust dominated era. These large fluctuations suddenly evaporate and are converted into fluctuations in a radiation fluid, which create strong pressure waves. These waves generate spacetime oscillations and are the source to induced GWs.

Let us explicitly show that ${\Phi }^{\prime }$ becomes very large right after reheating. To do that, we first consider that $\Phi$ has an arbitrary amplitude during the eMD of $\Phi ={\Phi }_{eMD}\left(k\right)$. Then, the eMD ends abruptly and we enter a late Radiation Domination (lRD). Matching the general solution (73) and its first derivative for $b=1$ and $b=0$ we arrive at

$$\begin{array}{c}\hfill {\Phi }_{lRD}={\left(c_{s}k\overline{\tau }\right)}^{-3/2}\left(C_1{J}_{3/2}\left(c_{s}k\overline{\tau }\right)+C_2{\gamma }_{3/2}\left(c_{s}k\overline{\tau }\right)\right),\end{array}$$

(142)

where $c_s=1/\sqrt{3}$ since we are in a radiation dominated universe, $\overline{\tau }$ is defined in (84) and comes from requiring continuity of a and $\mathcal{H}$ and we defined

$$\begin{array}{cc}\hfill C_1& =-{\Phi }_{\mathrm{eMD}}\frac{\pi }{2}{\left(\frac{c_{s}k{\overline{\tau }}_{\mathrm{rh}}}{2}\right)}^{5/2}{\gamma }_{5/2}\left(\frac{c_{s}k{\overline{\tau }}_{\mathrm{rh}}}{2}\right),\hfill \end{array}$$

(143)

$$\begin{array}{cc}\hfill C_2& ={\Phi }_{\mathrm{eMD}}\frac{\pi }{2}{\left(\frac{c_{s}k{\overline{\tau }}_{\mathrm{rh}}}{2}\right)}^{5/2}{J}_{5/2}\left(\frac{c_{s}k{\overline{\tau }}_{\mathrm{rh}}}{2}\right).\hfill \end{array}$$

(144)

If we now compute the variation of ${\Phi }_{\mathrm{lRD}}$ per Hubble time we find that the leading contribution for $x > x_{\mathrm{rh}}\gg 1$ is given by

$$\begin{array}{c}\hfill \frac{d{\Phi }_{\mathrm{lRD}}}{dln\overline{x}}\approx -{\Phi }_{\mathrm{eMD}}\frac{c_s{\overline{x}}_{\mathrm{rh}}^2}{\overline{x}}sin\left(c_s(\overline{x}-{\overline{x}}_{\mathrm{rh}})\right).\end{array}$$

(145)

Thus, we see that the amplitude of the time derivative is roughly

$$\begin{array}{c}\hfill \frac{d{\Phi }_{lRD}}{dln\overline{x}}\propto {\left(\frac{k}{k_{\mathrm{rh}}}\right)}^2.\end{array}$$

(146)

This means that if there are fluctuations on scales $k\gg k_{\mathrm{rh}}$ the amplitude of ${\Phi }_{\mathrm{lRD}}^{\prime }$ is very large and so it will produce induced GWs with a large amplitude.

Let us move on to the GWs induced after the sudden transition. For simplicity, we only focus on the largest contribution to the source (53), which comes from the time derivatives and reads

$$\begin{array}{cc}\hfill f(\tau ,q,|\mathbf{k}-\mathbf{q}\left|\right)\approx & \frac{1}{2}\frac{d{T}_{\Phi }\left(q\tau \right)}{dln\tau }\frac{d{T}_{\Phi }\left(\right|\mathbf{k}-\mathbf{q}\left|\tau \right)}{dln\tau }.\hfill \end{array}$$

(147)

Inserting Equations (145) and (147) into (60) we arrive at the kernel during the lRD, namely

$$\begin{array}{c}\hfill I_{lRD}(\overline{x},u,v)\approx {\Phi }_{\mathrm{eMD}}^2\frac{c_s^2{\overline{x}}_{\mathrm{rh}}^4}{2\overline{x}}uv{\int }_{{\overline{x}}_{\mathrm{rh}}}^{\overline{x}}d\tilde{\overline{x}}\frac{sin(\overline{x}-\tilde{\overline{x}})}{\tilde{\overline{x}}}sin\left(c_{s}v(\tilde{\overline{x}}-{\overline{x}}_{\mathrm{rh}})\right)sin\left(c_{s}u(\tilde{\overline{x}}-{\overline{x}}_{\mathrm{rh}})\right).\end{array}$$

(148)

The integral in (148) may be done analytically in terms of sine and cosine integrals but it is not so illuminating. The interested reader is referred to Ref. [66]. Using our experience of Section 4, we know that the integral often diverges when $c_s(u+v)\sim 1$. Thus, within our order of magnitude estimates, we shall focus solely on the neighbourhood of the resonance. In this case, the divergent piece is the cosine integral $\mathrm{Ci}\left[x\right]$ and the kernel is approximately given by

$$\begin{array}{c}\hfill I_{lRD,\mathrm{res}}(\overline{x},u,v)\approx {\Phi }_{\mathrm{eMD}}^2\frac{c_s^2{\overline{x}}_{\mathrm{rh}}^4}{8\overline{x}}uv\mathrm{Ci}\left[|1-c_s(u+v)|{\overline{x}}_{\mathrm{rh}}\right]sin\overline{x}.\end{array}$$

(149)

The calculation of the averaged kernel squared is straightforward and yields

$$\begin{array}{c}\hfill \overline{I_{lRD,\mathrm{res}}^2(\overline{x},u,v)}\approx {\Phi }_{\mathrm{eMD}}^4\frac{c_s^4{\overline{x}}_{\mathrm{rh}}^8}{128{\overline{x}}^2}{u}^2{v}^2{\mathrm{Ci}}^2\left[|1-c_s(u+v)|{\overline{x}}_{\mathrm{rh}}\right].\end{array}$$

(150)

For other approximations regarding the IR tail, such as the large momentum contribution $u\sim v\gg 1$, we refer the reader to Refs. [66,68].

With the analytic approximation for the kernel (150) we shall turn our attention to the dominant contribution to the tensor spectrum (58). For clarity, let us explicitly insert (150) into (58) to obtain the resonant part of the tensor spectrum, which leads us

$$\begin{array}{cc}\hfill \overline{{\mathcal{P}}_h^{\mathrm{res}}}\approx \frac{c_s^4{\overline{x}}_{\mathrm{rh}}^8}{16{\overline{x}}^2}{\int }_0^{\infty }dv{\int }_{|1-v|}^{1+v}du& {\left(\frac{4v^2-{(1-u^2+v^2)}^2}{4uv}\right)}^2{\mathcal{P}}_{\Phi }\left(ku\right){\mathcal{P}}_{\Phi }\left(kv\right)\hfill \\ & \times u^2{v}^2{\mathrm{Ci}}^2\left[|1-c_s(u+v)|{\overline{x}}_{\mathrm{rh}}\right].\hfill \end{array}$$

(151)

Before carrying out the integral, let us analyse some of the assumptions. First, we are assuming that (151) yields the dominant contribution. This will be the case if there is integration support for large k. This means that ${\mathcal{P}}_{\Phi }\left(k\right)$ should peak at large k. We also argued that there might be a limit to the magnitude of k due to non-linearities given by (139). We should use such a limit if there are no stronger motivations to consider $k > k_{NL}$. However, as we shall see, if PBH dominate the universe, there must be fluctuations above $k_{NL}$ [69]. For these reasons, let us parameterize the scalar power spectrum as power-law with a UV cut-off $k_{UV}$ given by [71]

$$\begin{array}{c}\hfill {\mathcal{P}}_{\Phi }={\mathcal{A}}_{\Phi }{\left(\frac{k}{k_{\mathrm{UV}}}\right)}^{-n}\Theta (k_{UV}-k),\end{array}$$

(152)

where ${\mathcal{A}}_{\Phi }$ is the amplitude and n the spectral tilt. Looking at (151) we see that the integrand peaks at large k if $n < 2$. This includes an almost scale invariant spectrum ($n\sim 0$) and the PBH density fluctuations ($n\sim 1$) [69,70].

With these assumptions we shall proceed to perform the integral in (151). We shall use the fact that $\mathrm{Ci}\left[z\right]$ is divergent for $z\sim 0$ to evaluate the integrand at the resonant point except for $\mathrm{Ci}\left[z\right]$. For simplicity, let us switch integration variables to

$$\begin{array}{c}\hfill z\equiv (1-c_s(u+v)){\overline{x}}_{\mathrm{rh}}\mathrm{and}s\equiv u-v,\end{array}$$

(153)

which have a Jacobian given by $\left|J\right|=1/\left(2c_s{\overline{x}}_{\mathrm{rh}}\right)$. Replacing $(u,v)$ for $(z,s)$ in (151), evaluating integrand at $z=0$ except for $\mathrm{Ci}\left[z\right]$ and maintaining the s dependence we arrive at

$$\begin{array}{c}\hfill \overline{{\mathcal{P}}_h^{\mathrm{res}}}\approx \frac{\pi {\overline{x}}_{\mathrm{rh}}^7}{512c_s{\overline{x}}^2}{\left(1-c_s^2\right)}^2{\int }_{-s_0\left(k\right)}^{s_0\left(k\right)}ds{(1-s^2)}^2{\mathcal{P}}_{\Phi }\left(k\frac{1+c_{s}s}{2c_s}\right){\mathcal{P}}_{\Phi }\left(k\frac{1-c_{s}s}{2c_s}\right),\end{array}$$

(154)

where we already integrated the cosine integral around the divergence19 and we have defined [66]

$$\begin{array}{c}\hfill s_0\left(k\right)=\left\{\begin{array}{ccc}\hfill & 1\hfill & \hfill k_{UV}/k\ge \frac{1+c_s^{-1}}{2}\\ \hfill & 2\frac{k_{UV}}{k}-c_s^{-1}\hfill & \hfill \frac{1+c_s^{-1}}{2}\ge \frac{k_{UV}}{k}\ge \frac{c_s^{-1}}{2}\\ \hfill & 0\hfill & \hfill \frac{c_s^{-1}}{2}\ge \frac{k_{UV}}{k}\end{array}\right..\end{array}$$

(155)

These bounds on s come from momentum conservation, i.e., $|1-v| < u < 1+v$, evaluated at the resonant point $z=0$. To have an order of magnitude estimate of the amplitude of the peak of the induced GW spectrum, we can evaluate (154) at $s=0$. By doing so, the amplitude of the induced GW spectrum at $k\sim k_{\mathrm{UV}}$ due to the sudden reheating is roughly given by

$$\begin{array}{c}\hfill {\Omega }_{\mathrm{GWs}}^{\mathrm{peak}}\approx \frac{\pi }{3\times 2^{12}{c}_s}{\left(1-c_s^2\right)}^2{\left(4c_s^2\right)}^n{\left(\frac{k_{\mathrm{UV}}}{k_{\mathrm{rh}}}\right)}^7{\mathcal{A}}_{\Phi }^2,\end{array}$$

(156)

where we used that ${\overline{x}}_{\mathrm{rh}}=k{\tau }_{\mathrm{rh}}/2=k/k_{\mathrm{rh}}$ and we used that the numerical evaluation of the integral in s roughly yields an extra factor $1/2$ [70]. Let us remind the reader that $c_s=1/\sqrt{3}$ since we are in radiation domination. From (156) we see how if $k_{\mathrm{UV}}\gg k_{\mathrm{rh}}$ the amplitude of the induced GWs is tremendously large by a power of 7. Thus, not to backreact, we either need ${\mathcal{A}}_{\Phi }$ to be extremely small or we restrict the value of $k_{\mathrm{UV}}$ not to be too large. Now, knowing the amplitude of the peak (156) we can approximate the induced GW spectrum by

$$\begin{array}{c}\hfill {\Omega }_{\mathrm{GWs},\mathrm{rh}}(k\sim k_{\mathrm{UV}})\approx {\Omega }_{\mathrm{GWs}}^{\mathrm{peak}}{\left(\frac{k}{k_{\mathrm{UV}}}\right)}^{7-2n}\Theta (k_{\mathrm{UV}}-k),\end{array}$$

(157)

where, although we used the cut-off in the induced GWs, which is at $k=2k_{\mathrm{UV}}$, the decay after $k\sim k_{\mathrm{UV}}$ is so fast that we can safely neglect it.

The resonant peak is not the only contribution after the sudden reheating induced GWs. There is also the IR tail which corresponds to the $u\sim v\gg 1$ region of integration. Nevertheless, this contribution is always suppressed by a factor $k_{\mathrm{rh}}/k_{\mathrm{UV}}$. The interested reader can find the corresponding formulas in Refs. [66,68]. In addition to that, there is the contribution of the GWs induced during eMD, i.e., using the kernel (141) into (58). In the case of the sudden transition, we can take the amplitude of the induced GWs as initial conditions for the subsequent radiation domination era. Then the total spectrum is well approximated by the sum of the two contributions [66]. This contribution will be small compared to the resonant part (157) but since it may present a large plateau [51,66,69] it may dominate on the very low frequency band. The case of the gradual transition requires a more careful treatment which can be found in Ref. [65]. We shall proceed with a very interesting application of (157).

6.2. PBH Dominated Era

An early matter dominated period could have been due to PBHs. Once formed, PBHs basically behave like a dust fluid, i.e., a fluid with no pressure and no propagation speed of fluctuations, like non-relativistic matter. Furthermore, their energy density practically redshifts as that of non-relativistic matter. You can see this from the following. Since the number density of PBH, $n_{\mathrm{PBH}}$, is conserved (unless they substantially merge or evaporate) we have that $n_{\mathrm{PBH}}\propto a^{-3}$. The total energy density contained in these PBHs is ${\rho }_{\mathrm{PBH}}=M_{\mathrm{PBH}}{n}_{\mathrm{PBH}}$. Then in periods where PBH evaporation is negligible we have that ${\rho }_{\mathrm{PBH}}\propto a^{-3}$. For concreteness, let us assume that PBH formed in a radiation dominated era, with an initial fraction of PBH given by

$$\begin{array}{c}\hfill \beta \equiv \frac{{\rho }_{\mathrm{PBH},i}}{{\rho }_{r,i}}\approx \frac{{\rho }_{\mathrm{PBH},i}}{3H_i^2{M}_{\mathrm{pl}}^2},\end{array}$$

(158)

where in the last step we used the first Friedmann Equation (A33). If $\beta$ is large enough, then PBH will eventually dominate the universe since we have that ${\rho }_r\propto a^{-4}$. Now, for simplicity, let us assume that these PBHs formed by the collapse of large primordial fluctuations with a very peaked primordial spectrum at a scale $k_p$. This has two implications. First, the PBH mass function is almost monochromatic. Then, we can use Equation (30) to estimate the peak mass in terms of the peak scale $k_p$ by using that $H_i\equiv k_p/a_i$. We refer the reader to Section 2.3 for more details. Second, the initial fraction of PBH $\beta$ is determined by the amplitude of the primordial spectrum. In this section, we take for convenience $M_{\mathrm{PBH}}$ and $\beta$ as the free parameters of the model. One can then relate the scale and the amplitude of the primordial spectrum to $M_{\mathrm{PBH}}$ and $\beta$, if necessary. Interestingly, in this situation the PBH evaporation occurs almost instantaneously [68] and, therefore, we can use the estimates we derived in Section 10.2 with some corrections that we describe below. It should be noted that the results we will derive here only apply to a sharply peaked PBH mass function. If we considered an extended PBH mass function, induced GWs would be very much suppressed [68].

A remarkable fact is that, as first realised by Papanikolaou, Vennin and Langlois [69], due to the inhomogeneous distribution, PBH themselves create density fluctuations which might later source induced GWs. To see this, note that PBH formation by the collapse of large fluctuations is a rare event [174] and so each PBH form uncorrelated of the others. Additionally, clustering is often negligible. This means that, as a good approximation PBHs will be randomly distributed uniformly in space. In other words, PBHs are distributed according to Poisson statistics. Such Poisson statistics are dictated by the mean inter-PBH co-moving separation. If PBH follow a monochromatic mass function, the mean inter-PBH co-moving separation at formation is given by

$$\begin{array}{c}\hfill d_i={\left(\frac{3}{4\pi n_{\mathrm{PBH},i}}\right)}^{1/3}\equiv k_{\mathrm{UV}}^{-1}.\end{array}$$

(159)

Note that in Equation (159) we identify the co-moving inter-PBH separation as the UV cut-off in Equation (152). This is because the fluid description of the PBH gas is valid only for $k < k_{\mathrm{UV}}$, i.e., at distances much larger than $d_i$. In this coarse grained regime the initial dimensionless power spectrum of PBH density fluctuations reads [69]

$$\begin{array}{c}\hfill {\mathcal{P}}_{\delta ,\mathrm{PBH}}=\frac{2}{3\pi }{\left(\frac{k}{k_{\mathrm{UV}}}\right)}^3\Theta (k_{\mathrm{UV}}-k).\end{array}$$

(160)

We have shown with Equation (160) that PBHs give rise to density fluctuations. However, in the current model, PBHs form in an early Radiation Dominated stage (eRD) and, therefore, such density fluctuations correspond to isocurvature fluctuations. Thus, they do not source induced GWs yet. Induced GWs are mainly generated by curvature fluctuations. If one is interested in the latter aspect, the general source term in a two fluid system can be found in Ref. [70]. The isocurvature nature of the PBH density fluctuations can be understood from the following [70]. PBH formation in the fluid picture may be regarded as a transition of a fraction of the homogeneous radiation fluid into dust. Yet, the total energy density of radiation remains homogeneous, while PBH are distributed randomly. Thus, the inhomogeneity due to PBHs must be isocurvature as the total energy density is homogeneous. Interestingly, such isocurvature is converted into curvature fluctuations when PBH dominate the universe. This type of transition from radiation to dust was studied analytically in detail by Kodama and Sasaki [285,286]. Thus, we can directly borrow their results which provide the transfer function for a vanishing initial curvature perturbation $\Phi$ and non-zero isocurvature deep inside the eMD era as

$$\begin{array}{c}\hfill T_{\mathrm{eRD}\to \mathrm{eMD}}^{\mathrm{iso}\to \mathrm{curv}}(k;a\gg a_{\mathrm{eq}})=\left\{\begin{array}{cc}\frac{1}{5}\hfill & \hfill k\ll k_{\mathrm{eq}}\\ \frac{3}{4}{\left(\frac{k_{\mathrm{eq}}}{k}\right)}^2\hfill & \hfill k\gg k_{\mathrm{eq}}\end{array}\right..\end{array}$$

(161)

The subscript “eq” refers to the early radiation-PBH equality. For an intuitive re-derivation, see Ref. [69]. Equation (161) tells us that the initial isocurvature has been transferred to a constant curvature perturbation. Modes which entered the horizon before the early radiation-PBH equality have decayed substantially and, hence, the suppression factor ${(k_{\mathrm{eq}}/k)}^2$. This means that the curvature power spectrum due to early PBH fluctuations in the eMD on the smallest scales, i.e., ($k_{\mathrm{UV}} > k\gg k_{\mathrm{eq}}$) is given by

$$\begin{array}{c}\hfill {\mathcal{P}}_{\Phi ,\mathrm{PBH}}=\frac{3}{8\pi }{\left(\frac{k_{\mathrm{eq}}}{k_{\mathrm{UV}}}\right)}^4\left(\frac{k_{\mathrm{UV}}}{k}\right)\Theta (k_{\mathrm{UV}}-k)\Theta (k-k_{\mathrm{eq}}).\end{array}$$

(162)

It is the curvature perturbation spectrum given by (162) which sources induced GWs during the eMD. However, we are most interested in the induced GWs generated right at the start of the lRD. As we shall see, since PBH evaporation has a finite duration even for the monochromatic case, we must take into account an important suppression. Then we shall use our estimate Equation (156). Before that, we need to understand and quantify the PBH domination a bit further.

Let us derive the conditions to have a PBH dominated era in the early universe allowed by current observations. We start with our PBHs given an initial mass $M_{\mathrm{PBH},i}$ (30) and an initial fraction $\beta$. After formation, these PBHs will evaporate emitting Hawking radiation. This means that $\beta$ cannot be too small, otherwise PBH evaporate before they dominate. To quantify this requirement it is enough to look at the evaporation time which for a non-spinning black hole reads [68,189,287]

$$\begin{array}{c}\hfill t_{\mathrm{eva}}\approx \frac{160}{3.8\pi g_H\left(T_{\mathrm{PBH}}\right)}\frac{M_{\mathrm{PBH},i}^3}{M_{\mathrm{pl}}^4},\end{array}$$

(163)

where $g_H\left(T_{\mathrm{PBH}}\right)$ are the spin-weighted degrees of freedom and

$$\begin{array}{c}\hfill T_{\mathrm{PBH}}\equiv M_{\mathrm{pl}}^2/M_{\mathrm{PBH},i}\approx 1.06\times {10}^9\mathrm{GeV}{\left(\frac{M_{\mathrm{PBH},\mathrm{f}}}{{10}^4\mathrm{g}}\right)}^{-1}.\end{array}$$

(164)

If we assume only the standard model of particle physics we have that for PBHs which evaporated before BBN then $g_H\left(T_{\mathrm{PBH}}\right)\approx 108$. If we include a possible PBH spin the evaporation time decreases, reaching a factor $1/2$ for the extremal BH case [288,289,290,291,292]. Thus, we can take into account the PBH spin by replacing $t_{\mathrm{eva}}\to t_{\mathrm{eva}}/2$. This effect is considered in the induced GWs by Ref. [71]. We will not consider PBH spin in what follows. From Equation (163) we see that the “reheating” temperature, i.e., the temperature of the radiation fluid which dominates the universe after PBH evaporate, is given solely in terms of $M_{\mathrm{PBH},i}$. Using the current cosmological parameters, which are given in Appendix A, we find that the evaporation temperature can be written as

$$\begin{array}{c}\hfill T_{\mathrm{eva}}\approx 2.76\times {10}^4\mathrm{GeV}{\left(\frac{M_{\mathrm{PBH},i}}{{10}^4\mathrm{g}}\right)}^{-3/2}{\left(\frac{g_H\left(T_{\mathrm{PBH}}\right)}{108}\right)}^{1/2}{\left(\frac{g_{*}\left(T_{\mathrm{eva}}\right)}{106.75}\right)}^{-1/4},\end{array}$$

(165)

where $g_{*}$ are the effective degrees of freedom in the energy density of radiation [293,294]. Here we find our first constraint on the model: PBH cannot evaporate too late if they are to reheat the universe. For a successful BBN the reheating temperature must be $T_{\mathrm{eva}} > 4\mathrm{MeV}$ [295,296,297,298]. This puts an upper bound on the PBH mass as

$$\begin{array}{c}\hfill M_{\mathrm{PBH},i}<5\times {10}^8\mathrm{g}.\end{array}$$

(166)

Using Equation (165), the co-moving scale corresponding to the horizon size at the time of evaporation is given by

$$\begin{array}{c}\hfill k_{\mathrm{eva}}\approx 4.7\times {10}^{11}{\mathrm{Mpc}}^{-1}{\left(\frac{g_H\left(T_{\mathrm{PBH}}\right)}{108}\right)}^{1/2}{\left(\frac{g_{*}\left(T_{\mathrm{eva}}\right)}{106.75}\right)}^{1/4}{\left(\frac{g_{*,s}\left(T_{\mathrm{eva}}\right)}{106.75}\right)}^{-1/3}{\left(\frac{M_{\mathrm{PBH},i}}{{10}^4\mathrm{g}}\right)}^{-3/2},\end{array}$$

(167)

where $g_{*,s}$ are the effective degrees of freedom in the entropy. With the above formulas we have that the end of the PBH domination era is fixed by the PBH mass. The second constraint on the model comes from requiring that PBHs dominate before they evaporate. This is a constrain on the initial fraction of PBHs in terms of the PBH mass, which is found to be [68,70]

$$\begin{array}{c}\hfill \beta > 6.35\times {10}^{-10}{\left(\frac{g_H\left(T_{\mathrm{PBH}}\right)}{108}\right)}^{1/2}{\left(\frac{\gamma }{0.2}\right)}^{-1/2}{\left(\frac{M_{\mathrm{PBH},f}}{{10}^4\mathrm{g}}\right)}^{-1}.\end{array}$$

(168)

Before we compute the final amplitude of the induced GW spectrum, it is very important to estimate the effects of the finite duration of the evaporation to the perturbations. Although the transition to a radiation dominated universe takes place within less than $1/4$ of e-fold, perturbations in very small scales are very much affected by the finite width of the transition, as was discovered by Inomata et al. [65,66,68]. The main reason is clear if we neglect the expansion of the universe, which is justified as we are interested in very subhorizon scales. As PBHs evaporate their energy density decays as [68]

$$\begin{array}{c}\hfill {\rho }_{\mathrm{PBH}}\propto M_{\mathrm{PBH}}\propto {\left(1-\frac{t}{t_{\mathrm{eva}}}\right)}^{1/3}.\end{array}$$

(169)

Now, we are interested in any effect that this decay may have on $\Phi$. The Newtonian potential is related to fluctuations by the Poisson equation which in Fourier space reads

$$\begin{array}{c}\hfill 2\frac{k^2}{a^2}\Phi \approx {\rho }_{\mathrm{PBH}}{\delta }_{\mathrm{PBH}}+{\rho }_r{\delta }_r,\end{array}$$

(170)

where ${\delta }_Q\equiv \delta {\rho }_Q/{\rho }_Q$ is the density contrast. We will do a big step here that can be easily confirmed by looking at the formulas in the appendix of Ref. [70] or in Ref. [68]. If we completely neglect the expansion of the universe we find that ${\delta }_{\mathrm{PBH}}\approx \mathrm{constant}$ is a solution to the equations of motion of PBH density fluctuations. With this solution, we can estimate the size of the fluctuations in the radiation fluid sourced by the evaporation. With the Green’s method to solve the equations of motion of radiation fluctuations, we find that

$$\begin{array}{c}\hfill {\delta }_r\sim \frac{{\rho }_{\mathrm{PBH}}}{{\rho }_r}\frac{a^2{\Gamma }^2}{k^2}{\delta }_{\mathrm{PBH}},\end{array}$$

(171)

where $\Gamma$ is the evaporation rate of PBHs and it is given by

$$\begin{array}{c}\hfill \Gamma \equiv -\frac{dln{M}_{\mathrm{PBH}}}{dt}\sim \frac{1}{t-t_{\mathrm{eva}}}.\end{array}$$

(172)

What (171) tell us is that fluctuations in the radiation fluid ${\delta }_r$ with $k/a\gg \Gamma$ are very suppressed with respect to ${\delta }_{\mathrm{PBH}}$ due to radiation pressure. This means that for such modes the dominant contribution to $\Phi$ is given by ${\delta }_{\mathrm{PBH}}$, even when PBH do not dominate the universe. Then by Equation (170) we see that $\Phi \propto {\rho }_{\mathrm{PBH}}$ until $\Gamma \sim k/a$, at which point the fluctuations in radiation dominate and perturbations behave as in a radiation dominated universe. The temporal decay of $\Phi$ due to the decay of ${\rho }_{\mathrm{PBH}}$ yields a relative suppression given by [68]

$$\begin{array}{c}\hfill S_{\Phi }(k/k_{\mathrm{eva}})\equiv \frac{\Phi }{{\Phi }_{\mathrm{instant}}}\approx {\left(\sqrt{\frac{2}{3}}\frac{k}{k_{\mathrm{eva}}}\right)}^{-1/3},\end{array}$$

(173)

where ${\Phi }_{\mathrm{instant}}$ refers to the value of $\Phi$ if the transition is instantaneous.

Now, we can use our estimate for an instantaneous transition (156) but taking into account the suppression of $\Phi$ until the transition to radiation domination is completely achieved. This yields a spectrum of curvature fluctuations at the start of the lRD given by

$$\begin{array}{cc}\hfill {\mathcal{P}}_{\Phi ,\mathrm{lRD}}& =S_{\Phi }^2(k/k_{\mathrm{eva}}){\mathcal{P}}_{\Phi ,\mathrm{PBH}}\hfill \\ & =\frac{3}{8\pi }{\left(\frac{2}{3}\right)}^{-1/3}{\left(\frac{k_{\mathrm{UV}}}{k_{\mathrm{eva}}}\right)}^{-2/3}{\left(\frac{k_{\mathrm{eq}}}{k_{\mathrm{UV}}}\right)}^4{\left(\frac{k}{k_{\mathrm{UV}}}\right)}^{-5/3}\Theta (k_{\mathrm{UV}}-k)\Theta (k-k_{\mathrm{eq}}).\hfill \end{array}$$

(174)

Comparing (174) with (152) we see that $n=5/3$ and

$$\begin{array}{c}\hfill {\mathcal{A}}_{\Phi }=\frac{3}{8\pi }{\left(\frac{2}{3}\right)}^{-1/3}{\left(\frac{k_{\mathrm{UV}}}{k_{\mathrm{eva}}}\right)}^{-2/3}{\left(\frac{k_{\mathrm{eq}}}{k_{\mathrm{UV}}}\right)}^4.\end{array}$$

(175)

The final ingredient to derive the induced GW spectrum is the hierarchy of the scales involved. First, note that the ratio between the cut-off $k_{\mathrm{UV}}$ and $k_{\mathrm{eva}}$ is independent of $\beta$ and it is given by [70]

$$\begin{array}{c}\hfill \frac{k_{\mathrm{UV}}}{k_{\mathrm{eva}}}\approx 2.3\times {10}^6{\left(\frac{g_H\left(T_{\mathrm{PBH}}\right)}{108}\right)}^{-1/3}{\left(\frac{M_{\mathrm{PBH},i}}{{10}^4\mathrm{g}}\right)}^{2/3}.\end{array}$$

(176)

The $\beta$ independence of the ratio is due to the fact that $H_{\mathrm{eva}}$ only depends on the PBH mass and $a_{\mathrm{eva}}/a_i\propto n_{\mathrm{PBH},i}^{1/3}$ because it is a dust dominated universe. The hierarchy is closed with the second and third relations respectively given by [70]

$$\begin{array}{c}\hfill \frac{k_{\mathrm{eq}}}{k_{\mathrm{UV}}}=\sqrt{2}{\gamma }^{1/3}{\beta }^{2/3}\mathrm{and}\frac{k_{\mathrm{eq}}}{k_p}=\sqrt{2}\beta .\end{array}$$

(177)

Now, inserting Equations (174), (176) and (177) into (157) yields an amplitude of induced GWs at the peak of approximately

$$\begin{array}{c}\hfill {\Omega }_{\mathrm{GWs},\mathrm{eva}}^{\mathrm{peak}}\approx 1.5\times {10}^{-2}{\left(\frac{\beta }{{10}^{-6}}\right)}^{16/3}{\left(\frac{\gamma }{0.2}\right)}^{8/3}{\left(\frac{g_H\left(T_{\mathrm{PBH}}\right)}{108}\right)}^{-17/9}{\left(\frac{M_{\mathrm{PBH},f}}{{10}^4\mathrm{g}}\right)}^{34/9}.\end{array}$$

(178)

Note that the estimate Equation (178) is the amplitude of induced GWs right after evaporation. It should also be noted that (178) is valid for very sharp PBH mass functions. As the mass function widens, the induced GW counterpart gets suppressed [68]. On top of that, we have that scales $k > k_{NL}$ (139), specially for $k\sim k_{UV}$, PBH density fluctuations ${\delta }_{\mathrm{PBH}}$ become $O\left(1\right)$ or larger. However, our estimate (178) has been derived in the linear regime and might receive corrections due to non-linear effects. Nevertheless, we expect it to be a crude order of magnitude estimate since the main source of induced GWs, the curvature perturbation, remains always smaller than unity,20 i.e., $\Phi \ll 1$. Numerical relativity simulations are needed to derive a preciser estimate of the induced GWs generated in the non-linear regime. We show the resulting induced GW spectrum in Figure 8. Using (167) and (176) we find that the peak of induced GWs today lies at a frequency given by

$$\begin{array}{c}\hfill f_{\mathrm{UV}}\approx 1.7\times {10}^3\mathrm{Hz}{\left(\frac{g_H\left(T_{\mathrm{PBH}}\right)}{108}\right)}^{1/6}{\left(\frac{g_{*}\left(T_{\mathrm{eva}}\right)}{106.75}\right)}^{1/4}{\left(\frac{g_{*,s}\left(T_{\mathrm{eva}}\right)}{106.75}\right)}^{-1/3}{\left(\frac{M_{\mathrm{PBH},i}}{{10}^4\mathrm{g}}\right)}^{-5/6}.\end{array}$$

(179)

Quite surprisingly, the frequency of the peak enters the observational range of LISA, DECIGO and LIGO for $5\times {10}^8\mathrm{g} > M_{\mathrm{PBH},\mathrm{i}} > 2\times {10}^4\mathrm{g}$. Thus, this scenario is testable in the future.

To compute the amount of induced GWs measured today we must consider the redshift of GWs until today. We did this in Section 5 with Equation (115). However, if we want to evaluate the amplitude of induced GWs at BBN, we only need to know the fraction of GWs that compose the total radiation after reheating. This can be done by considering the degrees of freedom of the standard model, which leads us to

$$\begin{array}{cc}\hfill {\Omega }_{\mathrm{GW},\mathrm{BBN}}& =0.39\left(\frac{g_{*}\left(T_{\mathrm{eva}}\right)}{106.75}\right){\left(\frac{g_{*s}\left(T_{\mathrm{eva}}\right)}{106.75}\right)}^{-4/3}{\Omega }_{\mathrm{GW},\mathrm{eva}}.\hfill \end{array}$$

(180)

Since these induced GWs act as additional radiation, they contribute to the effective number of relativistic species $N_{\mathrm{eff}}$ at BBN. The contribution of GWs to $N_{\mathrm{eff}}$ is by convention21 written as [29]

$$\begin{array}{c}\hfill {\Omega }_{\mathrm{GW},\mathrm{BBN}}=\frac{7}{8}{\left(\frac{4}{11}\right)}^{4/3}\Delta N_{\mathrm{eff}}<0.05,\end{array}$$

(181)

where in the last step we used that BBN [299,300] sets an upper bound $\Delta N_{\mathrm{eff}} < 0.2$.22 In this way, we can constrain the amplitude of GWs from existing BBN bounds. Note that to be precise the BBN bound is a constrain on the total energy density of relativistic particles. Thus, we should have actually integrated our induced GW spectrum over $lnk$ and then compare it with the bound (181). However, since the GW spectrum is very peaked and we are already working with order of magnitude estimates, it is justified to just look at the peak of the GW spectrum. If we use the current constraints from BBN, then we obtain that

$$\begin{array}{c}\hfill \beta <1.5\times {10}^{-6}{\left(\frac{M_{\mathrm{PBH},i}}{{10}^4\mathrm{g}}\right)}^{-17/24}.\end{array}$$

(182)

This is a new constraint on early epochs of PBH dominations which cannot be obtained by any other means.

Interestingly, PBH evaporation also emits “gravitons”, i.e., very high frequency GWs. While these gravitons are not observable by GW detectors, they might be seen as effective relativistic particles at BBN [288,289,290,291,292]. Future CMB probes such as CMB-S4 [303] will be able to improve current bounds down to $\Delta N_{\mathrm{eff}}\sim 0.02$. Unfortunately, the graviton contribution of non-spinning PBHs to $N_{\mathrm{eff}}$ is out of reach for future CMB probes. Nevertheless, if PBH have a large spin (and dominate the universe at some point), the graviton contribution to $N_{\mathrm{eff}}$ will be accessible to future experiments [291,292]. This signature from graviton emission together with induced GWs will be an interesting probe of PBH dominated eras. While the graviton contribution to $N_{\mathrm{eff}}$ tells us about PBH spin, the induced GW counterpart has information on the formation mechanism [71]. As a clarification, note that the graviton emission and induced GWs are two very different effects. While the former is emitted by Hawking evaporation, the latter is generated by the time dependence of density fluctuations. The composition of such density fluctuations is not important, e.g., if a fraction are gravitons, as long as it behaves as fluctuations of a radiation fluid.

Before ending this section, we should note that we also have the induced GWs generated in the PBH dominated stage by using the Kernel (141). This is studied in detail in Ref. [69]. However, as we discussed at the beginning of this section this contribution is sensitive to the detailed transition to radiation domination. In the present case, where the transition is almost instantaneous, we can continue the tensor mode amplitude build up during eMD to the lRD. Nevertheless, this contribution is always subdominant compared to the one after reheating [66]. Thus, our estimate (178) provides the dominant peak of the induced GW spectrum.

7. The Gauge Issue

Despite the large amount of research on induced GWs, there is still the open question of what we really observe. This was made explicit in 2017, when Hwang, Jeong and Noh [149] showed that the induced GW spectrum generated in a dust dominated universe depends very much (by orders of magnitude) on the gauge chosen. Although the gauge issue was mentioned earlier by Matarrese, Mollerach and Bruni [37], it was not until [149] that it received more attention. From then and on, there has been a big activity to understand the problem [150,151,152,153,154,155,156,157,158,159,160,161,162]. For example in Refs. [151,152] other cosmological backgrounds rather than dust were studied. They found that the gauge dependence is in general present. Later studies [153,154,155,156,157,158,159,160,161,162], proposed partial solutions and investigated in which situations the spectrum of induced GW might be well defined. We dedicate Section 7.1 to explain the origin of the gauge issue. Then, we discuss these approaches in detail in Section 7.2.

Main references: The discussion in Section 7.1 is mainly based on [150,161]. Then, in Section 7.2 we revisit the main arguments of [153,154] and explain the improved results of [161].

7.1. The Origin of the Issue

To understand the source of problem, we must distinguish between what we call tensor modes and GWs. On the theoretical side, we define tensor modes in a homogeneous and isotropic universe as the transverse-traceless fluctuations of the spatial metric.23 However, this definition must depend on how one chooses the spatial hypersurface. Although tensor modes are gauge invariant at linear order, they fail to be so at second order in cosmological perturbation theory. This is reflected by the fact that tensor, vector and scalar modes mix at second order [37]. This may pose a problem when we want to relate the tensor modes in a given slicing to the observational effect, which are the GWs. Note that this is not necessarily an “issue”. If the quantity related to observation is clear and calculations are done properly (perhaps in a reasonable coordinate system), one should get proper results. For example, in studies of the CMB the temperature anisotropies and the polarization patterns of the CMB are well defined at second order. The interested reader is referred to Refs. [306,307,308] and references therein.

The trouble with GW observations is that we only know the GW detector response to linear GWs. Only at first order the relation between the tensor modes and GWs in the linearised theory is clear, since tensor modes are gauge invariant at that order. There is a subtlety, though. A tensor mode whose wavelength is larger than the horizon is frozen, does not evolve. So it acts as a constant anisotropy or shear and, as such, it has nothing to do with actual GWs. Furthermore, since the wavelength of the tensor modes is larger than the curvature scale, the averaging procedure used to derive the (pseudo) energy momentum tensor of GWs in Section 2.2 does not apply. For these reasons, we will be only concerned about tensor modes on subhorizon scales. These tensor modes have a wavelength smaller than the horizon and behave as a wave, which are the GWs we could observe.

In the absence of the GW detector response at second order, we face the problem that we do not know the relation between tensor modes and the true GWs. In most situations though this should not really be a problem, unless one is choosing a very strange coordinate system to describe the observations. For instance, to describe the GWs measured by a GW detector today we usually do not need to consider second order terms as the amplitude of the GWs is already tiny enough. However, this issue is particularly relevant for induced GWs due to their secondary nature. Now, there are two important aspects of the problem. First, we have that the derived tensor spectrum depends on the gauge choice. Second, we expect that once the source term in (48) is absent, i.e., its amplitude decayed enough, such previously induced tensor modes behave as a freely propagating GWs. From that moment on we should be able to relate the induced tensor modes to GWs. Thus, we must reconcile our physical intuition with the fact that tensor modes are gauge dependent. To arrive to that, there are two points that we would like to make as clear as possible. The first one concerns the tensor modes $h_{ij}$ themselves and the second one to the actual observable, which in Section 2.2 we decided it would be the energy density of induced GWs (5). We present them in the following subsections and we later discuss the current state of the issue in Section 7.2.

7.1.1. The Definition of Tensor Modes Is Gauge Dependent

When working at second order in cosmological perturbation theory, it is very important to be clear about our convention. Here we will use the conformal decomposition from Equations (34) and (38) which in a general gauge reads

$$\begin{array}{cc}\hfill d{s}^2& =a^2\left(\tau \right)\left(-N^{2}d{\tau }^2+{\mathrm{e}}^{2\varphi }{{\rm Y}}_{ij}\left(d{x}^i+N^{i}d\tau \right)\left(d{x}^j+N^{j}d\tau \right)\right),\hfill \end{array}$$

(183)

where N and $N^i$ respectively are the lapse and shift vector. We then perturb the variables as follows

$$\begin{array}{cc}\hfill N& =1+\alpha ,N_i={\partial }_i\beta \hfill \end{array}$$

(184)

$$\begin{array}{cc}\hfill {\left[ln{\rm Y}\right]}_{ij}& =h_{ij}+2\left({\partial }_i{\partial }_j-\frac{1}{3}{\delta }_{ij}\Delta \right)E,\hfill \end{array}$$

(185)

where we only focused on the tensor and scalar components and $\varphi$ is already a perturbation. Our convention differs from the common expansions, which expands the metric linearly in the perturbation variables instead of exponentially. Nevertheless, our choice considerably simplifies calculations. To relate the two different conventions we refer the reader to the appendix of Ref. [161]. We also use the exponential mapping for a gauge transformation (e.g., see the review by Malik and Wands [164]) which naturally follows from the Hamiltonian flow (see e.g., [150] for the application to second order cosmological perturbations). To be more concrete, let us do an infinitesimal coordinate transformation by

$$\begin{array}{c}\hfill {\tilde{x}}^{\mu }=x^{\mu }-{\xi }^{\mu }\mathrm{with}{\xi }^{\mu }=(T,{\partial }^{i}L),\end{array}$$

(186)

where we focused only on the scalar component of the spatial vector ${\xi }^i$. Then a quantity Q in one gauge is related to the same quantity evaluated in another gauge $\tilde{Q}$ by the exponential mapping

$$\begin{array}{c}\hfill \tilde{Q}=e^{{\mathcal{E}}_{\xi }}Q,\end{array}$$

(187)

where ${\mathcal{E}}_{\xi }$ is the Lie derivative along $\xi$. With this convention, we find that the tensor modes at second order change according to

$$\begin{array}{cc}\hfill {\tilde{h}}_{ij}=h_{ij}-{\widehat{TT}}_{ij}{}^{ab}\{& -2{\partial }_a\left(\beta -E^{\prime }\right){\partial }_{b}T+2{\partial }_a{\partial }_{k}L{\partial }_k{\partial }_{b}E+{\partial }_{a}T{\partial }_{b}T+{\partial }_a{\partial }_{k}L{\partial }_b{\partial }_{k}L\}.\hfill \end{array}$$

(188)

Equation (188) is the first cause of the gauge issue for induced GWs. Since the source to induced tensor modes acts at second order in perturbation theory, the actual gauge in which we do the calculations might drastically change the form of the solutions for $h_{ij}$ [149]. The reader interested in the gauge transformations particularly related to induced GWs might go to Ref. [161]. General reviews of cosmological perturbation theory are listed at the Section 1.2.

7.1.2. The Definition of GW Energy Density Is Gauge Dependent

Once we showed that the tensor modes $h_{ij}$ are gauge dependent, let us look at the definition of the “observable” energy density of GWs. Using (2) we have that

$${\rho }_{\mathrm{GW}}=t_{00}^{\mathrm{GW}}=\frac{M_{pl}^2}{8a^2}〈{h\prime }^{ij}{h}_{ij}^{\prime }+{\partial }_k{h}^{ij}{\partial }^k{h}_{ij}〉.$$

(189)

Now, looking at (188) and (189) the first question that comes to mind is: $h_{ij}$ in (189) is evaluated in which gauge? Indeed, if we naively apply the gauge transformation (188) to (189) we see that

$${\tilde{\rho }}_{\mathrm{GW}}\sim \langle {\tilde{h}}^{\prime ij}{\tilde{h}}_{ij}^{\prime }\rangle ={\rho }_{\mathrm{GW}}+\langle {\partial }_{i}T{\partial }_{j}T{\partial }^{i}T{\partial }^{j}T\rangle +...,$$

(190)

where we neglected other terms to illustrate the main point, which is that the term quartic in T in (190) does not vanish in general. This shows something that we were expecting in Section 2: the energy density of GWs (189) in analogy of the Ricci flat result (2) is strictly speaking not well-defined on a cosmological background at second order in perturbation theory. Furthermore, it does not matter when you evaluate the energy density, e.g., today. In principle, we can find a coordinate transformation in which ${\rho }_{\mathrm{GW}}$ is very different from the usual estimate. For example, if we go to a frame where the detector is wildly oscillating, Equation (189) might not yield sensible results. This is the second and the most important cause for the gauge dependence of induced GWs.

7.2. Current “Solutions”

In view of the current lack of a description at second order in cosmological perturbation theory of the detector response to passing GWs, all we can do is to try to argue when Equation (189) yields sensible results. Before we start discussing the details, it is important to realize that a gauge invariant formulation of the equations of motion of induced GWs (48) does not help at all. There are infinitely many ways to define gauge invariant variables, each definition corresponding to a particular gauge choice. Thus, the same question remains, although with a slightly different language: which gauge invariant definition of the tensor modes goes into (189)? Similarly, if one tries to find a gauge invariant formulation of (189), which implies going to fourth order in cosmological perturbation theory (and it is clearly not the point), it is not clear what is the relation with the real observable. So, what can we do?

An interesting direction, first explored by [153] and later corrected by [154], is to argue what is the gauge (or coordinate frame) which best describes the GW detection. According to Ref. [153], based on the analogy with asymptotically flat spacetimes, the best gauge is the so-called transverse-traceless (TT) gauge where the metric reads

$$\begin{array}{c}\hfill d{s}^2=-d{t}^2+({\delta }_{ij}+h_{ij})d{x}^{i}d{x}^j.\end{array}$$

(191)

In Minkowski spacetime, there are only two degrees of freedom, the ones in $h_{ij}$, since there is no source to the energy momentum tensor. In this gauge and in the linearised theory, the GW detector is at rest, even for passing GWs, which makes the physical interpretation (and calculations) much easier [184,185]. The closest gauge to the TT gauge in a cosmological background is the synchronous gauge, which is a coordinate system spanned by a family of geodesics. In our convention the metric in the synchronous gauge is given by

$$\begin{array}{cc}\hfill d{s}^2& =a^2\left(\tau \right)\left(-d{\tau }^2+{\mathrm{e}}^{2\varphi }{{\rm Y}}_{ij}d{x}^{i}d{x}^j\right),\hfill \end{array}$$

(192)

where we set $\alpha =\beta =0$ in (183). Let us remind the reader that

$$\begin{array}{cc}\hfill {\left[ln{\rm Y}\right]}_{ij}& =h_{ij}+2\left({\partial }_i{\partial }_j-\frac{1}{3}{\delta }_{ij}\Delta \right)E,\hfill \end{array}$$

(193)

and therefore we have two more variables, $\varphi$ and E, than in (191). In the end, only one of them is dynamical. Thus, although at linear order (193) resembles the TT gauge (191) it is certainly not the same. At linear order and if we neglect the drift due to the expansion a, the metric (193) does describe a fixed detector. However, at second order the detector will not be fixed and will feel the influence of E. Another “problem” of the synchronous gauge is that it is not a completely fixed gauge. There is a residual gauge ambiguity which corresponds to a change of the reference geodesics. This gauge ambiguity turns out to be very important in the final spectrum of induced tensor modes. If one does not properly fix the synchronous gauge, then induced tensor modes could be very different than, e.g., in the Newtonian gauge [158]. This is because if we start with a family of geodesics defined in the very early universe and we use the same geodesics to describe the GW detection, it is very possible that these geodesics do not match the rest frame of the detector. The resulting induced tensor spectrum might be singular if those geodesics happen to have focusing singularities near the detector [161]. Nevertheless, by properly fixing the gauge one finds that the induced GW spectrum in the synchronous gauge and in the Newtonian gauge yield the same prediction in a radiation dominated universe [153,154,155,158]. This also holds for any cosmological background with a constant equation of state [161] if it is not dust, i.e., $c_s^2\ne 0$. This result provides evidence that perhaps the energy density of GWs (189) computed in the Newtonian gauge, which we do for simplicity of the calculations, yields meaningful results. However, it does not explain why these two predictions coincide. Furthermore, the original issue raised in [149] in a dust dominated universe persist.

Here is where the direction explored in Ref. [161] becomes relevant: when and why does the energy density of GWs (189) yield the same predictions and for which gauges? The physical intuition is that the analogy with Ricci flat spacetimes should work well if:

(i) 

we look at really subhorizon scales ($k\gg \mathcal{H}$), where cosmology should be less relevant;

(ii) 

there is (barely) no source of induced GWs (48), so that induced tensor modes are freely propagating GWs;

(iii) 

we use a coordinate system suitable for small distances calculations, so that we are not confused by strange coordinate artefacts.

As shown in [161], if points $\left(i\right)$–$\left(iii\right)$ are satisfied then (189) is gauge independent up to corrections of $O({\mathcal{H}}^2/k^2)$. Note that conditions $\left(ii\right)$ and $\left(iii\right)$ yield some restrictions respectively on the type of cosmological background and on the class of gauges. In order to show such approximate gauge independence, Ref. [161] argued that a good gauge choice to start with is the Newtonian gauge. The physical reason for that choice is that on small distances the Newtonian gauge reduces to the well-known Newtonian gravity. The practical reason is that calculations are simpler. Nevertheless, the argument provided in [161] does not strongly depend on the initial choice provided that it is suitable for subhorizon calculations. For example, one could have started with the synchronous gauge or the uniform Hubble gauge.

Starting from the Newtonian gauge we can relate the curvature perturbation $\Phi$ in any gauge, say ${\Phi }^G$, with that of the Newtonian gauge, say ${\Phi }^N$, by (187) which yields

$$\begin{array}{cc}\hfill {\Phi }_G& ={\Phi }_N+\mathcal{H}{T}_G+\frac{1}{3}\Delta L_G.\hfill \end{array}$$

(194)

Then, the class of suitable gauges for subhorizon physics is defined as those gauges in which

$$\begin{array}{c}\hfill {\Phi }_G(k\gg \mathcal{H})=O\left({\Phi }_N(k\gg \mathcal{H})\right).\end{array}$$

(195)

This means that the curvature perturbation need not be the same but have the same qualitative behaviour. We shall see that although this is a requirement on the scalar variable $\Phi$, it yields the gauge independence of $h_{ij}$. Using (195) into (194) we find conditions on the gauge parameters $T_G$ and $L_G$. For example, the gauge time parameter must be

$$\begin{array}{cc}\hfill T_G(k\gg \mathcal{H})& =O\left({\mathcal{H}}^{-1}{\Phi }_N(k\gg \mathcal{H})\right)\mathrm{or}\mathrm{higher}\mathrm{order}.\hfill \end{array}$$

(196)

Note that this abstract condition (195) actually includes most of the commonly used gauges such as the synchronous, the flat and the constant Hubble gauges [161]. However, it excludes the comoving slicing gauge [161]. The main reason for such exclusion is that the comoving slicing gauge is a choice where the spatial hypersurface is orthogonal to the flow lines of the fluid. This is a good choice for superhorizon scales where the flow lines are slowly changing with time, but it is a poor choice for subhorizon scales where fluid velocities oscillate. To convince ourselves of the gauge independence of (189), take the solution to ${\Phi }^N$ (74) and $h_{ij}^N$ (83) for constant w, which read

$$\begin{array}{c}\hfill {\Phi }_k^N\propto {\left(\frac{k}{\mathcal{H}}\right)}^{-2-b}\mathrm{and}{h}_k^N\propto {\left(\frac{k}{\mathcal{H}}\right)}^{-1-b}.\end{array}$$

(197)

For such solutions, the condition (196) translates into

$$\begin{array}{c}\hfill T_G=O\left(h_k^N\right)\mathrm{or}\mathrm{higher}\mathrm{order}.\end{array}$$

(198)

Then, using (188) with $L_G=0$ and in Fourier space we have

$$\begin{array}{c}\hfill h_k^G=h_k^N-\int \frac{d^{3}q}{{\left(2\pi \right)}^3}{e}^{ij}{q}_i{q}_j{T}_G\left(q\right)T_G\left(\right|\mathbf{k}-\mathbf{q}\left|\right),\end{array}$$

(199)

where we used that in the Newtonian gauge we have $E=\beta =0$. So when substituting (198) in (199) we conclude that for all gauges which satisfy (195) we have that

$$\begin{array}{c}\hfill h_k^G(k\gg \mathcal{H})=h_k^N(k\gg \mathcal{H})+O({\mathcal{H}}^2/k^2),\end{array}$$

(200)

which proves the approximate gauge independence of ${\rho }_{\mathrm{GW}}$ (189). It should be noted that while the work of [161] is not a solution per se, it shows when and why the predictions for induced GWs are well-defined and gives strong evidence that the calculations in this class of gauges are meaningful. In any case, the final word would be to find the GW detector response at second order.

Before ending this section, let us place some more attention to point $\left(ii\right)$, i.e., that the source term to induced GWs is not active. This means that the amplitude of scalar fluctuations has decayed enough so that they are no longer a significant source of induced GWs. This is often implied by taking the subhorizon limit $k\gg \mathcal{H}$ since induced GWs are mainly generated at horizon crossing, i.e., $k\sim \mathcal{H}$, and not afterwards. However, there is one exception: the dust dominated universe studied in Section 6. In a dust dominated universe, the Newtonian potential $\Phi$ is constant on all scales and, therefore, is constantly sourcing induced GWs. This is the reason why Ref. [149] finds such large gauge dependence of induced GWs in a dust dominated universe. Yet, we learned in Section 6 that the final spectrum of induced GWs is very sensitive to the transition from dust to radiation, where $\Phi$ changes from a constant to an oscillating function. Thus, one must follow the induced tensor modes generated in the dust dominated universe until they are deep inside the horizon in the radiation dominated era. It is from that moment on that the conditions $\left(i\right)$–$\left(iii\right)$ of [161] apply. Then we know that the spectrum of induced GWs is approximately gauge invariant.

8. Other GW Sources Related to PBH Formation

So far, we have focused on the GW induced by large primordial fluctuations. These large fluctuations collapse to form PBH if their rms amplitude is large enough. Once PBH form they may source other GWs independent of the induced GWs. Let us anticipate that the other GW counterparts are from PBHs mergers and graviton emission by Hawking evaporation. We briefly comment on these possibilities below.

On one hand, PBHs will form binaries in the early universe and eventually merge. The formation of such binaries is mainly due to the three body interaction with the nearest PBH [309]. These PBH binaries can be detected by the observation of resolved or unresolved merger events. In the former we would see the full GW waveform, e.g., as in the LIGO detections [310], and the latter will appear as a stochastic GW background. These make two additional possible GW counterparts to PBHs. Note that if PBHs are lighter than ${10}^{15}\mathrm{g}$ they evaporated by today and, therefore, we most likely cannot observe any resolved merger. PBHs heavier than ${10}^{15}\mathrm{g}$ have not yet evaporated and their binaries could be merging in the nearby universe. The GW spectrum from the unresolved mergers is studied in Refs. [142,311,312] while the estimates for the merger rates can be found in Refs. [174,309,313,314,315,316]. A very rough order of magnitude estimate of the frequency at which the GWs from the mergers of PBH binaries will show up is given by the frequency of the GWs emitted at the Innermost Stable Circular Orbit [185]

$$\begin{array}{c}\hfill f_{\mathrm{GW},\max }\approx 2f_{\mathrm{ISCO}}\approx 4.4\mathrm{kHz}\left(\frac{M_{\odot }}{M}\right),\end{array}$$

(201)

where M is the total mass of the binary and $M_{\odot }$ is a solar mass. For a monochromatic PBH mass function then $M=2M_{\mathrm{PBH},\mathrm{i}}$. This may change by a factor of a few depending on the redshift at the time of merger, so that $f_{\mathrm{GW},\max ,0}=f_{\mathrm{GW},\max }/(1+z)$. For PBHs with $M_{\mathrm{PBH},\mathrm{i}}\gg {10}^{15}\mathrm{g}$ we can assume that most of the mergers occur in the nearby universe and so the frequency (201) already gives a good intuition of the position of the peak in the SGWB from mergers and the chirp frequency of the GW waveform, see for example Ref. [142]. However, PBHs with $M_{\mathrm{PBH},\mathrm{i}}<5\times {10}^8\mathrm{g}$ have merged and evaporated long before BBN. In this case, we can assume that most of the mergers occur close to the evaporation time and, if PBH dominated the universe as in Section 6.2, we can take into account the redshift of the frequency (201) until today. A quick estimate yields

$$\begin{array}{c}\hfill f_{\mathrm{peak},0}^{\mathrm{merger}}\approx {10}^{15}\mathrm{Hz}{\left(\frac{M_{\mathrm{PBH},\mathrm{i}}}{{10}^4\mathrm{g}}\right)}^{1/2}{\left(\frac{g_{*}\left(T_{\mathrm{eva}}\right)}{106.75}\right)}^{1/4}{\left(\frac{g_{s*}\left(T_{\mathrm{eva}}\right)}{106.75}\right)}^{-1/3}{\left(\frac{g_H\left(T_{\mathrm{PBH}}\right)}{108}\right)}^{-1/2},\end{array}$$

(202)

which gives the right order of magnitude estimate for the position of the peak. For the detailed spectral shape in this case see Ref. [68]. The peak frequency of this SGWB is unfortunately too high to be observed by current laser interferometers. Note that the PBH mass range ${10}^9\mathrm{g}<M_{\mathrm{PBH},\mathrm{i}}<{10}^{14}\mathrm{g}$ is subject to tight constraints from entropy production, changes in the abundance of light elements, the extragalactic photon background and damping of CMB temperature anisotropies. For a detailed review see Ref. [176].

On the other hand, when PBHs evaporate they also emit “gravitons” by Hawking radiation. The typical frequency of such gravitons can be estimated by assuming that Hawking radiation follows a black body spectrum. Then, the peak frequency is roughly at $f_{\mathrm{peak}}\approx 2.82T_{\mathrm{PBH},\mathrm{i}}$ where $T_{\mathrm{PBH},\mathrm{i}}$ is given in Equation (164) in terms of the initial PBH mass $M_{\mathrm{PBH},\mathrm{i}}$. Note that the frequency $f_{\mathrm{peak}}$ is the peak frequency of the spectrum when PBH evaporate. Thus, we have to take into account that frequencies redshift proportional to $a^{-1}$ to calculate the frequency of these gravitons today. To do that, we will use the calculations of Section 6.2 which considers the case that PBH dominate the universe. Then, we can compute the evaporation temperature $T_{\mathrm{eva}}$ in terms of $M_{\mathrm{PBH},\mathrm{i}}$ by Equation (165). Inserting all numerical coefficients, we find that the graviton peak frequency today is given by

$$\begin{array}{c}\hfill f_{\mathrm{peak},0}^{\mathrm{graviton}}\approx {10}^{16}\mathrm{Hz}{\left(\frac{M_{\mathrm{PBH},\mathrm{i}}}{{10}^4\mathrm{g}}\right)}^{1/2}{\left(\frac{g_{*}\left(T_{\mathrm{eva}}\right)}{106.75}\right)}^{1/4}{\left(\frac{g_{s*}\left(T_{\mathrm{eva}}\right)}{106.75}\right)}^{-1/3}{\left(\frac{g_H\left(T_{PBH}\right)}{108}\right)}^{-1/2}.\end{array}$$

(203)

A more detailed derivation of the graviton spectrum and the frequency range can be found in Ref. [68]. Unfortunately, these are very high frequency GWs which are not observable by direct detection with the current capabilities. Nevertheless, they might be within range of future CMB experiments by looking at the contribution of these gravitons to the effective number of relativistic species at BBN [288,289,290,291,292]. Very interestingly if PBHs have a large spin, their imprint on the effective number is within range of the CMB-S4 experiment [291,292].

It is also important to notice that in some inflationary models involving axions coupled to gauge fields the enhancement in the primordial curvature perturbation is accompanied by an enhancement in the primordial tensor spectrum [88,100,317]. In this case the primordial GW signal is in general larger than the induced GW one. Thus, in this scenario the induced GW signal may be buried inside the primordial GWs. It is possible that a detailed observation of the spectrum and the PBH counterpart helps to distinguish such scenario.

9. Current and Future Observational Prospects

In this section, we briefly discuss the future observational prospects regarding induced GWs. We will not focus on the prospects for PBHs, although they are equally interesting. Nevertheless, let us mention the most interesting PBH windows as they provide good motivation to look for induced GWs in particular ranges. The first one is the LIGO binary black hole mergers [309,313,314,315]. Notably, although some models are already ruled out by current data [318] there seems to be evidence for a small fraction of PBHs in addition to astrophysical ones [319,320]. Interestingly, the corresponding induced GWs to the PBHs in the LIGO window fall in the PTA band.24 Coincidentally, the NANOGrav collaboration reported a possible SGWB in such range [18], but is yet to be confirmed. We discuss more on NANOGrav a bit later. Another interesting window is the planet mass PBHs [321] where the OGLE collaboration reported few microlensing events of planet mass objects [322]. Future microlensing observations such as with Subaru/HSC [323] will be able to confirm the results from OGLE. If we assume that radiation dominated early universe, planet mass PBHs have an induced GW counterpart that falls in between PTA and LISA. The last promising window is for PBHs in the asteroid mass range where there are practically no observational constraints yet [174,176] and therefore PBHs can make up for all dark matter. In this case, the induced GW counterpart falls right inside the LISA band. Thus, induced GW will be an essential probe of PBH as dark matter. For a review see Ref. [179]. Additionally, see Ref. [142] for a detailed discussion.

Let us now turn to the observational prospects of induced GWs. Although there are weak constraints from LIGO [324] on the induced SGWB [312,325,326], the most interesting part is the future prospects. On one hand, the absence of induced GWs directly translates into bounds on the primordial spectrum [30]. In this way, assuming radiation domination, we see from Ref. [30] that the bounds are roughly ${\mathcal{P}}_{\mathcal{R}}\sim {10}^{-5}$ for $k\sim {10}^6-{10}^8{\mathrm{Mpc}}^{-1}$ in the PTA band, ${\mathcal{P}}_{\mathcal{R}}\sim {10}^{-4}$ for $k\sim {10}^{11}-{10}^{14}{\mathrm{Mpc}}^{-1}$ in the LISA band and ${\mathcal{P}}_{\mathcal{R}}\sim {10}^{-3}$ for $k\sim {10}^{15}-{10}^{18}{\mathrm{Mpc}}^{-1}$ in the ET band. CMB spectral distortions will probe down to ${\mathcal{P}}_{\mathcal{R}}\sim {10}^{-8}$ for $k\sim 10-{10}^5{\mathrm{Mpc}}^{-1}$ [32,33]. However, one of the problems is that such future bounds depend on the shape of the induced GW spectrum. For example, if we only see part of the induced GW spectrum, e.g., the UV or IR tail but not the peak, then in order to extract information on the amplitude of primordial fluctuations we must extrapolate the results beyond the range of the GW detector by assuming some template [172]. On top of that, we have the issue that some induced GW templates might degenerate with other sources. As a good example for these potential problems, we proceed to discuss the various induced GW explanations of the NANOGrav results [18].

The NANOGrav collaboration reported a common process signal in the time residuals of the frequencies of pulsars [18]. While no quadrupolar nature of the signal has been found, it is interesting to consider the possibility that the signal measured is in fact the SGWB of induced GWs. The main features of the reported spectrum [18], in cosmologists terms, are that it has an amplitude of ${\Omega }_{\mathrm{GW},0}{h}^2\sim {10}^{-9}$ at a frequency of $2\times {10}^{-9}\mathrm{Hz}$ and a spectral index, that is ${\Omega }_{\mathrm{GW},0}\propto k^n$, that ranges from $n\sim 1/2$ to $n\sim -3/2$ within the 1-$\sigma$ bounds. This means that the NANOGrav data seems to prefer a rather flat GW spectrum. This is illustrated in Figure 9. Using the very rough estimate of Section 2, that is ${\Omega }_{\mathrm{GW},0}{h}^2\sim {10}^{-6}{\mathcal{P}}_{\mathcal{R}}^2$, we find that if the peak of the induced GW spectrum lies around $f\sim 2\times {10}^{-9}\mathrm{Hz}$ then the amplitude of primordial fluctuations that sourced such induced GWs is about ${\mathcal{P}}_{\mathcal{R}}\sim O\left({10}^{-2}\right)$. Such amplitude of primordial fluctuations is right at the value when PBH formation is interesting, and their masses are in the LIGO window.

The main complication is that the NANOGrav GW spectrum is rather flat and the induced GWs from a sharp peak are very peaked and decay as $k^3$ in the IR. This means that there are two possibilities: either the GWs were induced by a broad peak in the primordial spectrum (including a flat spectrum) [129,130,131,133,135] or the GWs were induced in non-radiation dominated universe with $w<1/3$ [132]. Below we list the main results of these works:

–

Ref. [129] found a negligible amount of PBH in LIGO region but the PBH could be the seeds of supermassive black holes.

–

Ref. [131] found a larger amount of PBH in the LIGO band than [129] and also pointed out that the SGWB from unresolved mergers could detectable in the future by LISA and ET.;

–

Ref. [130] showed that the induced GWs from a flat primordial spectrum could explain all dark matter;

–

Ref. [134] proposed a particular inflationary model with an axion-like curvaton which could explain the NANOGrav result as well as the LIGO events due to primordial non-gaussianities. Similar claims appear in Ref. [135] with a Higgs-like inflation with non-canonical kinetic terms;

–

In a different direction, Refs. [132,133] studied the induced GWs from a peaked primordial spectrum in a non-radiation dominated universe. They both find that a soft equation of state seems to be preferred. In particular, Ref. [132] use the IR tail of the induced GW spectrum from a very peaked primordial spectrum (see Section 5) to fit the NANOGrav results. They found that the 1$\sigma$ contours on the spectral tilt translate to bounds for the equation state parameter as $-0.1<w<0.05$;

–

Ref. [81] uses the UV tail of the induced GW spectrum from a broken power-law primordial spectrum. They concluded that the NANOGrav results imply a small non-gaussianity parameter $f_{NL}<0.7$;

It is important to note that the induced GWs are not the only explanations of the NANOGrav results. For example, among many others, possible candidates are cosmic strings [328,329,330,331] and first order phase transitions [332,333,334,335,336,337,338,339]. Within the next decades, we expect to have more information from PTA experiments. We also expect the GW detector LISA to launch. In order to distinguish these candidates to the NANOGrav data, we will need better resolution for a better Bayesian analysis [336] but crucially we will need other discriminators of cosmic GWs. In the case of induced GWs, if some of the LIGO binary black hole mergers are found to be PBHs, the NANOGrav signal would have a strong case to be induced GWs.

10. Summary of Main Formulas

We dedicate this section to recollect the main formulas scattered around the review and present them so that they are ready for use. We first give the GW spectrum measured today ${\Omega }_{\mathrm{GW},0}$ in terms of the spectral density of GWs at a pivot scale ${\Omega }_{\mathrm{GW},\mathrm{rh}}$. The subscript “rh” stands for reheating and refers to the time of instantaneous transition to the standard radiation dominated era. In the case where GWs are induced during radiation domination, the pivot time “rh” refers to the time when induced GWs of given wavenumber k are sufficiently inside the cosmological horizon to be treated as a radiation fluid in an expanding universe. With these clarifications, we can use the results of Section 4 and Section 5 to write

$$\begin{array}{c}\hfill {\Omega }_{\mathrm{GW},0}{h}^2=1.62\times {10}^{-5}\left(\frac{{\Omega }_{r,0}{h}^2}{4.18\times {10}^{-5}}\right)\left(\frac{g_{*}\left(T_{\mathrm{rh}}\right)}{106.75}\right){\left(\frac{g_{*s}\left(T_{\mathrm{rh}}\right)}{106.75}\right)}^{-4/3}{\Omega }_{\mathrm{GW},\mathrm{rh}}.\end{array}$$

(204)

We give the detailed expression of ${\Omega }_{\mathrm{GW},\mathrm{rh}}$ for a general equation of state in Section 10.1 and for the particular case of a dust universe with an instantaneous transition to radiation domination in Section 10.2.

10.1. General Equation of State

We studied in detail the derivation of ${\Omega }_{\mathrm{GW},\mathrm{rh}}$ for general expansion histories in Section 4. The main approximations used are the following: $\left(i\right)$ there is an epoch of constant w and $c_s$ after inflation, $\left(ii\right)$ the peak (or plateau) in the primordial spectrum enters the horizon during such era and $\left(iii\right)$ the transition to radiation domination is instantaneous. Note that in the case of standard radiation domination these assumptions are unnecessary. The presence of a transition to radiation domination introduces another relevant scale $k_{\mathrm{rh}}$, which is the last co-moving scale that entered the horizon at the transition. The GW spectrum for scales $k\gg k_{\mathrm{rh}}$ takes a general form given by

$${\Omega }_{\mathrm{GW},\mathrm{rh}}={\left(\frac{k}{k_{\mathrm{rh}}}\right)}^{-2b}{\int }_0^{\infty }dv{\int }_{|1-v|}^{1+v}du\mathcal{T}(u,v,b,c_s){\mathcal{P}}_{\mathcal{R}}\left(ku\right){\mathcal{P}}_{\mathcal{R}}\left(kv\right),$$

(205)

where we defined the transfer function as

$$\begin{array}{cc}\hfill \mathcal{T}(u,v,b,c_s)=& \mathcal{N}(b,c_s){\left(\frac{4v^2-{(1-u^2+v^2)}^2}{4u^2{v}^2}\right)}^2|1-y^2{|}^b\hfill \\ & \times \{{\left({\mathsf{P}}_b^{-b}\left(y\right)+\frac{b+2}{b+1}{\mathsf{P}}_{b+2}^{-b}\left(y\right)\right)}^2\Theta [c_s(u+v)-1]\hfill \\ & +\frac{4}{{\pi }^2}{\left({\mathsf{Q}}_b^{-b}\left(y\right)+\frac{b+2}{b+1}{\mathsf{Q}}_{b+2}^{-b}\left(y\right)\right)}^2\Theta [c_s(u+v)-1]\hfill \\ & +\frac{4}{{\pi }^2}{\left({\mathcal{Q}}_b^{-b}(-y)+2\frac{b+2}{b+1}{\mathcal{Q}}_{b+2}^{-b}(-y)\right)}^2\Theta [1-c_s(u+v)]\},\hfill \end{array}$$

(206)

and the numerical coefficient is given by

$$\begin{array}{c}\hfill \mathcal{N}(b,c_s)\equiv \frac{4^{2b}}{3c_s^4}{\Gamma }^4\left[b+\frac{3}{2}\right]{\left(\frac{b+2}{2b+3}\right)}^2{\left(1+b\right)}^{-2(1+b)}.\end{array}$$

(207)

In Equation (206) we have introduced the following notation for convenience:

$$b=\frac{1-3w}{1+3w}\mathrm{and}y=1-\frac{1-c_s^2{(u-v)}^2}{2c_s^{2}uv}.$$

(208)

These equations are general enough for $k\gg k_{\mathrm{rh}}$ that can be used for any shape of the primordial spectrum as long as the main feature, whether it is a peak or a plateau, enters during the $w=\mathrm{constant}$ epoch. If the transition to radiation domination is gradual, then we only expect a change in the IR tail of the spectrum for modes which entered during the relevant stages of the transition. Following with the instantaneous transition, it is important to note that for $b > 1$ the IR tail has a red tilt and the spectrum for $k\ll k_{\mathrm{rh}}$ becomes important. Nevertheless, a good approximation is to stop the spectrum at $k\sim k_{\mathrm{rh}}$ and match it to the typical IR scaling in radiation domination, i.e., ${\Omega }_{\mathrm{GW},\mathrm{rh}}\propto k^3$. Details on the spectrum for $k\ll k_{\mathrm{rh}}$ can be found in Section 4. It should be noted that (206) is not valid for $c_s=0$. The case of $c_s=0$ is studied in detail in Section 6 and a quick estimate is presented below in Section 10.2.

The transfer function (206) is expressed in terms of Associated Legendre functions which are provided in Appendix H in terms of Hypergeometric functions. This general form, while very useful, might be a hindrance to a quick implementation. For this reason, we present below particular examples of the transfer function (206). We choose the cases where $b\in \mathbb{Z}$ and $b\in \mathbb{Z}+\frac{1}{2}$ since then the associated Legendre functions (and the Hypergeometric functions) reduce to polynomials. Note that only when $b\in \mathbb{Z}$ logarithmic terms appear. We will first check the standard case of radiation domination ($b=0$, $w=1/3$) and we will later turn to other situations. We give the transfer function for a stiff fluid $b=-1/2$ ($w=1$), a soft fluid $b=1/2$ ($w=1/9$), a pressureless fluid $b=1$ ($w=0$) and negative equation of state fluid $b=2$ ($w=-1/9$). We keep the sound speed $c_s$ arbitrary. For an adiabatic perfect fluid, we have $c_s^2=w$ and for a canonical scalar field $c_s^2=1$.

10.1.1. Radiation Domination

The transfer function (206) for $b=0$ ($w=1/3$) reads

$$\begin{array}{cc}\hfill {\mathcal{T}}_{RD}(u,v,c_s,w=1/3)=& \frac{1y^2}{3c_s^4}{\left(\frac{4v^2-{(1-u^2+v^2)}^2}{4u^2{v}^2}\right)}^2\hfill \\ & \times \left\{\frac{{\pi }^2}{4}{y}^2\Theta [c_s(u+v)-1]+{\left(1-\frac{1}{2}yln\left|\frac{1+y}{1-y}\right|\right)}^2\right\},\hfill \end{array}$$

(209)

where we kept the variable y as it simplifies considerably the expressions. This recovers the result of Kohri and Terada [63] after taking into account all numerical factors in (7).

10.1.2. Stiff Fluid (Kinetic) Domination

Another notable case of interest is a stiff fluid with $b=-1/2$ ($w=1$). This is for example typical of quintessential inflation scenarios [340,341,342,343]. For a kinetic dominated period, the transfer function (206) is given by

$$\begin{array}{cc}\hfill & {\mathcal{T}}_{KD}(u,v,c_s,w=1)=\frac{4}{3\pi c_s^4|1-y^2|}{\left(\frac{4v^2-{(1-u^2+v^2)}^2}{4u^2{v}^2}\right)}^2\hfill \\ & \times \left\{\left(1+3y^2\right)\Theta [c_s(u+v)-1]+{\left(1-3y^2+3y\sqrt{|1-y^2|}\right)}^2\Theta [1-c_s(u+v)]\right\}.\hfill \end{array}$$

(210)

This coincides with the formulas derived in Ref. [64,74]. It has also been studied numerically in Ref. [75].

10.1.3. Soft Fluid Domination

The next case is that of a soft fluid with $b=1/2$ ($w=1/9$). This could be due to a scalar field rolling down an exponential potential. After some algebra, we find that

$$\begin{array}{cc}\hfill \mathcal{T}(u,v,c_s,w=1/9)& =\frac{2^8}{3^8\pi c_s^4}{\left(\frac{4v^2-{(1-u^2+v^2)}^2}{4u^2{v}^2}\right)}^2\hfill \\ & \times \{\left(4+45y^2\right)\Theta [c_s(u+v)-1]\hfill \\ & +{\left(y(3-10y^2)+(2+10y^2)\sqrt{|1-y^2|}\right)}^2\Theta [1-c_s(u+v)]\}.\hfill \end{array}$$

(211)

10.1.4. Pressure-Less Fluid Domination

We turn now to the case of a pressure-less fluid with $b=1$ ($w=0$). Note that this is not dust in the sense that $c_s\ne 0$. So that even though the fluid has no pressure, perturbations still propagate at a finite speed. This is the case of a scalar field rolling down an exponential potential and might also be close to the coherent oscillations of a scalar field around the bottom of the potential [344]. The kernel for $b=1$ is given by

$$\begin{array}{cc}\hfill \mathcal{T}(u,v,c_s,w=0)=& \frac{3^3{5}^2}{2^{14}{c}_s^4}{\left(\frac{4v^2-{(1-u^2+v^2)}^2}{4u^2{v}^2}\right)}^2\hfill \\ & \times \{\frac{{\pi }^2}{4}{(1-y^2)}^2{(1+3y^2)}^2\Theta [c_s(u+v)-1]\hfill \\ & +{\left(y(1-3y^2)-\frac{1}{2}(1+2y^2-3y^4)ln\left|\frac{1+y}{1-y}\right|\right)}^2\}.\hfill \end{array}$$

(212)

10.1.5. Negative EoS Fluid Domination

The last case we are interested in is a negative EoS parameter. When $-1/3<w<0$ the universe is decelerating very slowly, which yields to interesting slopes in the IR tail. For convenience we present the kernel for the case $b=2$ ($w=-1/9$), which reads

$$\begin{array}{cc}\hfill \mathcal{T}(u,v,c_s,w=-1/9)=& \frac{5^2{7}^2}{2^8{3}^9{c}_s^4}{\left(\frac{4v^2-{(1-u^2+v^2)}^2}{4u^2{v}^2}\right)}^2\hfill \\ & \times \{\frac{225{\pi }^2}{4}{(1-y^2)}^2{(1+y^2-2y^4)}^2\Theta [c_s(u+v)-1]\hfill \\ & +{\left(y(9+35y^2-30y^4)+\frac{15}{2}(1-3y^4+2y^6)ln\left|\frac{1+y}{1-y}\right|\right)}^2\}.\hfill \end{array}$$

(213)

10.2. Dust Domination

Here we provide a useful order of magnitude estimate for the following situation: $\left(i\right)$ the universe is dominated by dust $w=c_s^2=0$ when induced GWs are generated, $\left(ii\right)$ the transition to radiation domination is instantaneous and $\left(iii\right)$ the spectrum of fluctuations for the Newtonian potential $\Phi$ is given by

$$\begin{array}{c}\hfill {\mathcal{P}}_{\Phi }={\mathcal{A}}_{\Phi }{\left(\frac{k}{k_{\mathrm{UV}}}\right)}^{-n}\Theta (k_{UV}-k),\end{array}$$

(214)

where $n < 2$. This power spectrum includes a power-law generated during inflation and the case of PBH density fluctuations. The details of the derivation can be found in Section 6. With these assumptions, the spectrum of induced GWs is approximately given by

$$\begin{array}{c}\hfill {\Omega }_{\mathrm{GWs},\mathrm{rh}}(k\sim k_{\mathrm{UV}})\approx {\Omega }_{\mathrm{GWs}}^{\mathrm{peak}}{\left(\frac{k_{\mathrm{UV}}}{k_{\mathrm{rh}}}\right)}^{7-2n}\Theta (k_{\mathrm{UV}}-k),\end{array}$$

(215)

where the amplitude of the peak reads

$$\begin{array}{c}\hfill {\Omega }_{\mathrm{GWs}}^{\mathrm{peak}}\approx \frac{\pi }{9216\sqrt{3}}{\left(\frac{4}{3}\right)}^n{\left(\frac{k_{\mathrm{UV}}}{k_{\mathrm{rh}}}\right)}^7{\mathcal{A}}_{\Phi }^2.\end{array}$$

(216)

We can then use Equation (157) together with (204) to estimate the induced GW spectrum today. Note that if the UV cut-off $k_{\mathrm{UV}}$ is at very small scales, i.e., $k_{\mathrm{UV}}\gg k_{\mathrm{rh}}$, the amplitude of induced GW is very much enhanced and strong constraints apply. It is important to clarify that if $k_{\mathrm{UV}}$ is very large (larger than $k_{NL}$ in Equation (139)), some density fluctuations would enter the non-linear regime. Thus, the above result must be considered as a rough estimate. For detailed discussion we refer the reader to Section 6.

11. Conclusions

Gravitational waves induced by primordial fluctuations, so-called induced GWs, have been shown to be a very promising tool to complete our picture of inflation and the early universe. Although this effect was first pointed out 50 years ago, induced GWs together with primordial black holes are receiving serious attention since the first GW detection by LIGO [61]. They are even more interesting after the NANOGrav collaboration recently reported a possible stochastic GW signal [18]. Within the next couple of decades, LISA will launch, and we will have a new data on stochastic backgrounds of gravitational waves in a wide frequency range. Even in the worst-case scenario, where no conclusive induced GW signal is found, we will obtain new constraints on the physics of inflation and the early universe in completely unexplored regimes.

The field of cosmology with induced GWs is relatively new and still developing. In this review we have focused on the current analytical techniques and estimates to compute the induced GW spectrum. We believe they will be most useful in analysing future SGWB data. We have not dwelled on the details of PBH formation nor on the observational forecasts for future GW detectors, both of which deserve a separate review. In the next years there will be more studies on induced GWs; new and better analytical and numerical results might be developed. Until then, particularly interesting and important directions are: $\left(i\right)$ the impact of primordial non-gaussianities, not necessarily in the local shape, $\left(ii\right)$ possible discriminators of induced GWs to distinguish them from other sources, $\left(iii\right)$ the GW detector response at second order and $\left(iv\right)$ GWs from the non-linear regime in an early matter dominated stage. We hope this review will be useful to the next advancements in the cosmology of induced GWs.