Cosmological Gravitational Waves from Isocurvature Fluctuations

Gravitational waves induced by large primordial curvature fluctuations may result in a sizable stochastic gravitational wave background. Interestingly, curvature fluctuations are gradually generated by initial isocurvature fluctuations, which in turn induce gravitational waves. Initial isocurvature fluctuations commonly appear in multi-field models of inflation as well as in the formation of scattered compact objects in the very early universe, such as primordial black holes and solitons like oscillons and cosmic strings. Here we provide a review on isocurvature induced gravitational waves and its applications to dark matter and the primordial black hole dominated early universe.


I. 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 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.
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 [51,52]. 3 Isocurvature fluctuations generally occur in multi-field models of inflation [61,62], phase transitions [63] 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 [56,[64][65][66][67].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 CDM isocurvature fluctuations can be found in Ref. [51].
This work is organized as follows.In § II 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 discuss some applications.In § III we view GWs induced by CDM, based on the results of Ref. [52].In § IV we consider the isocurvature due to PBH number density fluctuations in the PBH reheating scenario.Here we will base our discussions on Refs.[65,[129][130][131].Lastly, we discuss more applications in § V. Some details of the calculations can be found in the appendices and in the aforementioned references.We work in reduced Planck units where c = ℏ = 1 and M pl = (8πG) −1/2 = 1.

II. 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 where g µν is the metric, ρ and p respectively are the energy density and pressure and u µ 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 = ρ r /3.For the metric g µν 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 where a is the scale factor, τ is conformal time, Ψ and Φ 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 App. A. In the same appendix we describe the notation for matter perturbations as well.
Most important to us is the fact that energy conservation (∇ µ T µν = 0) at the background level requires that ρ r ∝ a −4 and ρ m ∝ a −3 (also see App.A).From the different dilution of the energy densities, we see that and, therefore, matter fields eventually dominate the energy density of the universe, assuming they do not decay.Thus, if we call β to the initial fraction of matter, namely where "i" refers to the initial time, the universe will be matter dominated after a/a i > β −1 .Actually, there is an exact analytical solution for the scale factor in the radiation-matter universe.
It reads [7,132] a(τ ) where ( √ 2 − 1)τ * = τ eq .The subscript "eq" means matter-radiation equality, i.e. the time when ρ r = ρ m .It is easy to check that Eq. (2.6) goes from the radiation dominated universe where a ∼ τ to the matter dominated universe with a ∼ τ 2 .This is for the moment all we need to understand the source of induced gravitational waves.

A. 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 where G µν is the Einstein tensor and for simplicity we set M pl = (8π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 P ij ab the transverse-traceless projector (which can be found in, e.g., Refs.[48,130]), then we schematically have at linear order where H = a ′ /a, ′ = d/dτ and the superscript (1) refers to linear perturbations.We will use a superscript (2) to refer to second order perturbations.Eq. (2.8) basically tells us that since T µν 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 which after some simplifications the source term is given by [130] (2.10) Note that we selected the scalar component of the perturbation of the spatial velocity in Eq. (2.1), that is we took u i = a∂ i v. Also in Eq. (2.9) we considered the second order expansion of G 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 Ψ + Φ = 0 in the presence of no anisotropic stress (see App. A).
Let us discuss the secondary source to GWs, Eq. (2.10), within the big picture.First, we see that from G (2) ab we obtain gradients of Φ.So, one source of secondary gravitational waves are curvature (metric) fluctuations.Second, from T µν 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 cosmic 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 [133] (2.11) In terms of these variables we find that Eq. (2.10) is given by 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 Φ i = 0) must then be suppressed by a factor ρ m /ρ 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 ρ m /ρ, 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 § III for CDM fluctuations and the second in § IV for the PBH dominated early universe.We now proceed with the evolution of the curvature fluctuations and the formal solutions to induced GWs.

B. 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,134], Here we used ϕ to denote the spatial curvature perturbation which corresponds to Φ in the Newton gauge.In the uniform density slicing where δρ = 0 we have that ζ = −ϕ.The convenience of using ζ is that it is conserved on superhorizon scales.One may also define a curvature fluctuation for each fluid, namely where in our case y = {r, m}.And, from the individual definitions, the notion of isocurvature follows as the relative individual curvature fluctuations, namely where we used that ρ ′ y + 3H(ρ y + p y ) = 0 and that p r = 1/3ρ r .The physical intuition behind such formal definitions, Eqs.(2.13) and (2.15), is the following.Curvature fluctuations are metric fluctuations even in the time slicing where there are no total density fluctuations (δρ = 0).Adiabatic initial conditions then require no initial isocurvature, that is S i = 0, and ζ i ̸ = 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 ζ i = 0 but S i ̸ = 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 App. A), which reads where we imposed that initially Φ ′ 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 Φ i are then proportional to δρ i .Thus, adiabatic initial conditions correspond to S i = 0 and Φ i ∝ δρ i ̸ = 0 and isocurvature initial conditions to Φ i ∝ δρ i = 0 and S i ̸ = 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 ρ has also dominant δρ) while isocurvature initial conditions are mainly given by the sub-dominant field (since δρ m + δρ r = 0 and S depends inversely in ρ m ).
The closed system of equations for curvature-isocurvature fluctuations in Fourier modes is given by (see App.A or [52,130] for details) and where we defined as usual The relative velocity V rel can be computed from S by Let us note here that in the matter dominated universe where c 2 s → 0 the general "growing mode" solution to Eqs. (2.17) is Φ = 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 Φ i the initial values of isocurvature and curvature in the far past, formally when a → 0. The conclusions do not change if a has a non-zero value as the solutions are attractors [66].
a. Superhorizon fluctuations (k ≪ H): If we drop all terms containing k in Eqs.(2.17) and (2.18), and re-write them using a as a time variable, we find that the general solutions are given by [7,132,135] where for compactness we defined ξ = a/a eq .We see that while isocurvature remains constant on superhorizon scales, curvature fluctuations Φ are not.For the adiabatic component we have that (a ≫ a eq ) , where the factor 9/10 actually comes from the conservation on ζ and the relation ζ = 5+3w 3+3w Φ for constant w = p/ρ.For isocurvature fluctuations we instead have Thus, while we have vanishing initial curvature perturbation, i.e.Φ iso (a → 0) → 0, it grows as a/a eq and saturates to 1/5 of the initial isocurvature after matter-radiation equality.
b. Fluctuations during radiation domination (a ≪ a eq ): To study the evolution of general fluctuations during radiation domination, we shall take the limit τ ≪ τ eq (a ≪ a eq )in Eqs.(2.17) and (2.18).At leading order, they are given by [52] and where we defined for compactness In these new variables, the limit of interest is given by x ≪ κ (or k eq τ ≪ 1).Due to the length of the solutions we treat the initial adiabatic and initial isocurvature cases separately below.
For initial curvature fluctuations, Φ i ̸ = 0 and S i = 0, we solve the leading order terms in Eqs.(2.25) and (2.26), namely we first solve the homogeneous equation for Φ and plug it in in the equation for S.This yields where γ E ≈ 0.577 and Ci(x) is the cosine integral function.Note that Eq. (2.28) is the standard solution for adiabatic perturbations in the radiation universe.Φ is first constant and decays as x −2 once a given mode enters the sound horizon, that corresponds to c s x > 1.We also see that S is negligible for c s x < 1 as it is proportional to x 4 but grows as log(c s x) for c s x > 1.Such logarithmic grows is due to the fact that matter perturbations grow logarithmically during the radiation dominated universe [136].It should be noted that one may also use Eq.(2.29) to compute the next order solution to Φ.We refer the interested reader to Ref. [7] for general solutions.
For the initial isocruvature case, Φ i = 0 and S i ̸ = 0, we follow Ref. [52].In this case, we see that we can expand the solution of Φ and S in powers of κ −1 for κ ≫ 1 and in powers of x/κ for κ ≪ 1, as e.g. S = S i + S 1 + ... and Φ = Φ 1 + ... etcetera.The leading contribution to Eq. (2.25) Φ is a then a constant S and we may also compute the effects of the leading solution to Φ to the next leading solution to S. Doing so we find [52] Φ where Si(x) is the sine integral function.Looking at x ≪ 1 we see that initially the curvature perturbation grows as Φ iso ∝ x, reaches a maximum at around c s x ∼ 1 and then decays as x −2 .It is also interesting to note that for c s x ≫ 1 we have Φ ad ⊃ sin(c s x) and Φ iso ⊃ cos(c s x), recovering the well-known result that adiabatic and isocurvature initial conditions give an out of phase curvature fluctuations.Isocurvature S is constant for c s x < 1 and then grows with x for c s x > 1.Interestingly, it is possible that S reaches a high enough amplitude for PBHs to form [51].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 (2.32) At times τ > τ 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 τ 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 NL > 1) translates into an upper bound on the initial isocurvature, namely FIG. 1. Numerical solutions to the evolution of curvature fluctuations for adiabatic and isocurvature initial conditions in terms of x = kτ , respectively in purple and orange lines.We normalized the amplitude of the initial conditions to Φ i = 1 and S i = 1.On the left we show the behavior on scales with k ≪ k eq .For the numerical solution we chose k = 10 −3 k 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 ≫ k eq , in particular k = 500k eq .In this case the lines also saturate to a constant after matter-radiation equality.However, for k ≫ k eq the amplitude after τ eq is of the order of (k eq /k) 2 in both cases.
for a given isocurvature mode with wavenumber k.
c. Fluctuations during matter domination (a ≫ a 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 ≪ k eq its value is given by Eqs.(2.23) and (2.24) respectively for adiabatic and isocurvature initial conditions.For k ≫ k eq the amplitude has a suppression factor proportional to (k/k 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] Φ iso (a ≫ a eq ) We show the results of numerical integration in Fig. 1.In the numerical results we find that the coefficient for k ≫ k eq in Eq. (2.34) is close to 1.For easier comparison with the literature though, we maintain the coefficient of Eq. (2.34) as it only introduces small errors.

C. 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 (2.9) which is given by where S ij (2.12) using the linear Einstein Equations is given by (2.36) It will be more convenient to work in Fourier modes, where we give a primordial initial amplitude of S k (0) and/or Φ k (0) and describe the evolution of a given mode by a transfer function, say T Φ (kτ ).
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 Φ.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.(2.35) reads with where Φ q (0) and S q (0) respectively refer to initial adiabatic or initial isocurvature fluctuations, and we defined In Eq. (2.39) we introduced for later use the notation Also, in Eq. (2.39), T V rel (x) is the transfer function of V 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. (2.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 [52,137,138], namely where P S (k) is the initial dimensionless spectrum of isocurvature fluctuations, defined by For adiabatic initial conditions P S (k) should be replaced by P Φ (k) at the initial time.In Eq. (2.41) the notation x c refers to a time when the GW is well inside the horizon such that the spectral density is constant [139].In Eq. (2.41) we introduced the kernel which incorporates the information from the transfer functions and it is given by where G(x, 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 Then, we may split the kernel (2.43) into a sine and cosine terms as where This allows us to take the oscillation average of the spectral density, which is given by where we took the limit x → ∞ 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 (2.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 Ω GW,c (2.41) and the spectral density of GWs measured today.A straightforward relation can be found in the case when the radiation domination era where Ω GW,c is computed corresponds to the standard radiation dominated stage, i.e. prior to BBN and CDM domination.In that case, taking into account the redshifting of the GW energy density one finds that [48] Ω GW,0 h 2 ≈ 1.62 where Ω GW,0 is the energy density fraction of GWs evaluated today and h = H 0 /(100km/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.[139,140].

III. GRAVITATIONAL WAVES FROM CDM ISOCURVATURE
Let us consider that the matter component with isocurvature fluctuations is the CDM.This means that the initial fraction of CDM β is fixed its the current abundance, which according to Planck [4] is Ω CDM,0 h 2 ≈ 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 eq ≈ 3400 and that the comoving size of the universe at that time was k eq ≈ 0.01 Mpc −1 .To grasp the magnitude suppression factor k eq /k, let us write down the GW frequency in terms the relevant wavenumbers, that is 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 nHz and 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 τ NL .For a detailed discussion on the validity of perturbations see the appendix of Ref. [52].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. (2.30) during the radiation dominated epoch.Plugging in Eq. (2.30) into the integrals (2.46) we find that Source term (2.39) for induced GWs for adiabatic and isocurvature initial conditions, respectively in purple and orange, in terms of x = kτ .For simplicity we chose u = v = 1 but the main behavior is independent on such particular values.Here we consider τ ≪ τ 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.and where we took the limit x → ∞ for analytical simplicity.With the kernels Eqs.For illustration, we also compare the source term (2.39) between initially adiabatic and isocurvature fluctuations in Fig. 2.
It is interesting to note that the explicit κ dependence in the kernel (3.2) and (3.3) can be absorbed via the definition of an effective power spectrum of curvature perturbations, namely Although this is a mere redefinition, it illustrates the differences with the adiabatic initial conditions.For instance, after using P 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 [52].The detailed shape is of course different.Also note that the GW spectrum has a global suppression proportional to κ −4 , consistent with our expectations that Φ ∝ ρ m /ρ r × S ∝ κ −1 S for isocurvature initial conditions when evaluated at horizon crossing.We proceed to discuss some concrete examples.
A. GWs induced by a peaked primordial CDM 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 CDM isocurvature fluctuations peaks at given scale, say k p .For concreteness and to follow previous works we consider a log-normal peak given by where A S is the amplitude and ∆ is the logarithmic width of the peak.This is standard practice when dealing with GWs induced by primordial adiabatic fluctuations (see e.g.[141]) and it is well motivated from inflationary models [25,.Similar conclusions extend to the isocurvature case [167,168].We note that the formal limit of ∆ → 0 corresponds to a Dirac delta peak, namely A log-normal is considered to be "sharp" if ∆ ≲ 0.1 and broad otherwise [141].We take a similar ansatz for primordial curvature fluctuations Φ i .The Dirac delta case is convenient as it allows for an analytical expression for Ω GW,c simply by replacing u = v = k p /k in Eqs.(3.2), (3.3) and (2.41).Explicitly, we obtain where the sharp cut-off at k = 2k p comes from momentum conservation.From Eq. (3.7) we see that, as in the initially adiabatic case (see e.g. the discussions in Ref. [48]), the induced GW spectrum has a resonant peak at k = 2c s k p and decays as k 2 ln 2 k in the low frequency tail, namely for k ≪ k p [169][170][171].However, contrary to the adiabatic case, there is no destructive interference with vanishing GW spectrum at k = √ 2c s k p .This is because adiabatic initial conditions lead to coherent, order unity, oscillations of Φ while for isocurvature initial conditions the oscillations of Φ 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 ∼ 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 FIG. 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 A S = 1 and A Φ = 1 in Eq. (3.5) for P S as well as P Φ , 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 s k 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 ∆ = 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 ∆ as explained by Eq. (3.4).
This means that the amplitude of the initial power spectrum of isocurvature fluctuations must be very large to partly compensate for the suppression, roughly A S ∝ (k p /k eq ) 2 ≫ 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 ρ m /ρ r during the radiation dominated era.The product of initial isocurvature times ρ m /ρ 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 ∆ ≲ 0.01, except that for scales k < ∆ k p , where the slope transitions to a k 3 scaling, as in the adiabatic case [141].For ∆ ≳ 0.01 similar conclusions apply but because of the κ −2 dependence in P eff Φ (3.4), 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.

B. 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. (2.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 ∝ δρ m /ρ m , and while P S ≫ 1 is certainly possible, the probability distribution of S cannot be Gaussian as S > −1 (i.e.ρ m + δρ m > 0).As an illustration we chose A Φ = 0.01 and A S = 0.05 (k p /k eq ) 2 ∼ 10 18 to fit the PTA data by eye.The peak is chosen at f p ∼ 10 −7 Hz (k p ∼ 10 7 Mpc −1 ).Gray violins indicate the recent NANOGrav results [78] and in blue and orange we show power-law integrated sensitivity curves [172] for LISA and Taiji [120,121].See [126][127][128]173] for the sensitivity curves.The horizontal blue and purple dashed lines respectively show the current constraint from BBN [174][175][176] and future constraints from CMB-S4 experiments [175,177].We also include future sensitivity of µ-Ares [119] and Lunar Laser Ranging from Ref. [178] (see also [179]).GW detectors such as µ-Ares may tell apart whether the GWs were induced by adiabatic or isocurvature initial conditions.
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 ⟨S⟩ = 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. [52] 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 CDM isocurvature induced GWs.It would be interesting to study a concrete case and confirm these expectations.
It is also interesting to consider the hypothetical case where the possible GW background signal reported by the PTA collaborations [68][69][70][71][72][73][74][75][76][77] 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 CDM isocurvature fluctuations.From the analysis of GWs induced by primordial adiabatic fluctuations we have that, if given by a Dirac delta spectrum, one needs A Φ ∼ 10 −2 − 10 −1 and a peak position around f p ∼ 10 −7 Hz [78].This means that for the primordial CDM isocurvature one requires A S (k eq /k p ) 2 ∼ A Φ , which roughly corresponds to A S ∼ 10 18 −10 19 .While this value of A S is certainly big it is still within the validity range of our induced GW spectrum calculations as A S (k eq /k p ) 2 ≪ 1.With future detectors in the µHz window, such as µ-Ares [119] 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).

IV. 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 [180,181].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.[129,182], 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.

A. 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 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 f , via M PBH,f ∼ 4πM 2 pl /H 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.[129,130] for the details.
After formation, PBH make a fraction β of the total energy density at formation.It is interesting to relate the β to the number density of PBHs, namely where we used that ρ PBH,f = M PBH,f × n PBH,f and ρ r,f ≈ 3H 2 f M 2 pl .Thus, β 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 PBH ∝ a −3 , the PBH mass decreases due to Hawking radiation.One finds that the mass decays as [129] where t eva is the time of evaporation and is approximately given by where fs = 10 −15 s are femtoseconds.Here t is cosmic time related to conformal time by dt = adτ .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 eva ∼ 1/t eva ) we can estimate the temperature of the radiation fluid filling the universe right after evaporation.In doing so, we find that From the Hubble parameter we can also compute the comoving size of the Hubble horizon at evaporation, namely [130] 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 eva > 4 MeV and translates into M PBH,f < 5 × 10 8 g.We can also compute the minimum value of β such that PBHs eventually dominate the universe before evaporating, which is given by 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 Then, the energy density of the PBH fluid is gradually transferred to radiation and obeys where ˙≡ d/dt.The equation for the radiation energy density is similar but with the opposite sign in front of Γ, 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 ρ PBH and ρ r for a particular example.
The main point that we want to emphasize here is that fluctuations with k ≫ Γ 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 [129].This translates into a suppression factor for curvature fluctuations with k ≫ Γ after evaporation given by ), 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 β = 10 −4 for illustrative purposes.See how PBH domination can last for several e-folds and that evaporation is faster than one e-fold.
where "lRD" means late radiation domination and "instant" refers to the amplitude of Φ obtained by an instant matching from matter to radiation domination.Note that the exponent 1/3 in Eq. (4.9) is directly related to the exponent in the decay of the PBH mass (4.2).The relation between Φ and ρ PBH follows from the Poisson equation on very subhorizon scales, that is k 2 Φ ∼ a 2 ρ PBH δ PBH .We must take the suppression factor S Φ (k) (4.9) into account when computing any transfer function for the curvature perturbation.We suggest to read Ref. [129] for more details on the calculations.

B. 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 [65] On scales larger than d f we have a coarse grained fluid picture but on scales smaller than d 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 β 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 [65] δn where the factor a 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 inter-PBH comoving separation, namely where we used that k f = H 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 δρ PBH,f + δρ r,f = 0.This implies that where we used that initially ρ PBH ≪ ρ r .From Eq. (4.11), the dimensionless initial isocurvature power spectrum then reads This is the initial isocurvature that will eventually generate induced GWs [65,130].
C. 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 Φ.As we will be interested in the GWs sourced by curvature fluctuations after evaporation and on very small scales (since k UV ≫ k eva ), we just have to follow Φ 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. [65] 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 β, that is where k eq = H eq .It then follows that we have the resulting hierarchy: To put it into words, the initial fluctuations on scales similar to k UV are initially superhorizon (they are on scales larger than k f ).But, these scales around k UV enter the horizon well before matter-radiation equality and so k UV ≫ k 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. (2.34), namely the transfer function for isocurvature initial conditions, together with the suppression factor (4.9), which yields To compute the GWs induced after evaporation we have to continue the constant value of Φ during the PBH (matter) domination to the late radiation domination.After matching, we find that where Now, with the evolution of Φ 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. [129,130] for the details.
The most important point is to realize that, while the amplitude of Φ in Eq. (4.17) starts from a constant and quickly decays, the amplitude of its time derivative, Φ ′ , begins at zero (by continuity) and suddenly jumps to an amplitude proportional to k/k eva , which is huge for k ∼ k UV .The reason for this jump is that during the matter dominated era Φ remained constant but the relevant scales became more and more subhorizon.Then, at the late radiation domination, Φ resumes its decay but with very fast oscillations.Yet, since it has not decayed during the matter dominated phase, the amplitude of Φ ′ is suddenly very large by a factor a eva /a UV ∼ k UV /k eva (the supposed decay if there never were a matter dominated era).
Another explanation for the enhancement is given in Ref. [182].During the matter dominated phase PBH density fluctuations grow, as standard CDM 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 App.A V ∝ Φ ′ /H).Inspecting the source term (2.39) we see that the main contribution to the kernel (2.43) comes from Φ ′ and, basically, we can approximate the kernel by where G lRD (x, x1 ) is given by (2.44).After integration and selecting only those terms related to the resonant peak (those at c s (u + v) ∼ 1 where the frequency of the source term is equal to that of the tensor mode) we arrive at To compute the induced GW spectrum we plug in the averaged kernel (4.20) into Eq.(2.41) and integrate.However, for sufficiently peaked spectrum the integral is well approximated by a power-law with a sharp cut-off [129][130][131].A power spectrum is considered to be peaked enough if the integrand in Eq. (2.41) grows for growing u and v. Since it contains (uv) 2 ∼ v 4 , the effective spectral index of P Φ (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 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 where the peak amplitude of the GW spectrum is approximately given by .
Furthermore, using Eqs.(4.12) and (4.5), we find that the peak GW frequency is located at Thus, on one hand, we see that for 5 × 10 8 g > M PBH,f > 10 4 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 This is the best constraint we have on the fraction of PBHs that evaporated before BBN.We show the resulting bounds on β and an example of the GW spectrum from the PBH dominated universe in Fig. 6.

D. Remarks, issues and future work
In our calculations we have mainly focus on the curvature perturbation Φ right after evaporation and we trusted our linear perturbation theory since Φ ≪ 1.However, when we look into the M PBH,f = 10 6 g , β = 2 × 10 −8 FIG. 6.On the left we show in the shaded orange region the parameter space for the initial fraction of PBHs, β, such that PBHs dominate the early universe (β > β min ; magenta line) and that the induced GWs satisfies BBN bounds (β < β 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 PBH,f = 10 6 g and β = 2 × 10 −8 .We also show the power-law integrated sensitivity curves [172] for DECIGO, Einstein Telescope (ET), Cosmic Explorer (CE), Voyager and LIGO A+ experiments [126][127][128]173] as well as the upper bounds from the LIGO/Virgo/KAGRA collaboration [118].The horizontal lines show BBN [174][175][176] and CMB-S4 experiments constraints, respectively in blue and purple [175,177].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.evolution of matter density fluctuations, which on small scales are related to Φ by the Poisson equation 2k 2 Φ = a 2 δρ, we find that δρ/ρ at evaporation is larger than unity [65,130].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 nonlinear regime.But, the fact is that we do not expect the GW production to stop by non-linearities.Thus, we take our result, Eq. (4.23), as a rough order of magnitude estimate.For an interesting proposal in terms of turbulences in the radiation fluid after evaporation see Ref. [183].
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 [129].For works including both PBH isocurvature and adiabatic induced GWs see Refs.[184][185][186].For articles investigating the effect of a PBH dominated stage on the GW spectrum from cosmic strings see Refs.[187][188][189].Other applications of the GWs from the PBH dominated universe can be found in Refs.[190][191][192] 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 [131,175,[193][194][195][196][197].For a similar application in the case where inflaton oscillons dominate the universe after inflation see Ref. [208].
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 [131] (although see Ref. [198] 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 with some arbitrary cut-off at k UV , then the resulting induced GW spectrum after evaporation reads [131] Ω GWs,rh (k where It is also important to note that if transition is rather gradual, the enhancement of the induced GW spectrum disappears [199,200].This is particularly important if PBHs have a relatively broad mass function [129].

V. 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 CDM 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 ρ m /ρ 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 ≫ 1) is consistent with cosmological perturbation theory, it has the caveat that isocurvature fluctuations must be extremely non-Gaussian [52].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 CDM 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 (4.23) 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 [201] at the end of evaporation have for the induced GW signal.See refs.[202][203][204][205][206] for reviews and models on black hole remnants.Recently in Ref. [207], it has been shown that if such remnants constitute the totality of the CDM then the initial PBH mass must be 5 × 10 5 g and the induced GW spectrum must peak at around 100 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 CDM.
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.[66,67], 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 § IV.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 § III, and might give another contribution if solitons dominate the early universe before decaying.Such generic prediction has been named "Universal Gravitational Waves of Solitons" [66,67].
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.
(3.2) and(3.3)we are ready to compute the isocurvature induced GW spectrum via Eq.(2.41).Let us emphasize that Eqs.(3.2) and (3.3) 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 App.B.

FIG. 4 .
FIG.4.Example of induced GWs from a Dirac delta spectrum for primordial curvature and isocurvature fluctuations.As an illustration we chose A Φ = 0.01 and A S = 0.05 (k p /k eq ) 2 ∼ 10 18 to fit the PTA data by eye.The peak is chosen at f p ∼ 10 −7 Hz (k p ∼ 10 7 Mpc −1 ).Gray violins indicate the recent NANOGrav results[78] and in blue and orange we show power-law integrated sensitivity curves[172] for LISA and Taiji[120,121].See[126][127][128]173] for the sensitivity curves.The horizontal blue and purple dashed lines respectively show the current constraint from BBN[174][175][176] and future constraints from CMB-S4 experiments[175,177].We also include future sensitivity of µ-Ares[119] and Lunar Laser Ranging from Ref.[178] (see also[179]).GW detectors such as µ-Ares may tell apart whether the GWs were induced by adiabatic or isocurvature initial conditions.

FIG. 5 .
FIG. 5. Energy density fraction Ωy = ρ y /(3H 2 2pl ), 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 β = 10 −4 for illustrative purposes.See how PBH domination can last for several e-folds and that evaporation is faster than one e-fold.