Thermal History of the Early Universe and Primordial Gravitational Waves from Induced Scalar Perturbations

We study the induced primordial gravitational waves (GW) coming from the effect of scalar perturbation on the tensor perturbation at the second order of cosmological perturbation theory. We use the evolution of the standard model degrees of freedom with respect to temperature in the early Universe to compute the induced gravitational waves bakcground. Our result shows that the spectrum of the induced GW is affected differently by the standard model degrees of freedom than the GW coming from first order tensor perturbation. This phenomenon is due to the presence of scalar perturbations as a source for tensor perturbations and it is effective around the quark gluon deconfinement and electroweak transition. In case of considering a scalar spectral index larger than one at small scales or a non-Gaussian curvature power spectrum this effect can be observed by gravitational wave observatories.


I. INTRODUCTION
The first observation of gravitational wave (GW) event from binary black holes merger and also binary neutron stars by LIGO opened a new way to our understanding about the Universe [1][2][3]. Before using electromagnetic waves we could see the physical phenomena in the Universe, however now somehow we can listen to them in case they produce gravity waves. Beyond GW produced by astrophysical events at the recent era of cosmic evolution, it would be interesting to listen to the early moments of the Universe to detect the predicted GW by cosmological and particle physics models [4,5]. The inflationary scenario as a theory to explain the moments after the big bang and explain the current cosmological observations can produce that kind of primordial gravitational waves (PGW) [6,7]. In principle if we have detectors of high sensitivity, one can observe such PGW. However, their detection depends on the first order tensor to scalar ratio "r" which appears directly in the predicted observable relic density of tensor perturbation at the first order of cosmological perturbation theory. Since the cosmic microwave background (CMB) experiments like Planck have only put an upper bound on the tensor-to-scalar ratio and some theoretical predictions assume very tiny values close to zero for this quantity, it motivates to consider the second order terms in the cosmological perturbation theory induced by scalar perturbations. This has been done in the literature [8][9][10][11], and gained much interest recently [12][13][14][15][16][17][18][19][20][21][22][23][24][25] due to the first direct discovery of GW signals from binary mergers.
The predicted shape for the induced spectrum using the standard model (SM) equation of state at different scales or frequency and the evolution of induced tensor power spectrum should be roughly like the first order one. However, this has not been checked numerically and the goal of this paper is to do exactly this. Studying induced PGW is also interesting when there is non-Gaussianity on scales smaller than the cosmic microwave background (CMB) [20,54]. This can boost the spectrum of induced PGW to large values of the relic density which might be accessible by near future experiments [18][19][20] especially by pulsar timing arrays. Also, it would be one of the real chances to test particle physics models and early Universe cosmology based on our current data from CMB and standard model of particle physics and ΛCDM cosmology in the case that tensor-to-scalar perturbation ratio is much smaller than the current bound from Planck [26] or from future CMB experiments like CMBS4 [27].
Currently, GW experiments can mostly detect the GW spectrum from astrophysical events with a strong enough amplitude. However, future detectors will be able to probe much smaller signatures from possible phase transitions and the PGW in the early Universe cosmology. Some of these planned GW missions are the Laser Interferometer Space Antenna (LISA) [28], the Einstein Telescope (ET) [29], the Deci-hertz Interferometer Gravitational wave Observatory (DECIGO and B-DECIGO) [30,31], the Big Bang Observatory (BBO) [32], the Square Kilometer Array (SKA) telescope [33], the North American Nanohertz Observatory for Gravitational Waves (NanoGrav) [34], and the European Pulsar Timing Array (EPTA) [35]. Moreover, LIGO can be sensitive to the predicted GW for the early Universe depending on the model [36,37].
For a pre-BBN cosmology which is dominated by some matter with an equation of state different from radiation the spectrum of the first order and the induced PGW will be different as studied in refs. [9,12,13,15,16,[38][39][40][41]. Such models will not be the focus of the current paper.
The value of secondary produced PGW depends on the amplitude of scalar perturbations which is observed at the s (solid red curve) and equation of state parameter ω (blue dotted curve) with respect to temperature T including all SM degrees of freedom at relevant temperatures using the result of ref. [46]. The dashed line shows the ideal gas case c 2 CMB scale [26]. Since the scalar spectral index at large scales is highly constrained by CMB and is close to one, it is not expected to deviate significantly from a flat scalar power spectrum for the standard case. However, it might be different at smaller scales e.g. at pre-BBN era. This makes the induced PGW different from the first order one which has different free parameters i.e. the scalar-to-tensor ratio and the tensor spectral index [38,42,43].
In the next section (II) we investigate the thermal history of the early Universe including SM particles. In Sec. III the set of equations related to first order and second order perturbations and the induced PGW will be discussed. We compute the relic density of the induced PGW for scale invariant and scale dependent scalar perturbations in Sec. IV. Finally, we summarize our results in Sec. V.

II. THERMAL HISTORY OF THE EARLY UNIVERSE AND STANDARD MODEL PARTICLES
The effect of standard model degrees of freedom on the (first order) PGW especially around the quark gluon deconfinement and electroweak transition are studied in the literature [38,[42][43][44][45]. The evolution of degrees of freedom influences the relation between the scale factor and the temperature in the early Universe cosmology [45,46]. However, this can be interpreted as an evolution of the equation of state parameter. The equation of state parameter and speed of sound can be driven by other thermodynamics variables as where the energy density and entropy density are defined as and the effective relativistic degrees of freedom (DoF) for the energy and the entropy density are denoted by g eff and h eff . The contribution of every relativistic fermion (boson) to the effective DoF is 7/8 (1). For radiation or any ideal relativistic fluid c 2 s = ω = 1/3. We use the result of ref. [46] to compute the equation of state parameter and speed of sound for SM degrees of freedom as plotted in fig. 1. This figure shows that around temperatures of 150 MeV and 100 GeV, due to quantum chromodynamics (QCD) and electroweak crossover transitions ω and c 2 s will be smaller than 1/3. Due to the interaction between particles the thermal bath of the early Universe deviates from an ideal relativistic fluid. Every extra degree of freedom e.g. due to dark matter or any beyond the SM physics in the early Universe can also lead to small effects on these quantities. Especially, considering a highly populated sector like supersymmetry, can modify the result of fig. 1 and add a new valley at temperatures around 1 TeV or higher depending on the scale of supersymmetry (see refs. [43,47,48] for supersymmetric DoF and PGW).

III. TENSOR AND SCALAR PERTURBATION EQUATIONS
At second order of the cosmic perturbation theory one can consider the following metric [8,10,12] where a is the scale factor and η the conformal time. Assuming the metric does not contain anisotropic stress (Φ = Ψ) the scalar perturbation and tensor perturbation are denoted by Φ and h, respectively. The metric in eq. (3) is considered in Newtonian gauge. Here we mainly follow ref. [12] for the evolution equations of induced scalar and tensor perturbations. The evolution of the tensor power spectrum sourced from scalar perturbation can be obtained by solving the following equation in Fourier space [8][9][10]12] where the conformal Hubble parameter is denoted by H = aH and the derivative with respect to the conformal time by = d/dη. The Hubble parameter can be evaluated from the Friedmann equation H 2 = (8π/3M 2 P l )ρ tot , with M P l being the Planck mass. Also, the entropy conservation in the early Universe is considered 1 . To consider the whole evolution of DoF for finding a(η) one should solve a = a 2 H numerically.
The source term on the right hand side of eq. (4) assuming Φ = Ψ can be evaluated as [8][9][10]12] where the polarization tensor for GW is shown by ij ( k). In eq. (6) the terms depending on c 2 s cancel out so that the expression depends only on ω explicitly. We use the Green's function method to solve the differential equation for tensor perturbations. Then the solution in the comoving frame will be [12] where we make use of the Green's function as the solution of the following equation [8,12] For a fixed ω the scale factor evolves as a ∝ η 2 3ω+1 . Then eq. (8) will be given by To compute the Green's function we use the method of ref. [8]. If the functions G 1 and G 2 are two homogeneous solutions of eq. (8) for a mode k then the Green's function will be defined as where in practice we compute the numerical form of this function (assuming G 1 (η) = 0, G 1 (η) = 1; G 2 (η) = 1, G 2 (η) = 0) for each mode k or GW frequency f = 2π/k. The evolution of scalar perturbations in cosmology can be described with the differential equation [12,49] where δp = c 2 s δρ + τ δS. By assuming c 2 s = ω = const. and vanishing entropy perturbation in cosmology δS = 0 one gets [12] Φ k (η) + 6(1 + ω) 1 To find the evolution of temperature with respect to scale factor in the standard cosmology after reheating the entropy conservation should be assumed. This can be taken into account by the relation In the limit c 2 s → ω our numerical results match with the analytical solution of eq. (12) as given in refs. [8,12]. However, for a precise numerical calculation one should distinguish between the speed of sound c 2 s and the equation of state parameter ω as shown in fig. 4 and consider them separately in the computation. The correlation function for curvature perturbations is derived as shown in ref. [12] and gives where the superhorizon value of φ k is evaluated from the relation Φ k = Φ(kη)φ k . The transfer function Φ(kη) reaches one before the moment of horizon crossing (k = a hc H hc ). We will consider the evolution of the degrees of freedom of SM in the result of eq. (11) and other relevant equations in the following. The induced tensor power spectrum from scalar perturbations is determined via [8,10,12] where v = | l|/| k| and u = | k − l|/| k|. In eq. (14) the definition of the function I reads [10,12] where the source function f based on the equation of state parameter ω and transfer functions Φ is obtained from [10,12] f The terms depending on c 2 s for the case Φ = Ψ in the original equation of refs. [8,10] cancel out. Then only the terms on the right hand side of eq. (16) are left. By changing the variables u and v to u = (p + q + 1)/2 and v = (p − q + 1)/2 [12] the numerical solution of the power spectrum turns out to be computationally faster.
For different values of ω shown in fig. 3 we derive the Green's function and Φ from eqs. (10) and (12) and computed then the function f from eq. (16). Then, the function I is computed analytically from eq. (15) and used in eq. (14) for numerical integration. As an example, the Green's function for ω = 1/3 considering a ∝ η and a = 0 from eq. (9) can be evaluated as G k (η,η) = sin(η −η)/k [12]. Also, the scalar transfer function for a radiation like fluid by solving eq. (12) and assuming (Φ(kη → 0) → 1) can be shown to be [12,50] Φ(kη) = 9 (kη) 2 The explicit form of the function I is complicated and long and will not be presented here (see refs. [10,12]). In ref. [12] the limits of the function I in a radiation dominated Universe are found to be I(u, v, x → 0) ≈ x 2 /2 and I(u, v, x → ∞) ≈ I(u, v)/x 2 . Using the latter limit the evolution of the tensor power spectrum is estimated to be P(η, k)/A 2 R ≈ 19.7/x 2 [12]. This expression is vastly used afterwards to estimate the value of the GW relic density of today. However, the value 19.7 does not match the results for ω = 1/3 shown by the solid black line in fig. 3 based on a complete numerical calculation. This will be explained in more details in the next section.

IV. PRIMORDIAL GRAVITATIONAL WAVES AND EXPERIMENTAL CONSTRAINTS
The scale invariant scalar power spectrum defined in eq. (13) can be written as [51] Moreover, the power law scalar power spectrum can be defined as [51] with the pivot scalek = 0.05 Mpc −1 and the scalar spectral index n s . From Planck satellite observations the value of the scalar spectral index has been determined to be n s = 0.965 ± 0.004 [26]. Also, the value of the curvature perturbation amplitude at the CMB scale is A R 2.1 × 10 −9 [26]. In addition, the scale invariant tensor power spectrum at first order is defined as [51] P (1)  As ω → 0 the tensor power spectrum diverges as the upper limit of first integral in eq. (14) becomes infinite [12].
Then the tensor to scalar ratio at first order is [51] The current limit from the Planck collaboration is r 0.07 [26]. In fig. 2 we assume that the parameter ω is fixed and constant in the calculation of eqs. (8), (12), and (14). The scaled tensor power spectrum (scaled by the square of the scale invariant curvature power spectrum: P T (η hc , k)/A 2 R ) at horizon crossing is computed for the range of 0 < ω < 1 where the minimum of the tensor power spectrum appears at ω min = 0.175. As fig. 2 shows and as it is analytically calculated in refs. [9,12,41] if ω → 0 the tensor power spectrum diverges P T (η ≥ η hc , k) → ∞. The final result for P T (η hc , k)/A 2 R in the range 0 ≤ kη ≤ 10 and different ω's is shown in fig. 3. The function I inside the integral of eq. (14) is highly oscillatory. It gets its maximum value slightly after horizon crossing (2 kη 4 depending on the value of ω) then decays, as ∝ 1/(kη) 2 during radiation domination [12].
For the induced tensor perturbation in the computation of the tensor power spectrum after horizon crossing until today the WKB method can not be used to match the result at some point and then extend it until today, since (to the best of our knowledge) there is not any well defined shape of the transfer function during radiation domination for the induced PGW like A sin(δ + kη)/kη as assumed in refs. [38,42,43] for the first order PGW.
The power spectrum, shown in fig. 3, increases after horizon crossing to a maximum value then dilutes with the expansion of the Universe. The numerical calculation of such a power spectrum is difficult and time consuming since it includes various numerical integrations over oscillatory functions. In practice it is very difficult to solve for a large range of kη from zero to today kη 0 . So one should calculate until horizon crossing or until the peak for ω = 1/3 in fig. 3 then, use 1/(kη) 2 as the dilution factor due to the expansion of the Universe.
The tensor power spectrum can also be defined in terms of the transfer function P(η hc , k) ∝ T (η hc , k) 2 = k 2 T (η hc , k) 2 . It scales like 1/a 2 after horizon crossing. For the radiation domination case (without the change of DoF) the scaled tensor power spectrum in fig. 3 at horizon crossing and at its first peak (x 1st,peak 2.17) are P T (η hc , k)/A 2 R 0.87 and P T (η 1st,peak , k)/A 2 R 1.75, respectively (P T (η 1st,peak , k)/P T (η hc , k) 2). It is expected that the evolution of the DoF in the early Universe plays a role in the computation of eqs. (8), (11) and (14) so we checked it. We have included all DoF of the SM to find the scaled tensor power spectrum P T (η hc , k)/A 2 R . Then the shape of the spectrum will have a correction factor, shown in fig. 4, due to the evolution of different functions inside the integral of eq. (14) by a change of the DoF over time. This fact originates from the combination of all functions in the integral a(η)/a(η), f (u, v, kη), and G(η,η) due to carrying the retarded effects from previous times. In practice, an analytical formula to show the explicit form of fig. 4 has not been found. However, it roughly behaves like the evolution of dω/dT including the retardation effects due to the Green's function. The inflection points in the correction factor of secondary PGW spectrum roughly show the points of the crossover transition having the smallest values of c 2 s and ω at the relevant frequencies (3 × 10 −9 Hz and 10 −6 Hz). In the case that one adds the supersymmetric DoF or any highly populated sector around the TeV scale, there will be another inflection point at higher frequencies around 10 −4 Hz.
The density of PGW per frequency (wave number) from tensor perturbation at first order in the case of eq. (4) is source free and can be estimated as [38,42,43] where T (η hc , k)/(kT (1) (η hc , k)) 2 for the first order PGW when T (1) (η hc , k) 2 is the 1st order tensor perturbation transfer function. Here we focus on frequencies larger than 3 × 10 −11 Hz which are equivalent to temperatures larger than 4 MeV. Consequently, we do not consider the effects of neutrinos and photons free streaming which appear at smaller frequencies [42,52,53] as a new source term on the right hand side of eq. (4). The relic PGW at today from induced scalar perturbation can be obtained from [50] Ω GW (η 0 , k)h 2 = 1 24 the correction factor is defined as Z(k) = P T (η hc , k)/P T (η hc , k high ) shown in fig. 4 as blue solid line (P T (η hc , k high )/A 2 R 0.87). The shape of this function completely depends on the details of the functions inside the integral of eq. (14). The function Z(k) already corrects the spectrum on the order of 10% especially around the QCD epoch which is equivalent to a frequency of 3 × 10 −9 Hz. As mentioned earlier the evolution of DoF of SM particles is taken from the results of ref. [46]. Using the numerical factor 2 × (2.17) 2 in eq. (23) from the fact thatx 1st,peak 2.17, the tensor power spectrum from scale invariant curvature power spectrum during radiation domination evolves as [P T (η 1st,peak , k)x 2 1st,peak ]/[P T (η hc , k)x 2 hc ] × 1/x 2 2 × (2.17) 2 /x 2 . The results for the first order and second order PGW are shown in fig. 5 in addition to the different experimental constraints [28][29][30][31][32][33][34][35][36][37] and the BBN bound [12] as outlined in the legend of the plot (the suffix number for the SKA experiment shows the constraints based on the number of years in operation). In case of scale invariance the relic of the first order PGW can have any value from ∼ 10 −16 down to the lower limit by the induced PGW of ∼ 10 −23 depending on the value of r. If r 10 −9 then the scalar induced contribution on the PGW will be dominating. The scale dependent induced PGW for scalar spectral indices of n s = 0.96 and 1.5 are also plotted in fig. 5. For the case n s > 1 the spectrum can be constrained by current and future experiments.
The curvature power spectrum can deviate from a flat shape at scales smaller than the CMB scale especially if one assumes that there is some Gaussian or non-Gaussian source for it from the curvature bispectrum [18,20,54]. This can enhance the spectrum by some peaks of the relic densities at frequencies accessible by pulsar timing arrays or other experiments. The effect of the evolution of the degrees of freedom in the early Universe also appears on the secondary PGW if there is some source of non-Gaussianity [18,20]. To represent this effect we assume first that the following Gaussian curvature power spectrum (for k in the logarithmic scale) around a specific modek with width σ [18][19][20] is parameterised by where A G is the Gaussian amplitude that can depend on the characteristic scalek and other model dependent variables. Here we simply assume it to be constant to test for its observational consequence. By considering the presence of non-Gaussianity in the curvature power spectrum one gets [18][19][20] where the nonlinear factor for the non-Gaussian curvature perturbation is denoted by F N L . We choose different benchmark points fork and F N L to see how its consequence will be on the PGW. Also, we consider the evolution of DoF to check how it influences the induced PGW by non-Gaussian power spectrum. This is also shown in fig. 4 for the scaled tensor power spectrum assuming F 2 N L = 10, A G = 10 −2 , andσ = 1 for differentk's. As we expect from previous discussion the effect of QCD and electroweak transitions is also apparent in this case for different peaks around modesk/2π =f = 1.6 × 10 −9 , 1.6 × 10 −8 , 1.6 × 10 −7 , and 1.6 × 10 −6 Hz. The value of scaled tensor power spectrum in for the mentioned modes is scaled to its value at a high frequencyf = 10 −4 Hz. We should emphasize again that this effect is due to the changes of the equation of state parameter and DoF during cosmic transitions (see fig. 1) in the calculation of eq. (14).
Non-Gaussianity at the CMB scale is constrained by Planck data [55]. However, there are less constraints on non-Gaussianity at small scales. To predict a realistic signal for the presence of non-Gaussianity we assume A G F 2 N L 1 [20]. This causes that the second term on the right hand side of eq. (25) is smaller than the first term which means that the non-Gaussian part of the curvature power spectrum is smaller than the Gaussian one. However, the contribution of non-Gaussianity on the second order PGW peak can be larger than the Gaussian part depending on the values of F N L andσ [20]. In fig. 6 the non-Gaussian induced PGW for [F 2 N L , A G ] = [10, 10 −2 ] and [60, 10 −4 ] for different modes is plotted using eqs. (23), (24), and (25). For these peaks pulsar timing array experiments or LISA can be sensitive enough to observe such peaks on the GW background which can be distinguished from other prediction, e.g. by first order phase transitions, due to the specific shape of the spectrum. Consequently, the effect of the thermal history of the Universe and the induced tensor spectra from scalar perturbations will be small at the peaks and can provide information about pre-BBN cosmology.
The shape of the tensor spectrum is shown in fig. 3 for different equation of state parameters for 0 ≤ kη ≤ 10. Moreover, the effect of the degrees of freedom of thermal bath particles evolving with temperature on the spectrum of induced PGWs at different frequencies is studied and shown by a correction function Z(k) (Z(T hc )) in eq. (23) and fig. 4. This function includes the evolution of the tensor perturbation and the retarded effects from the scalar Green's function and is particularly important around the QCD and electroweak transitions. It is possible that future experiments do not observe the first order PGW for a scale independent tensor spectral index and tiny values of the tensor to scalar ratio ∼ 10 −9 at first order, based on the current knowledge of ΛCDM cosmology and observations. Then it is expected that the induced PGW should be observable when GW experiments can probe the PGW relic down to ∼ 5 × 10 −23 (fig. 5). Larger values due to non-Gaussianity ( fig. 6) are measurable by the recently planned missions, and are distinguishable from astrophysical backgrounds. In a sense what is calculated here is the combined prediction of general relativity based on just the SM of particle physics using the current CMB observation for the early Universe cosmology. Principally, one should be able to observe the impact of thermal history of the early Universe on the induced PGW background from scalar perturbations. Since we have already observed the scalar perturbation at the CMB scale and their effect on structure formation, it would be reasonable to see its indirect impact on the production of induced PGW at smaller scales. Otherwise, there should be some theoretical model dependent explanation for it from alternative approaches to the big bang cosmology and inflation or an unknown effect at early epochs that has not been investigated before [8,10,12,22].
Considering the effect of the standard model DoF on the induced PGW will be influential for the future searches of gravitational waves with cosmological origin, since it helps to improve the theoretical prediction due to different effects for the GW spectra from the early Universe. In addition, studying details of the induced PGW affects our understanding about the impact of thermal evolution of the SM in the early Universe and pre-BBN cosmology on the GW background. Also, searching for the PGW can be useful as a cosmic laboratory for the indirect probe of new particles being present in the thermal bath of the early Universe, and new physics even without observing any boosted effects on the GW spectrum from the predicted cosmic phase transitions and nonstandard cosmological scenarios.