1 Introduction

The standard model of cosmology requires adiabatic primordial fluctuations as initial conditions set in the very early universe [1,2,3]. Fluctuations are adiabatic when all fields filling the universe agreed to share a slicing of the spacetime where each of its separate energy density fluctuations vanish. In such choice of coordinates the energy density of those fields is homogeneous (by definition) and only the metric contains fluctuations. These are the so-called adiabatic curvature fluctuations. The reason for this sort of coherent initial conditions may simply be that one “primordial” field is responsible for all the fluctuations, e.g., because it decays into all the other fields. This is precisely the prediction from the simplest models of cosmic inflation and what Cosmic Microwave Background (CMB) observations confirmed [4].

If the initial conditions for primordial fluctuations are not adiabatic, they are said to be isocurvature (see, e.g., [5,6,7,8,9]). As the name implies, isocurvature requires that there is no initial adiabatic curvature perturbation. This means that there is no slicing of spacetime where each of the energy density fluctuations separately vanish (although there is always the slicing where the total energy density fluctuation is zero). And so, isocurvature initial conditions are related to relative energy density fluctuations (to be precise relative number density fluctuations). On scales larger than \(10\,\text {Mpc}\), CMB tells us that isocurvature primordial fluctuations may not account for more than \(1-10\%\) of the total fluctuations [2]. Since the measured amplitude of the power spectrum of primordial adiabatic fluctuations is about \(10^{-9}\), the power spectrum of isocurvature fluctuations may have an amplitude of less than \(10^{-10}\) on large scales.

The story changes when we consider scales smaller than \(1\,\,\text {Mpc}\), where CMB constrains no longer apply. For scales between \(1\,\text {Mpc}\) and \(1\,\text {pc}\), future CMB spectral distortions might be able to test isocurvature fluctuations [10, 11].Footnote 1 For scales smaller than a parsec our best bet to constrain the amplitude and nature of primordial fluctuations are primordial black holes (PBHs) and induced gravitational waves (GWs).Footnote 2 PBHs form from the collapse of large primordial fluctuations [16,17,18,19,20,21,22] and could explain the dark matter (DM) [23,24,25,26,27,28,29], some of the LIGO/VIRGO/KAGRA black hole mergers [30,31,32,33] and be the seeds of supermassive black holes [34, 35]. Induced GWs are a consequence of the non-linear nature of gravity: density fluctuations eventually lead to metric fluctuations, which include GWs as a secondary effect [36,37,38,39,40,41,42]. In some sense, the evolution of primordial fluctuations (e.g., the resulting density waves) yield an anisotropic stress which sources the secondary, or induced, GWs. See Refs. [43,44,45,46,47,48] and [49, 50] for recent reviews on PBHs and induced GWs, respectively. See also Ref. [51] for lecture notes on the collection of GW signatures of PBHs.

By extension of the CMB, the standard approach is to consider adiabatic initial conditions for PBHs and induced GWs. But, both PBHs and induced GWs may also be generated from isocurvature initial conditions [52, 53].Footnote 3 Isocurvature fluctuations generally occur in multi-field models of inflation [62, 63] and phase transitions [64] and may be consequence of the Poisson noise (or clustering) from the formation of compact structures in the very early universe, such as PBHs and solitonic structures like oscillons and cosmic strings [57, 65,66,67,68]. This opens the door to probes of the nature of primordial fluctuations and the formation of compact objects in the early universe. In this review, we will focus on isocurvature induced GWs. Details on PBHs formed by the collapse of primordial DM isocurvature fluctuations can be found in Ref. [52].

We are in an exciting time for cosmology with GWs. Recently, pulsar timing array (PTA) collaborations around the globe have announced tentative evidence of a GW background at nHz frequencies [69,70,71,72,73,74,75,76,77,78]. On plausible explanation is that they are induced GWs from primordial fluctuations [79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105,106,107,108,109,110,111,112,113,114,115,116] (or the merger of supermassive PBHs [117,118,119]). As a new application, here we will also explore the possibility that the PTA results are explained by DM isocurvature induced GWs. The PTA results are and will be complemented by other GW detectors at higher frequencies such as the LIGO/Virgo/KAGRA collaboration [120] and future GW detectors like \(\mu\)-Ares [121], LISA, Taiji [122, 123], TianQin [124], DECIGO [125, 126], Einstein Telescope (ET) [127], Cosmic Explorer (CE) [128], and Voyager [129, 130].

This work is organized as follows. In Section 2, we present an overview of the basic equations and the notion of adiabatic and isocurvature initial conditions. The content of this section is based on Refs. [7, 8]. Then, we discus, we review GWs induced by DM isocurvature, based on the results of Ref. [53]. In Section 4, we consider the isocurvature due to PBH number density fluctuations in the PBH reheating scenario. Here, we will base our discussions on Refs. [66, 131,132,133]. Lastly, we discuss more applications in Section 5. Some details of the calculations can be found in the appendices and in the aforementioned references. We work in reduced Planck units where \(c=\hbar =1\) and \(M_{\text {pl}}=(8\pi G)^{-1/2}=1\).

2 Basic equations

We start by deriving the general formalism for isocurvature induced GWs. Since isocurvature fluctuations result from relative number density fluctuations, we should at least consider two fluids filling the primordial universe. For simplicity, we assume that after cosmic inflation the universe is dominated by relativistic particles, so-called radiation, and that there is a small but non-zero fraction of non-relativistic particles, let us call it matter. The energy momentum tensors of radiation and matter are respectively given by

\(\begin{aligned} T_{r\mu \nu }&=(\rho _r+p_r)u_{r\mu }u_{r\nu }+p_r g_{\mu \nu }\,,\end{aligned}\)

\(\begin{aligned} T_{m\mu \nu }&=\rho _m u_{m\mu }u_{m\nu }\,, \end{aligned}\)

where \(g_{\mu \nu }\) is the metric, \(\rho\) and p respectively are the energy density and pressure and \(u_{\mu }\) the fluid 4-velocity. The subscripts “r” and “m” respectively refer to radiation and matter. In particular, we have that \(p_m=0\) and \(p_r=\rho _r/3\). For the metric \(g_{\mu \nu }\), we take a perturbed Friedmann-Lemaître-Robertson-Walker (FLRW) universe. For convenience (the equations for induced GWs are simplest), we work in the Newtonian gauge in which the metric reads

\(\begin{aligned} ds^2=a^2(\tau )\left[ -(1+2\Psi )d\tau ^2+((1+2\Phi )\delta _{ij}+h_{ij})dx^idx^j\right] \,, \end{aligned}\)

where a is the scale factor, \(\tau\) is conformal time, \(\Psi\) and \(\Phi\) are the gravitational potentials and \(h_{ij}\) the tensor perturbations (or we may say GWs). The dynamics of the scale factor are dictated by the Friedmann equations, which we present in Appendix 1. In the same appendix, we describe the notation for matter perturbations as well.

Most important to us is the fact that energy conservation (\(\nabla ^\mu T_{\mu \nu }=0\)) at the background level requires that \(\rho _r\propto a^{-4}\) and \(\rho _m\propto a^{-3}\) (also see Appendix  1). From the different dilution of the energy densities, we see that

\(\begin{aligned} \frac{\rho _m}{\rho _r}\propto a\,, \end{aligned}\)

and, therefore, matter fields eventually dominate the energy density of the universe, assuming they do not decay. Thus, if we call \(\beta\) to the initial fraction of matter, namely

\(\begin{aligned} \beta =\frac{\rho _{m,i}}{\rho _{r,i}}<1\,, \end{aligned}\)

where “i” refers to the initial time, the universe will be matter dominated after \(a/a_i>\beta ^{-1}\). Actually, there is an exact analytical solution for the scale factor in the radiation-matter universe. It reads [7, 134]

\(\begin{aligned} \frac{a(\tau )}{a_{\text {eq}}}=2\left( \frac{\tau }{\tau _*}\right) +\left( \frac{\tau }{\tau _*}\right) ^2\,, \end{aligned}\)

where \((\sqrt{2}-1)\tau _*=\tau _{\text {eq}}\). The subscript “eq” means matter-radiation equality, i.e.,  the time when \(\rho _r=\rho _m\). It is easy to check that Eq. (6) goes from the radiation dominated universe where \(a\sim \tau\) to the matter dominated universe with \(a\sim \tau ^2\). This is for the moment all we need to understand the source of induced gravitational waves.

2.1 What sources secondary gravitational waves?

Before entering into the computational details, let us qualitatively understand what is the main source of secondary gravitational waves (at least in the Newton gauge). Let us formally start with Einstein equations, that is

\(\begin{aligned} G_{\mu \nu }=T_{m\mu \nu }+T_{r\mu \nu }\,, \end{aligned}\)

where \(G_{\mu \nu }\) is the Einstein tensor and for simplicity we set \(M_{\text {pl}}=(8\pi G)^{-1/2}=1\). The linear equations of motion for the tensor modes correspond to the transverse-traceless projection of the spatial-spatial components. If we call \({\mathcal {P}}_{ij}\,^{ab}\) the transverse-traceless projector (which can be found in, e.g., Refs. [49, 132]), then we schematically have at linear order

\(\begin{aligned} {\mathcal {P}}_{ij}\,^{ab}G^{(1)}_{ab}={\mathcal {P}}_{ij}\,^{ab}\left( T^{(1)}_{mab}+T^{(1)}_{rab}\right) \quad \Rightarrow \quad h_{ij}^{\prime \prime }+2{\mathcal {H}}h'_{ij}+\Delta h_{ij}=0\,, \end{aligned}\)

where \({\mathcal {H}}=a'/a\), \('=d/d\tau\) and the superscript (1) refers to linear perturbations. We will use a superscript (2) to refer to second order perturbations. Equation (8) basically tells us that since \(T_{\mu \nu }\) has no tensor component at linear order (because we assumed a perfect fluid with no anisotropic stress), gravitational waves propagate freely. If we include second order scalar terms though, we find that

\(\begin{aligned} {\mathcal {P}}_{ij}\,^{ab}G^{(1)}_{ab}={\mathcal {P}}_{ij}\,^{ab}\left( -G^{(2)}_{ab}+T^{(2)}_{mab}+T^{(2)}_{rab}\right) \quad \Rightarrow \quad h_{ij}^{\prime \prime }+2{\mathcal {H}}h'_{ij}+\Delta h_{ij}={\mathcal {P}}_{ij}\,^{ab}S_{ab}\,, \end{aligned}\)

which after some simplifications the source term is given by [132]

\(\begin{aligned} S_{ij}&=4\partial _i\Phi \partial _j\Phi +2a^2\rho _m \partial _iv_{m}\partial _jv_{m}+2a^2(\rho _r+p_r)\partial _iv_{r}\partial _iv_{r}\,. \end{aligned}\)

Note that we selected the scalar component of the perturbation of the spatial velocity in Eq. (1), that is we took \(u_i=a\partial _iv\). Also in Eq. (9), we considered the second order expansion of \(G^{(2)}_{ab}\) as a source (or backreaction) to the linear equations and, as such, we moved it to the right hand side. We then used that \(\Psi +\Phi =0\) in the presence of no anisotropic stress (see Appendix 1).

Let us discuss the secondary source to GWs, Eq. (10), within the big picture. First, we see that from \(G^{(2)}_{ab}\) we obtain gradients of \(\Phi\). So, one source of secondary gravitational waves are curvature (metric) fluctuations. Second, from \(T_{\mu \nu }\) only the spatial component of the fluid velocity contributes. So velocity flows also generate GWs (this is the main, intuitive, source of GWs inside the Hubble horizon). However, this does not tells us much yet about isocurvature initial conditions. To make it more intuitive, let us introduce the total spatial velocity (the one seen by the linear Einstein Equations and so linear metric fluctuations) and the relative velocity which respectively read [135]

\(\begin{aligned} \rho V=(\rho _m+\rho _r) V=\rho _m v_m+(\rho _r+p_r) v_r\quad {\text {and}}\quad V_{\text {rel}}=v_m-v_r \,. \end{aligned}\)

In terms of these variables we find that Eq. (10) is given by

\(\begin{aligned} S_{ij}&=4\partial _i\Phi \partial _j\Phi +2a^2\frac{\rho ^2}{\rho +p}\left( \partial _i V\partial _j V+\frac{\rho _m(\rho _r+p_r)}{\rho ^2}\partial _i V_{\text {rel}}\partial _j V_{\text {rel}}\right) \,, \end{aligned}\)

where \(p=p_m+p_r=p_r\). Now, it is clear that the contribution from the relative velocity is always suppressed by the energy density of the subdominant field. So unless we are considering scales that enter the horizon close to the matter-radiation equality, we may neglect the relative velocities. With this, we conclude that secondary gravitational waves are mainly sourced by the dominant fluid, which is the main source of curvature fluctuations. Any contribution from isocurvature initial conditions (= no initial curvature or \(\Phi _i=0\)) must then be suppressed by a factor \(\rho _m/\rho\) in the early stages.

A natural question then arises: why should we consider isocurvature induced GWs? But this question misses that the point that isocurvature is a matter of initial conditions. As the system evolves, isocurvature fluctuations are transferred into curvature fluctuations, the transfer being complete after matter-radiation equality [7]. So the answer is that isocurvature induced GWs can be important when:

  • Isocurvature fluctuations are large enough to compensate for the suppression factor \(\rho _m/\rho\), or,

  • Non-relativistic particles dominate the universe in an early matter dominated era with isocurvature induced curvature fluctuations.

We will study the first case in Section 3 for DM fluctuations and the second in Section 4 for the PBH dominated early universe. Note that for simplicity we consider the case of standard cold DM. However, as far as the GW computations are concerned, the results may also apply to the case of scalar dark matter, like ultra-light (or wave) dark matter (see Refs. [136, 137] for recent reviews on ultra-light dark matter). We briefly discuss this further in Section 3.2. We now proceed with the evolution of the curvature fluctuations and the formal solutions to induced GWs.

2.2 Evolution of curvature perturbation: adiabatic vs isocurvature initial conditions

The general notion of curvature and isocurvature fluctuations is better understood using the gauge invariant definition of the curvature perturbation on uniform density slices, which is given by, see, e.g., [8, 138],

\(\begin{aligned} \zeta =\phi -{\mathcal {H}}\frac{\delta \rho }{\rho '}\,. \end{aligned}\)

Here, we used \(\phi\) to denote the spatial curvature perturbation which corresponds to \(\Phi\) in the Newton gauge. In the uniform density slicing where \(\delta \rho =0\) we have that \(\zeta =\phi\). The convenience of using \(\zeta\) is that it is conserved on superhorizon scales. One may also define a curvature fluctuation for each fluid, namely

\(\begin{aligned} \zeta _y=\phi -{\mathcal {H}}\frac{\delta \rho _y}{\rho '_y}\,, \end{aligned}\)

where in our case, \(y=\{r,m\}\). And, from the individual definitions, the notion of isocurvature follows as the relative individual curvature fluctuations, namely

\(\begin{aligned} S=3(\zeta _m-\zeta _r)=\frac{\delta \rho _m}{\rho _m}-\frac{\delta \rho _r}{\rho _r+p_r}=\frac{\delta \rho _m}{\rho _m}-\frac{3}{4}\frac{\delta \rho _r}{\rho _r}\,, \end{aligned}\)

where we used that \(\rho '_y+3{\mathcal {H}}(\rho _y+p_y)=0\) and that \(p_r=1/3 \rho _r\). The physical intuition behind such formal definitions, Eqs. (13) and (15), is the following. Curvature fluctuations are metric fluctuations even in the time slicing where there are no total density fluctuations (\(\delta \rho =0\)). Adiabatic initial conditions then require no initial isocurvature, that is \(S_i=0\), and \(\zeta _i\ne 0\). Isocurvature initial conditions correspond to the case when there are no such curvature fluctuations but there are relative density fluctuations (i.e., we may have \(\zeta _i=0\) but \(S_i\ne 0\)).

To understand this also in terms of the Newton gauge variables, it is convenient to look at the 0-0 Einstein Equation on superhorizon scales (see Appendix 1), which reads

\(\begin{aligned} 6{\mathcal {H}}_i^2\Phi _i=a_i^2\left( \delta \rho _{m,i}+\delta \rho _{r,i}\right) =a_{i}^2\delta \rho _i\,, \end{aligned}\)

where we imposed that initially \(\Phi '_i=0\), we neglected gradient terms and we used the subscript “i” to denote evaluation at the initial time. In the Newton gauge, initial curvature fluctuations \(\Phi _i\) are then proportional to \(\delta \rho _i\). Thus, adiabatic initial conditions correspond to \(S_i=0\) and \(\Phi _i\propto \delta \rho _i\ne 0\) and isocurvature initial conditions to \(\Phi _i\propto \delta \rho _i=0\) and \(S_i\ne 0\). An interesting perspective from these definitions is that adiabatic initial conditions are set by the fluid with dominant energy density (since \(S=0\) the dominant \(\rho\) has also dominant \(\delta \rho\)) while isocurvature initial conditions are mainly given by the sub-dominant field (since \(\delta \rho _m+\delta \rho _r=0\) and S depends inversely in \(\rho _m\)).

The closed system of equations for curvature-isocurvature fluctuations in Fourier modes is given by (see Appendix 1 or [53, 132] for details)

\(\begin{aligned} \Phi ^{\prime \prime }+3{\mathcal {H}}(1+c_s^2)\Phi '+({\mathcal {H}}^2(1+3c_s^2)+2{\mathcal {H}}')\Phi +c_s^2k^2\Phi =\frac{a^2}{2}\rho _mc_s^2S\,, \end{aligned}\)


\(\begin{aligned} S^{\prime \prime }&+ 3{\mathcal {H}}c_s^2S'-\frac{3}{2a^2\rho _r}c_s^2{k^4\Phi }+\frac{3\rho _m}{4\rho _r}c_s^2k^2 S=0\,, \end{aligned}\)

where we defined as usual

\(\begin{aligned} c_s^2\equiv \frac{4}{9}\frac{\rho _r}{\rho _m+4\rho _r/3}\,. \end{aligned}\)

The relative velocity \(V_{\text {rel}}\) can be computed from S by

\(\begin{aligned} V_{\text {rel}}=S'/k^2\,. \end{aligned}\)

Let us note here that in the matter dominated universe where \(c_s^2\rightarrow 0\) the general “growing mode” solution to Eq. (17) is \(\Phi ={\text {constant}}\). The precise value of the constant will be set by the evolution during the radiation dominated universe. We present first the general solutions for superhorizons scales and later the solutions for general k in the radiation dominated universe. We will denote with \(S_i\) and \(\Phi _i\) the initial values of isocurvature and curvature in the far past, formally when \(a\rightarrow 0\). The conclusions do not change if a has a non-zero value as the solutions are attractors [67].

Superhorizon fluctuations (\(k\ll {\mathcal {H}}\)):

If we drop all terms containing k in Eqs. (17) and (18), and re-write them using a as a time variable, we find that the general solutions are given by [7, 134, 139]

\(\begin{aligned} S(\xi )&=S_i\,,\end{aligned}\)

\(\begin{aligned} \Phi (\xi )&=\Phi _i \left( \frac{8}{5\xi ^3}\left( \sqrt{1+\xi }-1\right) - \frac{4}{5\xi ^2} + \frac{1}{5\xi }+\frac{9}{10}\right) \nonumber \\&+S_i\left( \frac{16}{5\xi ^3} \left( 1-\sqrt{1+\xi }\right) +\frac{8}{5 \xi ^2} -\frac{2}{5\xi }+\frac{1}{5}\right) \,, \end{aligned}\)

where for compactness we defined \(\xi =a/a_{\text {eq}}\). We see that while isocurvature remains constant on superhorizon scales, curvature fluctuations \(\Phi\) are not. For the adiabatic component we have that

\(\begin{aligned} \Phi _{\text {ad}}(a)\approx \Phi _i\times \left\{ \begin{array}{ll} 1 &{}(a\ll a_{\text {eq}})\\ \frac{9}{10} &{}(a\gg a_{\text {eq}}) \end{array} \right. \,, \end{aligned}\)

where the factor 9/10 actually comes from the conservation on \(\zeta\) and the relation \(\zeta =\tfrac{5+3w}{3+3w}\Phi\) for constant \(w=p/\rho\). For isocurvature fluctuations, we instead have

\(\begin{aligned} \Phi _{\text {iso}}(a)\approx S_i\times \left\{ \begin{array}{ll} \frac{1}{8}\frac{a}{a_{\text {eq}}} &{}(a\ll a_{\text {eq}})\\ \frac{1}{5} &{}(a\gg a_{\text {eq}}) \end{array} \right. \,. \end{aligned}\)

Thus, while we have vanishing initial curvature perturbation, i.e., \(\Phi _{\text {iso}}(a\rightarrow 0)\rightarrow 0\), it grows as \({a}/a_{\text {eq}}\) and saturates to 1/5 of the initial isocurvature after matter-radiation equality.

Fluctuations during radiation domination (\(a\ll a_{\text {eq}}\)):

To study the evolution of general fluctuations during radiation domination, we shall take the limit \(\tau \ll \tau _{\text {eq}}\) (\(a\ll a_{\text {eq}})\)in Eqs. (17) and(18). At leading order, they are given by [53]

\(\begin{aligned} \frac{d^2\Phi }{dx^2}+\frac{4}{x}\frac{d\Phi }{dx}+\frac{1}{3}\Phi +\frac{1}{4\sqrt{2}\kappa x}\left( x\frac{d\Phi }{dx}+(1-x^2)\Phi -2S\right) \simeq 0\,, \end{aligned}\)


\(\begin{aligned} \frac{d^2S}{dx^2}&+\frac{1}{x}\frac{dS}{dx}-\frac{x^2}{6}\Phi -\frac{1}{2\sqrt{2}\kappa }\left( \frac{dS}{dx}-\frac{x}{2}S-\frac{x^3}{12}\Phi \right) \simeq 0\,, \end{aligned}\)

where we defined for compactness

\(\begin{aligned} x=k\tau \quad {\text {and}}\quad \kappa =\frac{k}{k_{\text {eq}}}\,. \end{aligned}\)

In these new variables, the limit of interest is given by \(x\ll \kappa\) (or \(k_{\text {eq}}\tau \ll 1\)). Due to the length of the solutions we treat the initial adiabatic and initial isocurvature cases separately below.

For initial curvature fluctuations, \(\Phi _i\ne 0\) and \(S_i=0\), we solve the leading order terms in Eqs. (25) and (26), namely we first solve the homogeneous equation for \(\Phi\) and plug it in in the equation for S. This yields

\(\begin{aligned} \Phi _{\text {ad}}(x/\kappa \ll 1)&\approx 3{\Phi _i} \frac{j_1(c_s x)}{c_s x}+\mathcal {O}\left( x/\kappa \right) \,,\end{aligned}\)

\(\begin{aligned} S_{\text {ad}}(x/\kappa \ll 1)&\approx 9\Phi _i\left( \gamma _E-\frac{1}{2}+\frac{1}{2}\cos (c_sx)-{\text {Ci}}(c_sx)+\log (c_sx)\right) +\mathcal {O}\left( x/\kappa \right) \,, \end{aligned}\)

where \(\gamma _E\approx 0.577\) and \({\text {Ci}}(x)\) is the cosine integral function. Note that Eq. (28) is the standard solution for adiabatic perturbations in the radiation universe. \(\Phi\) is first constant and decays as \(x^{-2}\) once a given mode enters the sound horizon, that corresponds to \(c_sx>1\). We also see that S is negligible for \(c_sx<1\) as it is proportional to \(x^4\) but grows as \(\log (c_s x)\) for \(c_sx>1\). Such logarithmic grows is due to the fact that matter perturbations grow logarithmically during the radiation dominated universe [140]. It should be noted that one may also use Eq. (29) to compute the next order solution to \(\Phi\). We refer the interested reader to Ref. [7] for general solutions.

For the initial isocruvature case, \(\Phi _i=0\) and \(S_i\ne 0\), we follow Ref. [53]. In this case, we see that we can expand the solution of \(\Phi\) and S in powers of \(\kappa ^{-1}\) for \(\kappa \gg 1\) and in powers of \(x/\kappa\) for \(\kappa \ll 1\), as, e.g., \(S=S_i+S_1+...\) and \(\Phi =\Phi _1+...\) etcetera. The leading contribution to Eq. (25) \(\Phi\) is a then a constant S and we may also compute the effects of the leading solution to \(\Phi\) to the next leading solution to S. Doing so we find [53]

\(\begin{aligned} \Phi _{\text {iso}}(x/\kappa \ll 1)&\approx \frac{3S_i}{2\sqrt{2}\kappa } \frac{1}{x^3} \left( 6+x^2-2\sqrt{3}x\sin (c_sx)-6\cos (c_sx)\right) +\mathcal {O}\left( x/\kappa \right) ^{2}\,,\end{aligned}\)

\(\begin{aligned} S_{\text {iso}}(x/\kappa \ll 1)&\approx S_i+\frac{3S_i}{2\sqrt{2}\kappa }\left( x+\sqrt{3}\sin (c_sx)-2\sqrt{3}{\text {Si}}(c_sx)\right) +\mathcal {O}\left( x/\kappa \right) ^{2}\,, \end{aligned}\)

where \({\text {Si}}(x)\) is the sine integral function. Looking at \(x\ll 1\) we see that initially the curvature perturbation grows as \(\Phi _{\text {iso}}\propto x\), reaches a maximum at around \(c_sx\sim 1\) and then decays as \(x^{-2}\). It is also interesting to note that for \(c_sx\gg 1\) we have \(\Phi _{\text {ad}}\supset \sin (c_sx)\) and \(\Phi _{\text {iso}}\supset \cos (c_sx)\), recovering the well-known result that adiabatic and isocurvature initial conditions give an out of phase curvature fluctuations. Isocurvature S is constant for \(c_sx<1\) and then grows with x for \(c_sx>1\). Interestingly, it is possible that S reaches a high enough amplitude for PBHs to form [52]. Although we will not explore this possibility in this work, let us write down the time when the local density of matter is larger than that of radiation. This happens at

\(\begin{aligned} \tau _{\text {NL}}=\frac{\sqrt{2}}{k_{\text {eq}} S_i}\,. \end{aligned}\)

At times \(\tau >\tau _{\text {NL}}\) we can no longer trust of linear solutions. Nevertheless, this only occurs when \(S_i\) is very large. In most situations regarding induced GWs we may consider that \(\tau _{\text {NL}}\) is late enough such that it does not affect the production of GWs as curvature perturbations already decayed significantly. For example, requiring that the non-linear regime occurs deep inside the horizon (i.e., \(x_{\text {NL}}>1\)) translates into an upper bound on the initial isocurvature, namely

\(\begin{aligned} S_i<\sqrt{2}\kappa \,, \end{aligned}\)

for a given isocurvature mode with wavenumber k.

Fluctuations during matter domination (\(a\gg a_{\text {eq}}\)):

Well inside matter domination curvature fluctuations on all scales becomes constant. The amplitude of such fluctuations is then determined by whether they were superhorizon or subhorizon during the radiation domination epoch. For \(k\ll k_{\text {eq}}\) its value is given by Eqs. (23) and (24), respectively, for adiabatic and isocurvature initial conditions. For \(k\gg k_{\text {eq}}\) the amplitude has a suppression factor proportional to \((k/k_{\text {eq}})^{-2}\), which comes from the \(x^{-2}\) decay during the radiation dominated phase. For the analytical approximations, we refer the reader to the work of Kodama and Sasaki [7], although in our simplified set up the numerical prefactors could be refined by matching at matter-radiation equality. As it is not crucial for our purposes we leave it as an exercise. Most important for the present review is the result of the curvature perturbation for isocurvature initial conditions at matter domination, which is given by [7]

\(\begin{aligned} \Phi _{\text {iso}}(a\gg a_{\text {eq}})\approx S_i\times \left\{ \begin{array}{ll} \frac{1}{5} &{}(k\ll k_{\text {eq}})\\ \frac{3}{4}\left( \frac{k}{k_{\text {eq}}}\right) ^{-2} &{}(k\gg k_{\text {eq}}) \end{array} \right. \,. \end{aligned}\)

We show the results of numerical integration in Fig. 1. In the numerical results we find that the coefficient for \(k\gg k_{\text {eq}}\) in Eq. (34) is close to 1. For easier comparison with the literature though, we maintain the coefficient of Eq. (34) as it only introduces small errors.

Fig. 1
figure 1

Numerical solutions to the evolution of curvature fluctuations for adiabatic and isocurvature initial conditions in terms of \(x=k\tau\), respectively in purple and orange lines. We normalized the amplitude of the initial conditions to \(\Phi _i=1\) and \(S_i=1\). On the left, we show the behavior on scales with \(k\ll k_{\text {eq}}\). For the numerical solution, we chose \(k=10^{-3}k_{\text {eq}}\). See how both lines saturate after matter-radiation equality to a constant. For the adiabatic case, the constant is 9/10 and for isocurvature 1/5. On the right, we instead show \(k\gg k_{\text {eq}}\), in particular \(k=500k_{\text {eq}}\). In this case the lines also saturate to a constant after matter-radiation equality. However, for \(k\gg k_{\text {eq}}\), the amplitude after \(\tau _{\text {eq}}\) is of the order of \((k_{\text {eq}}/k)^2\) in both cases

2.3 General formulation of isocurvature induced gravitational waves

We proceed with the formal solution to induced gravitational waves. Our starting point is the equations of motion for secondary GWs (9) which is given by

\(\begin{aligned} h_{ij}^{\prime \prime }+2{\mathcal {H}}h'_{ij}+\Delta h_{ij}={\mathcal {P}}_{ij}\,^{ab}S_{ab}\,, \end{aligned}\)

where \(S_{ij}\) (12) using the linear Einstein Equations is given by

\(\begin{aligned} S_{ij}&=4\partial _i\Phi \partial _j\Phi +6c_s^2\frac{\rho }{\rho _r}\partial _i\left( \frac{\Phi '}{{\mathcal {H}}} +\Phi \right) \partial _j\left( \frac{\Phi '}{{\mathcal {H}}}+\Phi \right) +6a^2c_s^2\rho _m\partial _iV_{\text {rel}}\partial _jV_{\text {rel}}\,. \end{aligned}\)

It will be more convenient to work in Fourier modes, where we give a primordial initial amplitude of \(S_k(0)\) and/or \(\Phi _k(0)\) and describe the evolution of a given mode by a transfer function, say \(T_\Phi (k\tau )\). In the notation of the previous subsections, \(S_k(0)\) should be understood as \(S_i\) for a given mode k. The same applies to \(\Phi\). We will also only consider the case of either adiabatic or isocurvature initial condition and, therefore, we neglect any cross term coming from the gradient squared terms. Doing so Eq. (35) reads

\(\begin{aligned} h^{\prime \prime }_{\textbf{k},\lambda }+2{\mathcal {H}} h'_{\textbf{k},\lambda }+k^2 h_{\textbf{k},\lambda }={\mathcal {S}}_{\textbf{k},\lambda }\,, \end{aligned}\)


\(\begin{aligned} {\mathcal {S}}_{\textbf{k},\lambda }=4\int \frac{d^3q}{(2\pi )^3} e_\lambda ^{ij}(k)q_iq_j \left\{ \begin{array}{c} \Phi _{\textbf{q}}(0)\Phi _{|\textbf{k}-\textbf{q}|}(0)\\ S_{\textbf{q}}(0)S_{|\textbf{k}-\textbf{q}|}(0) \end{array} \right\} f(\tau , q,|\textbf{k}-\textbf{q}|)\,, \end{aligned}\)

where \(\Phi _{\textbf{q}}(0)\) and \(S_{\textbf{q}}(0)\) respectively refer to initial adiabatic or initial isocurvature fluctuations, and we defined

\(\begin{aligned} f(\tau , q,|\textbf{k}-\textbf{q}|)=&f(x,k, u,v)=T_\Phi (v x)T_\Phi (u x)+ \frac{3}{2}c_s^2 a^2\rho _m T_{V_{\text {rel}}}(v x )T_{V_{\text {rel}}}(u x)\nonumber \\&+\frac{3}{2}c_{s}^{2}\left( 1+\frac{\rho _{m}}{\rho _{r}}\right) \left( T_{\Phi }(v x)+\frac{T'_{\Phi }(v x)}{\mathcal {H}}\right) \left( T_{\Phi }(u x)+\frac{T'_{\Phi }(u x)}{\mathcal {H}}\right) \,. \end{aligned}\)

In Eq. (39), we introduced for later use the notation

\(\begin{aligned} v=q/k\quad ,\quad u=|\textbf{k}-\textbf{q}|/k\,. \end{aligned}\)

Also, in Eq. (39), \(T_{V_{\text {rel}}}(x)\) is the transfer function of \(V_{\text {rel}}\) which is not important in the regimes of interest, so we will not consider it in the following sections. Note that the transfer functions that appear in Eq. (39) have to be properly chosen according to whether we are considering adiabatic or isocurvature initial conditions.

Assuming that the isocurvature fluctuations are Gaussian we arrive at a compact expression for the spectral density of induced GWs [53, 141, 142], namely

\(\begin{aligned} \Omega _{\text {GW,c}}(k)=\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(x_c,k,u,v)}{{\mathcal {P}}_{S}(ku)}{{\mathcal {P}}_{S}(kv)}\,, \end{aligned}\)

where \({\mathcal {P}}_{S}(k)\) is the initial dimensionless spectrum of isocurvature fluctuations, defined by

\(\begin{aligned} \langle S_\textbf{k}(0)S_{\textbf{k}'}(0)\rangle =\frac{2\pi ^2}{k^3}{\mathcal {P}}_S(k)\times (2\pi )^3\delta ^{(3)}\left( \textbf{k}+\textbf{k}'\right) \,. \end{aligned}\)

For adiabatic initial conditions \({\mathcal {P}}_{S}(k)\) should be replaced by \({\mathcal {P}}_{\Phi }(k)\) at the initial time. In Eq. (41) the notation \(x_c\) refers to a time when the GW is well inside the horizon such that the spectral density is constant [143]. In Eq. (41) we introduced the kernel which incorporates the information from the transfer functions and it is given by

\(\begin{aligned} I(x,k,u,v)\equiv x \,\int _{x_i}^{x}d\tilde{x} \,G(x,\tilde{x})f(\tilde{x},k,u,v)\,, \end{aligned}\)

where \(G(x,\tilde{x})\) is the tensor mode’s Green’s function.

As we will be mostly interested in the induced GWs generated during an radiation dominated phase, either the one preceding matter domination or the one after an early matter domination, we restrict ourselves to the radiation dominated universe. In particular, the Green’s function in the radiation dominated universe reads

\(\begin{aligned} G(x,\tilde{x})= \frac{a(\tilde{x})}{a( x )}\left( \sin x \cos \tilde{x}-\cos x \sin \tilde{x} \right) \,. \end{aligned}\)

Then, we may split the kernel (43) into a sine and cosine terms as

\(\begin{aligned} I(x,k,u,v)=I_c(x,k,u,v)\sin x-I_s(x,k,u,v)\cos x\,, \end{aligned}\)


\(\begin{aligned} I_{c/s}(x,k,u,v)\equiv \int _{0}^{x}d\tilde{x} \, \tilde{x}\,\left\{ \begin{array}{c} \cos \tilde{x}\\ \sin \tilde{x} \end{array} \right\} f(\tilde{x},k,u,v)\,. \end{aligned}\)

This allows us to take the oscillation average of the spectral density, which is given by

\(\begin{aligned} \overline{I^2(x\rightarrow \infty , k, u, v)} \simeq \frac{1}{2}\left( I_{c,\infty }^2(k,u,v)+I_{s,\infty }^2(k,u,v)\right) \,, \end{aligned}\)

where we took the limit \(x\rightarrow \infty\) as we are interested in GW frequencies which are well within the horizon. It is interesting to note that at this point the main difference between the isocurvature and adiabatic initial condition cases is the different behavior of the transfer functions. The integrals \(I_s\) and \(I_c\) (46) can be analytically carried out as they are integrals of trigonometric functions. However, the expressions are considerably long and so we will only present them in simplified situations. In the next sections we show in more detail the different GW spectrum from isocurvature and adiabatic initial conditions.

Before moving on, we write down the explicit relation between \(\Omega _{\text {GW,c}}\) (41) and the spectral density of GWs measured today. A straightforward relation can be found in the case when the radiation domination era where \(\Omega _{\text {GW,c}}\) is computed corresponds to the standard radiation dominated stage, i.e., prior to BBN and DM domination. In that case, taking into account the redshifting of the GW energy density one finds that [49]

\(\begin{aligned} \Omega _{\text {GW,0}}h^2\approx 1.62\times 10^{-5}\,\Omega _{\text {GW,c}}(k)\,, \end{aligned}\)

where \(\Omega _{\text {GW,0}}\) is the energy density fraction of GWs evaluated today and \(h=H_0/(100 \text {km/s/Mpc})\). \(H_0\) is the Hubble parameter evaluated today. For simplicity, we assumed the standard model of particle physics. In the general case there is a dependence on the effective degrees of freedom, see, e.g., [143, 144].

3 Gravitational waves from DM isocurvature

Let us consider that the matter component with isocurvature fluctuations is the DM. This means that the initial fraction of DM \(\beta\) is fixed its the current abundance, which according to Planck [4] is \(\Omega _{\text {DM,0}}h^2\approx 0.12\). For our practical purposes though, we just need to use that the time of matter-radiation equality, normalized to today, is at \(a^{-1}_{\text {eq}}\approx 3400\) and that the comoving size of the universe at that time was \(k_{\text {eq}}\approx 0.01 \,{\text {Mpc}}^{-1}\). To grasp the magnitude suppression factor \(k_{\text {eq}}/k\), let us write down the GW frequency in terms the relevant wavenumbers, that is

\(\begin{aligned} f_{\text {GW}}\sim 2\times \frac{k}{2\pi }\approx 3\,{\text {Hz}}\left( \frac{k}{10^{15}\,{\text {Mpc}}^{-1}}\right) \,, \end{aligned}\)

where the first factor 2 comes from the fact that two scalar modes source one induced tensor mode. Thus, for the frequencies of interest, say between \(\text {nHz}\) and \(\text {kHz}\), we have suppression factors respectively between \(10^{-8}\) and \(10^{-19}\). This means that initial isocurvature has to be roughly of the inverse order of magnitude so that isocurvature induced GWs are detectable. It should be noted, though, that such large values of initial isocurvature are compatible with cosmological perturbation theory up to the time \(\tau _{\text {NL}}\). For a detailed discussion on the validity of perturbations, see the appendix of Ref. [53]. The current challenge is then not the validity of perturbations but to find a model with such large initial isocurvature. At the end of this section we will discuss some interesting cases where isocurvature need not be that large.

Regarding isocurvature induced GWs, we are mostly interest in the small scale power spectrum. This constitutes fluctuations on scales which entered the horizon much before matter-radiation equality. Therefore, it is a very good approximation to use our solution Eq. (30) during the radiation dominated epoch. Plugging in Eq. (30) into the integrals (46), we find that

\(\begin{aligned} I_{c,\infty }(k,u,v)&=\frac{9}{32u^4v^4\kappa ^{2}}\Bigg \{-3u^2v^2+\left( -3+u^2\right) \left( -3+u^2+2v^2\right) \ln \left| 1-\frac{u^2}{{3}}\right| \nonumber \\&+\left( -3+v^2\right) \left( -3+v^2+2u^2\right) \ln \left| 1-\frac{v^2}{{3}}\right| \nonumber \\&-\frac{1}{2}\left( -3+v^2+u^2\right) ^2\ln \left[ \left| 1-\frac{(u+v)^2}{{3}}\right| \left| 1-\frac{(u-v)^2}{{3}}\right| \right] \Bigg \}\,, \end{aligned}\)


\(\begin{aligned} I_{s,\infty }(k,u,v)&=\frac{9\pi }{32u^4v^4\kappa ^{2}}\Bigg \{9-6v^2-6u^2+2u^2v^2+\left( 3-u^2\right) \left( -3+u^2+2v^2\right) \Theta \left( 1-\frac{u}{\sqrt{3}}\right) \nonumber \\&+\left( 3-v^2\right) \left( -3+v^2+2u^2\right) \Theta \left( 1-\frac{v}{\sqrt{3}}\right) \nonumber \\&+\frac{1}{2}\left( -3+v^2+u^2\right) ^2\left[ \Theta \left( 1-\frac{u+v}{\sqrt{3}}\right) +\Theta \left( 1+\frac{u-v}{\sqrt{3}}\right) \right] \Bigg \}\,, \end{aligned}\)

where we took the limit \(x\rightarrow \infty\) for analytical simplicity. With the kernels Eqs. (50) and (51), we are ready to compute the isocurvature induced GW spectrum via Eq. (41). Let us emphasize that Eqs. (50) and (51) are valid for any primordial spectrum of isocurvature fluctuations. The only requirement is that the relevant fluctuations enter the horizon much before matter-radiation equality. The averaged kernel for GWs induced by adiabatic fluctuations can be found in Appendix 2. For illustration, we also compare the source term (39) between initially adiabatic and isocurvature fluctuations in Fig. 2.

Fig. 2
figure 2

Source term (39) for induced GWs for adiabatic and isocurvature initial conditions, respectively in purple and orange, in terms of \(x=k\tau\). For simplicity we chose \(u=v=1\) but the main behavior is independent on such particular values. Here we consider \(\tau \ll \tau _{\text {eq}}\) and, therefore, the source term only decays after a given mode enters the horizon during radiation domination. See how for adiabatic initial conditions the source term is initial constant and then decays with large amplitude oscillations. Instead for isocurvature initial conditions the source term initially grows and then decays with smaller amplitude oscillations. Note that the frequency of the oscillations is the same in both cases but are out of phase

It is interesting to note that the explicit \(\kappa\) dependence in the kernel (50) and (51) can be absorbed via the definition of an effective power spectrum of curvature perturbations, namely

\(\begin{aligned} {\mathcal {P}}_\Phi ^{\text {eff}}(k)\equiv \kappa ^{-2}{\mathcal {P}}_{\mathcal {S}}(k)\,. \end{aligned}\)

Although this is a mere redefinition, it illustrates the differences with the adiabatic initial conditions. For instance, after using \({\mathcal {P}}_\Phi ^{\text {eff}}(k)\) the averaged kernel resembles that of the adiabatic initial conditions case, in the sense that they have the same resonance and same scaling in the limits of integration [53]. The detailed shape is of course different. Also note that the GW spectrum has a global suppression proportional to \(\kappa ^{-4}\), consistent with our expectations that \(\Phi \propto \rho _m/\rho _r \times S\propto \kappa ^{-1} S\) for isocurvature initial conditions when evaluated at horizon crossing. We proceed to discuss some concrete examples.

3.1 GWs induced by a peaked primordial DM isocurvature spectrum

To have an idea of the shape of the isocurvature induced GW spectrum and its difference with the standard (curvature) induced GWs, let us consider that the primordial spectrum of DM isocurvature fluctuations peaks at given scale, say \(k_p\). For concreteness and to follow previous works we consider a log-normal peak given by

\(\begin{aligned} \mathcal {P}_{S}(k)=\frac{\mathcal {A_S}}{\sqrt{2\pi }\Delta }\exp \left[ -\frac{\ln ^2(k/k_p)}{2\Delta ^2}\right] , \end{aligned}\)

where \(A_S\) is the amplitude and \(\Delta\) is the logarithmic width of the peak. This is standard practice when dealing with GWs induced by primordial adiabatic fluctuations (see, e.g., [145]) and it is well motivated from inflationary models [25, 25, 146,147,148,149,150,151,152,153,154,155,156,157,158,159,160,161,162,163,164,165,166,167,168,169,170]. Similar conclusions extend to the isocurvature case [171, 172]. We note that the formal limit of \(\Delta \rightarrow 0\) corresponds to a Dirac delta peak, namely

\(\begin{aligned} \mathcal {P}_{S}(k)=\mathcal {A_S}\delta (\ln (k/k_p))\,. \end{aligned}\)

A log-normal is considered to be “sharp” if \(\Delta \lesssim 0.1\) and broad otherwise [145]. We take a similar ansatz for primordial curvature fluctuations \(\Phi _i\).

The Dirac delta case is convenient as it allows for an analytical expression for \(\Omega _{\text {GW,c}}\) simply by replacing \(u=v=k_p/k\) in Eqs. (50), (51), and (41). Explicitly, we obtain

\(\begin{aligned} \Omega _{\text {GW,c}}(k)&= \frac{{\mathcal {A}}_S^2}{3} \left( \frac{k}{k_p}\right) ^{-2} \left( 1-\frac{k^2}{4k_p^2}\right) ^2 \left( I^2_{c,\infty }\left( \frac{k_p}{k},\frac{k_p}{k}\right) + I^2_{s,\infty }\left( \frac{k_p}{k},\frac{k_p}{k}\right) \right) \Theta \left( 2 k_p-k\right) , \end{aligned}\)

where the sharp cut-off at \(k=2k_p\) comes from momentum conservation. From Eq. (55), we see that, as in the initially adiabatic case (see, e.g., the discussions in Ref. [49]), the induced GW spectrum has a resonant peak at \(k=2c_sk_p\) and decays as \(k^2\ln ^2k\) in the low frequency tail, namely for \(k\ll k_p\) [173,174,175]. However, contrary to the adiabatic case, there is no destructive interference with vanishing GW spectrum at \(k=\sqrt{2}c_sk_p\). This is because adiabatic initial conditions lead to coherent, order unity, oscillations of \(\Phi\) while for isocurvature initial conditions the oscillations of \(\Phi\) are modulations with a different phase. Thus, one can in principle distinguish the adiabatic and isocurvature cases by the shape of the GW spectrum around \(k\sim k_p\). One of the important differences, though, is that the isocurvature induced GW spectrum has a suppressed amplitude, namely the peak amplitude is given by

\(\begin{aligned} \Omega ^{\text {iso,peak}}_{\text {GW,c}}&\propto {{\mathcal {A}}_S^2}\times \left( \frac{k_{\text {eq}}}{k_p}\right) ^4\,. \end{aligned}\)

This means that the amplitude of the initial power spectrum of isocurvature fluctuations must be very large to partly compensate for the suppression, roughly \({{\mathcal {A}}_S}\propto \left( {k_{p}}/{k_{\text {eq}}}\right) ^2\gg 1\). Such large values of initial isocurvature fluctuations are not in contradiction with the validity of cosmological perturbation. The reason is that any effect of isocurvature is accompanied by a factor \(\rho _m/\rho _r\) during the radiation dominated era. The product of initial isocurvature times \(\rho _m/\rho _r\) is always smaller than unity during the relevant times for our calculations.

For the log-normal peak, we recover similar results as the Dirac delta for \(\Delta \lesssim 0.01\), except that for scales \(k<\Delta \, k_p\), where the slope transitions to a \(k^3\) scaling, as in the adiabatic case [145]. For \(\Delta \gtrsim 0.01\), similar conclusions apply but because of the \(\kappa ^{-2}\) dependence in \({\mathcal {P}}_\Phi ^{\text {eff}}\) (52), coming from the fact that modes that enter earlier or more suppressed, the peak of the GW spectrum moves to lower values of k. We show the numerical results for the Dirac delta and log-normal peak in Fig. 3.

Fig. 3
figure 3

Induced GW spectrum from a log-normal primordial spectrum in terms of comoving wavenumber k normalized to the peak wavenumber \(k_p\). The amplitude of the GW spectrum is normalized such that \({\mathcal {A}}_{S}=1\) and \({\mathcal {A}}_{\Phi }=1\) in Eq. (53) for \(\mathcal {P}_{S}\) as well as \(\mathcal {P}_{\Phi }\), i.e., primordial isocurvature and curvature respectively. On the left we respectively show the induced GWs from primordial curvature and primordial isocurvature in purple and orange, in the case of a Dirac delta power spectrum. See how the overall shape is similar, that is a resonant peak at \(k=2c_sk_p\) and a \(k^2\) low frequency tail with a logarithmic running. However, the shape around the peak and the dip is different which breaks the degeneracy of the GW signals. On the right we show the isocurvature induced GW spectrum from a log-normal peak with \(\Delta =0.01,0.3,1\) respectively in red, green and blue. The orange dotted line show the Dirac delta result of the left figure. See how the peak of the induced GW spectrum shifts to lower k for larger \(\Delta\) as explained by Eq. (52)

3.2 Remarks, issues, and future work

So far we have discussed that a large amplitude of primordial isocurvature fluctuations is compatible with perturbation theory, as long as they satisfy Eq. (33). We also assumed that the initial isocurvature fluctuations can be considered as Gaussian. Such assumption, while expected to give a rough approximation, is strictly speaking not correct. The reason is that isocurvature fluctuations are mostly density fluctuations of the matter fluid, i.e., \(S\propto \delta \rho _m/\rho _m\), and while \({\mathcal {P}}_S\gg 1\) is certainly possible, the probability distribution of S cannot be Gaussian as \(S>-1\) (i.e., \(\rho _m+\delta \rho _m>0\)). Thus, a correct distribution function would be highly skewed, starting at \(S=-1\) and with large tails for large S, so as to keep the average \(\langle S \rangle =0\). The main issue here is that in the absence of a concrete realization it is difficult to parametrize the distribution function as it cannot be expanded in terms of Gaussian distributions. Some estimates are given in the appendix of Ref. [53] but the general expectation is that large tails lead to large 4-point functions and a larger amplitude of the induced GWs. Thus, what we computed in the previous section might well be a lower bound on the amplitude of DM isocurvature induced GWs. It would be interesting to study a concrete case and confirm these expectations.

It should be noted that, although we focused on the case of cold DM as an exact pressure-less fluid, our discussion may be extended to the case of isocurvature from scalar dark matter, such as ultra-light dark matter[136, 137, 176]. The reason is that most of the isocurvature induced GWs are generated when the relevant curvature perturbation enters the sound horizon (see, e.g., Fig. 2 and the discussion around Eq. (30)). But, any difference with respect to the case of cold DM is expected to occur only for small scale DM density fluctuations once they are well inside the Hubble horizon. Thus, our results should provide a good approximation to the case of light scalar field dark matter.Footnote 4 Furthermore, we expect that similar results apply to the case of isocurvature from vector dark matter, as discussed in Ref. [177]. As future work, it would be interesting to study in more detail induced GWs by isocurvature from ultra-light and vector dark matter.

It is also interesting to consider the hypothetical case where the possible GW background signal reported by the PTA collaborations [69,70,71,72,73,74,75,76,77,78] is due to isocurvature induced GWs. Although we will not carry out a detailed data analysis, we may infer a good order of magnitude estimate for the required amplitude for the primordial spectrum of DM isocurvature fluctuations. From the analysis of GWs induced by primordial adiabatic fluctuations we have that, if given by a Dirac delta spectrum, one needs \({\mathcal {A}}_\Phi \sim 10^{-2}-10^{-1}\) and a peak position around \(f_p\sim 10^{-7}\,\text {Hz}\) [79]. This means that for the primordial DM isocurvature one requires \({\mathcal {A}}_S (k_{\text {eq}}/k_p)^2\sim {\mathcal {A}}_\Phi\), which roughly corresponds to \({\mathcal {A}}_S\sim 10^{18}-10^{19}\). While this value of \({\mathcal {A}}_S\) is certainly big it is still within the validity range of our induced GW spectrum calculations as \({\mathcal {A}}_S (k_{\text {eq}}/k_p)^2\ll 1\). With future detectors in the \(\mu\)Hz window, such as \(\mu\)-Ares [121] it may be possible to see the peak of the GW spectrum and to distinguish the nature of the initial conditions. We present an example in Fig. 4.

Fig. 4
figure 4

Example of induced GWs from a Dirac delta spectrum for primordial curvature and isocurvature fluctuations. As an illustration we chose \(A_\Phi =0.01\) and \(A_S=0.05 \,(k_{p}/k_{\text {eq}})^2\sim 10^{18}\) to fit the PTA data by eye. The peak is chosen at \(f_p\sim 10^{-7}\,\text {Hz}\) (\(k_p\sim 10^7\,{\text {Mpc}}^{-1}\)). Gray violins indicate the recent NANOGrav results [79] and in blue and orange we show power-law integrated sensitivity curves [178] for LISA and Taiji [122, 123]. See [128,129,130, 179] for the sensitivity curves. The horizontal blue and purple dashed lines respectively show the current constraint from BBN [180,181,182] and future constraints from CMB-S4 experiments [181, 183]. We also include future sensitivity of \(\mu\)-Ares [121] and Lunar Laser Ranging from Ref. [184] (see also [185]). GW detectors such as \(\mu\)-Ares may tell apart whether the GWs were induced by adiabatic or isocurvature initial conditions

4 Gravitational waves from the PBH dominated early universe

An interesting model which includes an early stage of matter domination ended by a fast transition to radiation domination is the PBH reheating scenario [186, 187]. In this scenario, tiny PBHs are formed after inflation with a fraction large enough such that they dominate the early universe before evaporating via Hawking radiation. The most interesting part is that, as shown by Refs. [131, 188], a sudden transition from matter to radiation domination greatly enhances the spectrum of induced GWs. In our case, this generates a distinct GW signal of a PBH dominated early universe.

4.1 PBH formation and evaporation

Consider that initially we have a universe filled with a homogeneous radiation fluid. Then, at some moment in time, PBHs form roughly at the same time and with the same mass \(M_{\text {PBH}}\) (this is often called a monochromatic mass function). We assume that PBHs formed by the collapse of large primordial fluctuations (although it does not have to be necessarily the case) and so their mass at formation is related to the Hubble parameter at the time of formation, say \(H_{\text {f}}\), via \(M_{\text {PBH,f}}\sim 4\pi M_{\text {pl}}^2/H_{\text {f}}\). Note that fact that PBHs formed from adiabatic primordial fluctuations does not affect our simplified picture of an initial homogeneous radiation fluid, as we shall see. We also assume from now on that the particle content of the universe after evaporation is given by the standard model of particles particle. Additional particles might change the precise coefficients through the effective degrees of freedom. We refer the interested reader to Refs. [131, 132] for the details.

After formation, PBH make a fraction \(\beta\) of the total energy density at formation. It is interesting to relate the \(\beta\) to the number density of PBHs, namely

\(\begin{aligned} \beta =\frac{\rho _{\text {PBH,f}}}{\rho _{r,\text {f}}}\sim \frac{4\pi }{3H_{\text {f}}^3}n_{\text {PBH,f}}\,, \end{aligned}\)

where we used that \(\rho _{\text {PBH,f}}=M_{\text {PBH,f}}\times n_{\text {PBH,f}}\) and \(\rho _{r, \text {f}}\approx 3H_{\text {f}}^2M_{\text {pl}}^2\). Thus, \(\beta\) can be interpreted as the mean number of PBHs for Hubble volume. We emphasize that it refers to the mean number because in general there will be statistical number density fluctuations.

Now, while the number density of PBHs is conserved, which in an expanding universe means that \(n_{\text {PBH}}\propto a^{-3}\), the PBH mass decreases due to Hawking radiation. One finds that the mass decays as [131]

\(\begin{aligned} M_{\text {PBH}}(t)\approx M_{{\text {PBH}},f}\left( 1-\frac{t}{t_{\text {eva}}}\right) ^{1/3}\,, \end{aligned}\)

where \(t_{\text {eva}}\) is the time of evaporation and is approximately given by

\(\begin{aligned} t_{\text {eva}}\sim \frac{M_{\text {PBH,f}}^3}{M_{\text {pl}}^4}\sim 0.4 \,{\text {fs}}\times \left( \frac{M_{\text {PBH,f}}}{10^4\,{\text {g}}}\right) ^{3}\,, \end{aligned}\)

where \({\text {fs}}=10^{-15}\,{\text {s}}\) are femtoseconds. Here t is cosmic time related to conformal time by \(dt=ad\tau\). Since we assumed that all PBHs have the same mass and that formed roughly at the same time, they all evaporate at the same time as well. From the Hubble parameter at the time of evaporation (\(H_{\text {eva}}\sim 1/t_{\text {eva}}\)) we can estimate the temperature of the radiation fluid filling the universe right after evaporation. In doing so, we find that

\(\begin{aligned} T_{\text {eva}}\approx 2.76\times 10^4\,{\text {GeV}}\left( \frac{M_{{\text {PBH}},f}}{10^4{\text {g}}}\right) ^{-3/2}\,. \end{aligned}\)

From the Hubble parameter we can also compute the comoving size of the Hubble horizon at evaporation, namely [132]

\(\begin{aligned} k_{\text {eva}}={\mathcal {H}}_{\text {eva}}\approx 4.7\times 10^{11}{\text {Mpc}}^{-1}\left( \frac{M_{{\text {PBH}},f}}{10^4{\text {g}}}\right) ^{-3/2}\,. \end{aligned}\)

As expected, the larger the PBH mass the later the reheating and the larger the comoving horizon at evaporation. It should be noted that PBHs must reheat the universe well before BBN, which imposes that \(T_{\text {eva}}>4\,{\text {MeV}}\) and translates into \(M_{{\text {PBH}},f}<5\times 10^8\,{\text {g}}\). We can also compute the minimum value of \(\beta\) such that PBHs eventually dominate the universe before evaporating, which is given by

\(\begin{aligned} \beta >6\times 10^{-10}\left( \frac{M_{{\text {PBH}},f}}{10^4{\text {g}}}\right) ^{-1}\,. \end{aligned}\)

When dealing with fluctuations, we treat the collection of PBHs as an almost pressureless matter fluid. The evaporation is then considered as the decay of “matter” into radiation via a decay rate given by

\(\begin{aligned} \Gamma \equiv -\frac{d\ln M_{\text {PBH}}}{dt}\,. \end{aligned}\)

Then, the energy density of the PBH fluid is gradually transferred to radiation and obeys

\(\begin{aligned} \dot{\rho }_{\text {PBH}}+\left( 3H+\Gamma \right) \rho _{\text {PBH}}=0\,, \end{aligned}\)

where \(\dot{\,}\equiv d/dt\). The equation for the radiation energy density is similar but with the opposite sign in front of \(\Gamma\), as to satisfy energy conservation. One can check that the total evaporation of PBHs is not instantaneous but takes about a quarter of an e-fold. We show in Fig. 5 the evolution of \(\rho _{\text {PBH}}\) and \(\rho _r\) for a particular example.

Fig. 5
figure 5

Energy density fraction \(\Omega _y=\rho _y/(3H^2M_{\text {pl}}^2)\), with \(y={m,r}\) respectively for PBH (matter) in blue and radiation in red, in terms of the scale factor normalized at formation. We chose the initial PBH fraction to be \(\beta =10^{-4}\) for illustrative purposes. See how PBH domination can last for several e-folds and that evaporation is faster than one e-fold

The main point that we want to emphasize here is that fluctuations with \(k\gg \Gamma\) have a typical time scale much larger than the evaporation rate and, as such, they are very much affected by the finite duration of evaporation [131]. This translates into a suppression factor for curvature fluctuations with \(k\gg \Gamma\) after evaporation given by

\(\begin{aligned} {\mathcal {S}}_\Phi (k)\equiv \frac{\Phi _{\text {lRD}}}{\Phi ^{\text {instant}}_{\text {lRD}}}\approx \left( \frac{k}{k_{\text {eva}}}\right) ^{-1/3} \,, \end{aligned}\)

where “lRD” means late radiation domination and “instant” refers to the amplitude of \(\Phi\) obtained by an instant matching from matter to radiation domination. Note that the exponent 1/3 in Eq. (65) is directly related to the exponent in the decay of the PBH mass (58). The relation between \(\Phi\) and \(\rho _{\text {PBH}}\) follows from the Poisson equation on very subhorizon scales, that is \(k^2\Phi \sim a^2\rho _{\text {PBH}}\delta _{\text {PBH}}\). We must take the suppression factor \({\mathcal {S}}_\Phi (k)\) (65) into account when computing any transfer function for the curvature perturbation. We suggest to read Ref. [131] for more details on the calculations.

4.2 PBH number density fluctuations and initial isocurvature

While on average PBH are homogeneously distributed, the fact that PBHs are discrete objects introduces inhomogeneities. For instance, the mean PBH separation is given by [66]

\(\begin{aligned} d_{\text {f}}\equiv \left( \frac{3}{4\pi n_{{\text {PBH}},{\text {f}}}}\right) ^{1/3}\approx \beta ^{-1/3}{H_{\text {f}}}^{-1}\,. \end{aligned}\)

On scales larger than \(d_{\text {f}}\), we have a coarse grained fluid picture but on scales smaller than \(d_{\text {f}}\) we either see a PBH or not. Now, since PBH formation is a rare event (the number of PBH per Hubble volume at formation \(\beta\) is very small), we can approximate the process of formation as PBHs appearing randomly and uniformly distributed in space (in other words, everywhere has the same probability of hosting a PBH). And, this means that the probability of having n PBHs in a given volume is described by Poisson statistics. This allows us to compute the variance of number density fluctuations, which for Poisson fluctuations is constant in wavenumber and given by [66]

\(\begin{aligned} \left\langle \frac{\delta n_{\text {PBH,f}}(k)}{n_{\text {PBH,f}}}\frac{\delta n_{\text {PBH,f}}(k')}{n_{\text {PBH,f}}}\right\rangle =\frac{4\pi }{3}\frac{d_{\text {f}}^3}{a_{\text {f}}^3}\,\delta (k+k')\,, \end{aligned}\)

where the factor \(a_{\text {f}}\) comes from the fact that we are working with comoving scales. In terms of dimensionless quantities, a constant physical variance implies a dimensionless variance proportional to \(k^3\). Thus, the variance of fluctuations increases with wavenumber until the fluid picture is no longer valid, which occurs at an “ultra-violet” cut-off given by the mean inter-PBH comoving separation, namely

\(\begin{aligned} k_{\text {UV}}=a_{\text {f}}/d_{\text {f}}\approx \beta ^{1/3}{{\mathcal {H}}_{\text {f}}}=\beta ^{1/3}{{k}_{\text {f}}}\,, \end{aligned}\)

where we used that \(k_{\text {f}}={\mathcal {H}}_{\text {f}}\).

Such initial PBH number density fluctuations are in fact initial isocurvature fluctuations. Simply put, by energy conservation any missing part of the radiation fluid (that ended up in a PBH) is compensated by the PBHs themselves, such that at formation, we have \(\delta \rho _{\text {PBH,f}}+\delta \rho _{r,\text {f}}=0\). This implies that

\(\begin{aligned} S_f=\frac{\delta \rho _{\text {PBH,f}}}{\rho _{\text {PBH,f}}}-\frac{3}{4}\frac{\delta \rho _{r,\text {f}}}{\rho _{r,\text {f}}}\approx \frac{\delta \rho _{{\text {PBH}},f}}{\rho _{{\text {PBH}},f}}\approx \frac{\delta n_{{\text {PBH}},f}}{n_{{\text {PBH}},f}}\,, \end{aligned}\)

where we used that initially \(\rho _{\text {PBH}}\ll \rho _r\). From Eq. (67), the dimensionless initial isocurvature power spectrum then reads

\(\begin{aligned} {\mathcal {P}}_{S}(k)=\frac{2}{3\pi }\left( \frac{k}{k_{\text {UV}}}\right) ^3\,. \end{aligned}\)

This is the initial isocurvature that will eventually generate induced GWs [66, 132].

4.3 GWs from PBH isocurvature fluctuations after PBH evaporation

PBH number density fluctuations provide isocurvature initial conditions for cosmological fluctuations. However, as we discussed, induced GWs are mainly sourced by curvature fluctuations. Thus, before computing the induced GWs, we have to compute the transfer functions for \(\Phi\). As we will be interested in the GWs sourced by curvature fluctuations after evaporation and on very small scales (since \(k_{\text {UV}}\gg k_{\text {eva}}\)), we just have to follow \(\Phi\) until the end of the PBH dominated era. There are of course GWs induced during the early radiation domination and the early matter domination phase. However, they turn out to be subdominant and, therefore, we neglect them. See Ref. [66] for the GWs induced during the PBH dominated phase.

To compute the transfer functions, we start noting that the comoving Hubble parameter at equality is proportional to \(\beta\), that is

\(\begin{aligned} {k_{\text {eq}}}={\sqrt{2}\beta }\,k_{\text {f}}\,, \end{aligned}\)

where \(k_{\text {eq}}={\mathcal {H}}_{\text {eq}}\). It then follows that we have the resulting hierarchy: \(k_{\text {f}}\gg k_{\text {UV}}\gg k_{\text {eq}}\gg k_{\text {eva}}\). To put it into words, the initial fluctuations on scales similar to \(k_{\text {UV}}\) are initially superhorizon (they are on scales larger than \(k_{\text {f}}\)). But, these scales around \(k_{\text {UV}}\) enter the horizon well before matter-radiation equality and so \(k_{\text {UV}}\gg k_{\text {eq}}\). Then at evaporation all fluctuations of interest are largely subhorizon. This means that we can use the second line of the transfer function in Eq. (34), namely the transfer function for isocurvature initial conditions, together with the suppression factor (65), which yields

\(\begin{aligned} T_{\Phi }(t_{\text {eva}}; k\gg k_{\text {eq}})=T_{\Phi _{\text {iso}}}(k){\mathcal {S}}_\Phi (k)=\frac{3}{4}\left( \frac{k}{k_{\text {eq}}}\right) ^{-2}\left( \frac{k}{k_{\text {eva}}}\right) ^{-1/3}\,. \end{aligned}\)

To compute the GWs induced after evaporation, we have to continue the constant value of \(\Phi\) during the PBH (matter) domination to the late radiation domination. After matching, we find that

\(\begin{aligned} \Phi _{\text {lRD}}(k\tau )=\frac{1}{c_sk\bar{\tau }}\left( C_1j_1(c_sk\bar{\tau })+C_2y_1(c_sk\bar{\tau })\right) \quad ; \quad \bar{\tau }\equiv \tau -\tau _{\text {eva}}/2, \end{aligned}\)


\(\begin{aligned} C_1=-\Phi _{\text {eva}}(k)\left( c_sk\tau _{\text {eva}}/{2}\right) ^{2}\cos (c_sk\tau _{\text {eva}}/2)\quad , \quad C_2=C_1 \tan (c_sk\tau _{\text {eva}}/2)\,. \end{aligned}\)

Now, with the evolution of \(\Phi\) we can re-compute the kernel for the induced GWs. Here, we only provide the main steps of the calculations and we refer the reader to Ref. [131, 132] for the details.

The most important point is to realize that, while the amplitude of \(\Phi\) in Eq. (73) starts from a constant and quickly decays, the amplitude of its time derivative, \(\Phi '\), begins at zero (by continuity) and suddenly jumps to an amplitude proportional to \(k/k_{\text {eva}}\), which is huge for \(k\sim k_{\text {UV}}\). The reason for this jump is that during the matter dominated era \(\Phi\) remained constant but the relevant scales became more and more subhorizon. Then, at the late radiation domination, \(\Phi\) resumes its decay but with very fast oscillations. Yet, since it has not decayed during the matter dominated phase, the amplitude of \(\Phi '\) is suddenly very large by a factor \(a_{\text {eva}}/a_{\text {UV}}\sim k_{\text {UV}}/k_{\text {eva}}\) (the supposed decay if there never were a matter dominated era).

Another explanation for the enhancement is given in Ref. [188]. During the matter dominated phase PBH density fluctuations grow, as standard DM fluctuations do. But suddenly all those density fluctuations are converted into radiation which wants to propagate. The large density fluctuations become density waves with a huge velocity (note that by Einstein Equations in Appendix  1 \(V\propto \Phi '/{\mathcal {H}}\)). Inspecting the source term (39), we see that the main contribution to the kernel (43) comes from \(\Phi '\) and, basically, we can approximate the kernel by

\(\begin{aligned} I_{\text {lRD}}(x,u,v,x_{\text {eva}})\approx \frac{\bar{x}}{2}uv\int _{x_{\text {eva}}/2}^{\bar{x}} d{ {\bar{x}}_1}\,{ {\bar{x}}_1}^2\, G^{\text {lRD}}(\bar{x},{ {\bar{x}}_1})\frac{dT_\Phi ^{\text {lRD}}(u{ {\bar{x}}_1})}{d(u{ {\bar{x}}_1})}\frac{dT_\Phi ^{\text {lRD}}(v{ {\bar{x}}_1})}{d(v{ {\bar{x}}_1})}\,, \end{aligned}\)

where \(G^{\text {lRD}}(\bar{x},{ {\bar{x}}_1})\) is given by (44). After integration and selecting only those terms related to the resonant peak (those at \(c_s(u+v)\sim 1\) where the frequency of the source term is equal to that of the tensor mode), we arrive at

\(\begin{aligned} \overline{I^2_{\text {lRD,res}}}&(u+v\sim c_s^{-1},x_{\text {eva}})\approx \frac{c_s^4u^2v^2}{2^{15}}x_{\text {eva}}^8T_{\Phi }^2(t_{\text {eva}}; vk)T_{\Phi }^2(t_{\text {eva}}; uk){\text {Ci}}^2(|1-(u+v)c_s|x_{\text {eva}}/2)\,. \end{aligned}\)

To compute the induced GW spectrum, we plug in the averaged kernel (76) into Eq. (41) and integrate. However, for sufficiently peaked spectrum the integral is well approximated by a power-law with a sharp cut-off [131,132,133]. A power spectrum is considered to be peaked enough if the integrand in Eq. (41) grows for growing u and v. Since it contains \((u v)^2\sim v^4\), the effective spectral index of \(P_\Phi (k)\) should be larger that \(-5/2\). Otherwise the integral has to be carried out numerically. Our effective power spectrum of curvature fluctuations at evaporation reads

\(\begin{aligned} {\mathcal {P}}_\Phi ^{\text {eva}}(k\gg k_{\text {eq}})\approx T^2_{\Phi _{\text {iso}}}(k)\,{\mathcal {S}}^2_\Phi (k)\,{\mathcal {P}}_{S}(k)=\frac{3}{8\pi }\left( \frac{k}{k_{\text {eq}}}\right) ^{-4}\left( \frac{k}{k_{\text {eva}}}\right) ^{-2/3}\left( \frac{k}{k_{\text {UV}}}\right) ^{3}\,, \end{aligned}\)

which has an effective spectral index of \(-5/3\). We see that our effective power spectrum is well within the power-law approximation for the integral.

With the aforementioned approximations, the PBH isocurvature induced GW spectrum is roughly given by

\(\begin{aligned} \Omega _{\text {GW,c}}\approx \Omega _{\text {GW,res}}\left( \frac{k}{k_{\text {UV}}}\right) ^{11/3}\Theta (k_{\text {UV}}-k)\,, \end{aligned}\)

where the peak amplitude of the GW spectrum is approximately given by

\(\begin{aligned} \Omega _{\text {GW,res}}(k\sim k_{\text {UV}})&\approx \frac{1}{24576\pi \,2^{1/3}\sqrt{3}}\left( \frac{k_{\text {UV}}}{k_{\text {eva}}}\right) ^{17/3} \left( \frac{k_{\text {eq}}}{k_{\text {UV}}}\right) ^8 \approx 10^{30}\beta ^{16/3} \left( \frac{M_{{\text {PBH}},f}}{10^4{\text {g}}}\right) ^{34/9}\,. \end{aligned}\)

Furthermore, using Eqs. (68) and (61), we find that the peak GW frequency is located at

\(\begin{aligned} f_{\text {UV}} \approx 1700\, {\text {Hz}}\left( \frac{M_{\text {PBH,f}}}{10^4\,{\text {g}}}\right) ^{-5/6}\,. \end{aligned}\)

Thus, on one hand, we see that for \(5\times 10^8{\text {g}}>M_{\text {PBH,f}}>10^4\,{\text {g}}\) the peak of the GW spectrum enters in the observable range of LIGO/VIRGO/KAGRA and future detectors such as ET and CE. It is interesting to note that the peak frequency only depends on the PBH mass at formation providing a clear probe of the PBH mass. On the other hand, the peak amplitude of the GW spectrum can be used to probe the initial PBH fraction. For instance, requiring that the peak amplitude is not in contradiction with current BBN constraints imposes that

\(\begin{aligned} \beta < 10^{-6}\left( \frac{M_{\text {PBH,f}}}{10^4\,{\text {g}}}\right) ^{-17/24}\,. \end{aligned}\)

This is the best constraint we have on the fraction of PBHs that evaporated before BBN. We show the resulting bounds on \(\beta\) and an example of the GW spectrum from the PBH dominated universe in Fig. 6.

Fig. 6
figure 6

On the left we show in the shaded orange region the parameter space for the initial fraction of PBHs, \(\beta\), such that PBHs dominate the early universe (\(\beta >\beta _{\text {min}}\); magenta line) and that the induced GWs satisfies BBN bounds (\(\beta <\beta _{\text {max}}\); purple line), as a function of PBH mass. On the right we present an example of induced GWs from the PBH dominated early universe with \(M_{\text {PBH,f}}=10^6\,{\text {g}}\) and \(\beta =2\times 10^{-8}\). We also show the power-law integrated sensitivity curves [178] for DECIGO, Einstein Telescope (ET), Cosmic Explorer (CE), Voyager, and LIGO A+ experiments [128,129,130, 179] as well as the upper bounds from the LIGO/Virgo/KAGRA collaboration [120]. The horizontal lines show BBN [180,181,182] and CMB-S4 experiments constraints, respectively in blue and purple [181, 183]. Finding a very peaked GW signal with a slope of \(k^{11/3}\) above Hz frequencies might be an indication of the PBH reheating scenario

4.4 Remarks, issues, and future work

In our calculations, we have mainly focus on the curvature perturbation \(\Phi\) right after evaporation and we trusted our linear perturbation theory since \(\Phi \ll 1\). However, when we look into the evolution of matter density fluctuations, which on small scales are related to \(\Phi\) by the Poisson equation \(2k^2\Phi =a^2\delta \rho\), we find that \(\delta \rho /\rho\) at evaporation is larger than unity [66, 132]. Depending on the PBH mass and fraction, the amplitude of density fluctuations might be of a few orders of magnitude. This signals the on-set of non-linear physics, such as halo formation and PBH mergers. Unfortunately, one would need numerical simulations to explore such non-linear regimes. And, although it is an interesting question, it is out of the scope of this review. Note that we could also be conservative and stop our calculations of induced GWs when density fluctuations enter the non-linear regime. But, the fact is that we do not expect the GW production to stop by non-linearities. Thus, we take our result, Eq. (79), as a rough order of magnitude estimate. For an interesting proposal in terms of turbulences in the radiation fluid after evaporation, see Ref. [189].

In this work, we have also neglected the contribution for the adiabatic fluctuations that formed the PBHs, which for large PBH masses might enter the sensitivity of ET in the low frequency tail, and any primordial adiabatic fluctuations enhanced by the PBH dominated stage [131]. For works including both PBH isocurvature and adiabatic induced GWs, see Refs. [190,191,192]. For articles investigating the effect of a PBH dominated stage on the GW spectrum from cosmic strings, see Refs. [193,194,195]. Other applications of the GWs from the PBH dominated universe can be found in Refs. [196,197,198,199,200] in the context of modified gravity. Another interesting signal of the PBH dominated universe are GWs from Hawking evaporation itself, which could be important for highly spinning PBHs [133, 181, 201,202,203,204,205]. For a similar application in the case where inflaton oscillons dominate the universe after inflation see Ref. [206].

We end this section by emphasizing that the formalism employed is also applicable to any early matter dominated stage followed by a transition to radiation domination, as discussed in [133] (although see Ref. [207] for a transition to kinetic domination). For instance, if at the start of the early matter domination we have a spectrum of curvature fluctuations given by a power-law, for instance

\(\begin{aligned} {\mathcal {P}}_{\Phi }={\mathcal {A}}_\Phi \left( \frac{k}{k_{\text {UV}}}\right) ^{-n}\Theta (k_{\text {UV}}-k)\,, \end{aligned}\)

with some arbitrary cut-off at \(k_{\text {UV}}\), then the resulting induced GW spectrum after evaporation reads [133]

\(\begin{aligned} \Omega _{\text {GWs,rh}}(k\sim k_{\text {UV}})\approx \Omega ^{\text {peak}}_{\text {GWs}}\left( \frac{k}{k_{\text {UV}}}\right) ^{7-2n}\Theta (k_{\text {UV}}-k)\,, \end{aligned}\)


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

It is also important to note that if transition is rather gradual, the enhancement of the induced GW spectrum disappears [208, 209]. This is particularly important if PBHs have a relatively broad mass function [131]. We note, however, that this does not exclude the possibility of forming larger mass PBHs, such as those that could explain the DM [23,24,25,26,27,28,29], some of the LIGO events [30,31,32,33] and the seeds of supermassive BHs [34, 35]. In that case, the PBH mass function would present different peaks with at least one sharp peak for very light PBHs and another peak at the other scale of interest.

5 Discussions and conclusions

Initial isocurvature fluctuations may source induced GWs, in addition to any initial adiabatic fluctuations. One possibility is that isocurvature fluctuations do so via an induced curvature fluctuation, which is for example the case of initial DM isocurvature. The revelant GWs are then induced during the standard radiation dominated era. In that situation, the amplitude of the induced GWs is roughly suppressed by a factor \(\rho _m/\rho _r\) evaluated at horizon crossing for a given k mode. This results in a large amplitude of initial isocurvature fluctuations, if such GWs enter the observable window of current and future GW detectors. Although a large amplitude (\(S\gg 1\)) is consistent with cosmological perturbation theory, it has the caveat that isocurvature fluctuations must be extremely non-Gaussian [53]. We expect that such non-Gaussianities may in fact enhance the production of induced GWs. As an example, we considered in Fig. 4 the possibility that DM isocurvature fluctuations explain the tentative GW signal reported by PTAs.

Another possibility is that the field mainly responsible for initial isocurvature fluctuations dominates the very early universe. Initial isocurvature fluctuations are then converted into curvature fluctuations which directly source induced GWs. As an example of this situation we considered the PBH dominated universe and the number density fluctuations coming from their Poisson statistics. When the PBH mass function is monochromatic, the final PBH evaporation is almost instantaneous, resulting in a large production of induced GWs right after evaporation. One may use this signal to probe the PBH reheating scenario and place constraints on the initial fraction of PBHs, as we showed in Fig. 6. Although our estimate of the induced GW spectrum (79) might be susceptible to non-linearities that occur during the PBH dominated phase and to the assumption of a monochromatic mass function, it provides hopes that we may be able to probe the PBH reheating scenario.

As an extension of our result, we may wonder what implications would Planck remnants [210] at the end of evaporation have for the induced GW signal. See refs. [211,212,213,214,215] for reviews and models on black hole remnants. Recently in Ref. [216], it has been shown that if such remnants constitute the totality of the DM then the initial PBH mass must be \(5\times 10^5{\text {g}}\) and the induced GW spectrum must peak at around \(100\,{\text {Hz}}\). Although it is a speculative possibility, such rather precise prediction for the peak frequency of the induced GWs might be one of the unique ways to test the PBH remnant scenario as DM.

We would like to end this review with one interesting consequence of the formalism for early isocurvature induced GWs presented in this review. As argued in Refs. [67, 68], the formation of compact solitonic structures, such as oscillons, monopoles, domain walls, Q-balls, cosmic strings, etcetera, leads to a Poissonian distribution on scales much larger than the mean inter-soliton separation, just in the PBH scenario of Section 4. This means that statistical fluctuations of the number density of compact structures will lead to a dimensionless power spectrum proportional to \(k^3\). Such power spectrum of isocurvature fluctuations first induces GW during the radiation dominated phase, as in Section 3, and might give another contribution if solitons dominate the early universe before decaying. Such generic prediction has been named “Universal Gravitational Waves of Solitons” [67, 68].

We conclude by emphasizing that, while the consequences of primordial adiabatic fluctuations have been extensively studied, initial isocurvature fluctuations present new opportunities to test the physics of the unexplored very early universe.