Constraints on primordial gravitational waves from the Cosmic Microwave Background

Searches for primordial gravitational waves have resulted in constraints in a large frequency range from a variety of sources. The standard Cosmic Microwave Background (CMB) technique is to parameterise the tensor power spectrum in terms of the tensor-to-scalar ratio, $r$, and spectral index, $n_{\rm t}$, and constrain these using measurements of the temperature and polarization power spectra. Another method, applicable to modes well inside the cosmological horizon at recombination, uses the shortwave approximation, under which gravitational waves behave as an effective neutrino species. In this paper we give model-independent CMB constraints on the energy density of gravitational waves, $\Omega_\text{gw} h^2$, for the entire range of observable frequencies. On large scales, $f \lesssim 10^{-16}\, \text{Hz}$, we reconstruct the initial tensor power spectrum in logarithmic frequency bins, finding maximal sensitivity for scales close to the horizon size at recombination. On small scales, $f \gtrsim10^{-15}\,\mbox{Hz}$, we use the shortwave approximation, finding $\Omega_\text{gw} h^2<1.7 \times10^{-6}$ for adiabatic initial conditions and $\Omega_\text{gw} h^2<2.8 \times10^{-7}$ for homogeneous initial conditions (both $2\sigma$ upper limits). For scales close to the horizon size at recombination, we use second-order perturbation theory to calculate the back-reaction from gravitational waves, finding $\Omega_\text{gw} h^2<1.0 \times10^{-6}$.


Introduction
Primordial gravitational waves (PGWs) offer a revelatory breakthrough for our knowledge of the physics of the early universe, but are currently unobserved. There are two possible sources for them: those produced during inflation, and those produced between the end of inflation and Big Bang Nucleosynthesis (BBN). For standard models of inflation, metric perturbations give rise to an almost scale invariant spectrum of PGWs, directly related to the energy scale of inflation. These result in a characteristic B-mode polarization signal in the Cosmic Microwave Background (CMB), which has been constrained by the Planck and BICEP2/Keck experiments [1,2].
A number of post inflationary mechanisms could also result in the production of PGWs (see e.g. [3] for a recent review). Some of these processes include: (1) During the reheating phase, the non-perturbative excitation of fields can result in a stochastic background of gravitational waves, with a well defined peak at f ∼ 10 7 − 10 8 Hz (see [4] for examples of GW production during preheating in a series of inflationary models); (2) If the curvature power spectrum has large, broad peaks on small scales, the production and merger of Primordial Black Holes (PBHs) can lead to a background within the range of direct detection experiments (for a recent review see [5]); (3) A network of cosmic strings (or other topological defects) can give rise to a background at lower frequencies. There are two sources from strings: the irreducible emission from the time-evolution of the energy momentum tensor during scaling, and the production and subsequent decay of cosmic string loops. There is still some debate in the literature as to the bounds placed on the string parameters, as they vary depending how the network evolution is modelled. However, the key property that is constrained is the combination of the string tension µ and Newton's constant G. For the most recent bounds arising from a search in the 10 2 Hz region for an isotropic stochastic background of GWs from the LIGO/VIRGO collaboration see [6]. Depending on the model, the bound varies between Gµ/c 2 ≤ 1.1 × 10 −6 to Gµ/c 2 ≤ 2.1 × 10 −14 . A different complementary set of bounds can be obtained using the pulsar timing limits which give for the two string models, Gµ/c 2 ≤ 1.6 × 10 −11 to Gµ/c 2 ≤ 6.2 × 10 −12 [7]. The wide variety of processes means it is important to constrain the energy density of gravitational waves, Ω gw h 2 , for the entire range of observable frequencies. Limits can be obtained using observations from BBN, pulsar timing, gravitational wave interferometers and the CMB.
Measurements of the CMB temperature and polarization power spectra can be used to constrain Ω gw h 2 . For single-field slow-roll models of inflation, the initial spectrum of tensor perturbations is well approximated by a power-law, parameterised by the tensor-to-scalar ratio, r, and spectral index, n t , and is directly related to Ω gw h 2 . This provides a limit on PGWs in the frequency range f 10 −16 Hz. At higher frequencies, the expected signal can be extrapolated, assuming the power-law is valid across many decades in scale. In this work we did not assume a power-law, instead directly reconstructing the tensor power spectrum in logarithmic frequency bins using the latest data from Planck and BICEP2/Keck. Primordial gravitational waves also have an effect on small-scale CMB anisotropies, through the back-reaction of tensor fluctuations on the cosmological background. The standard approach relies on the so-called shortwave approximation, which is valid for modes well inside the cosmological horizon [8][9][10][11][12][13]. Under this approximation, the energy-momentum tensor has an equation of state w = 1/3, and so acts as an effective relativistic neutrino species. In this paper we provide the most recent constraints using the shortwave approximation, considering both adiabatic and homogeneous initial conditions for the PGW perturbations.
For scales close to the horizon size at recombination, the shortwave approximation is no longer valid. No constraints on PGWs from the back-reaction of tensor fluctuations currently exist for these scales. In this work we use the approach of [14], which shows that the effective energy-momentum tensor of super-Hubble modes has the form of a fluid with equation of state w = −1/3. There is a calculable transition period between the w = −1/3 super-Hubble regime and the sub-Hubble w = 1/3 regime. Consequently a constraint can be found that isn't restricted to sub-Hubble gravitational waves.
The structure of the paper is as follows. In section 2 we review previous constraints on PGWs from CMB polarisation and from CMB temperature anisotropies. In section 3 we give model independent limits on the low-frequency reconstruction of the tensor power spectrum. Section 4 updates previous constraints that use the shortwave approximation using Planck 2018 data. In section 5 a constraint is found in the intermediate regime using an expression for the gravitational wave energy momentum tensor that is valid on all scales. 1 In section 6 we give conclusions.

Constraints from B-mode polarisation
Inflation is predicted to produce gravitational waves with both E and B-mode polarisation, but primordial density perturbations do not result in B-modes. Searches for inflationary gravitational waves have therefore focused on detecting B-mode polarisation of the CMB (see [15] and references therein). However, there are still significant contaminants to B-mode observations, such as gravitational lensing along the line of sight and galactic foregrounds, which need to be accurately modelled before constraints on cosmological parameters can be found.
When the scalar and tensor primordial power spectra have conventional power-law parameterisations, the tensor-to-scalar ratio r k is defined as the ratio of amplitudes evaluated at scale k. This is used preferentially to the tensor amplitude A t by convention, but because the scalar amplitude A s is well determined, by for example Planck [16], they can be interchanged easily. The current best constraint on the tensor-to-scalar ratio is r 0.002 < 0.056 at 95% confidence level when combining Planck 2018, BICEP2/Keck data and Baryon Acoustic Oscillations (BAO) [1,2].
A constraint on the tensor-to-scalar ratio r can be converted to the gravitational wave density parameter using eq. (140) of [3] which corrects the more commonly seen expression of [7]. The gravitational wave density parameter as a function of frequency is Here h contains the uncertainty in the Hubble parameter, H 0 = 100 h km s −1 Mpc −1 , Ω r is the density of relativistic species and k eq = √ 2H 0 Ω m / √ Ω r , for matter density Ω m , is the wavenumber of a mode that enters the horizon at matter-radiation equality.
For single-field slow-roll models of inflation the tensor primordial power spectrum is well approximated by where the standard value for the pivot scale, k * = 0.05 Mpc −1 and n t is the tensor spectral tilt. In slow-roll models there is a consistency relation, n t = −r/8, so inflation predicts a slightly red-tilted spectrum with −0.007 < n t < 0 (95% confidence). Conversions between frequencies and wavenumbers are done using, where the numerical factor comes from the speed of light and the definition of a parsec. Using the recent Planck and BICEP2/Keck constraint of r 0.002 < 0.056 , Ω gw h 2 ranges from 1.9 × 10 −13 to 7.4 × 10 −17 for frequencies 3.4 × 10 −19 Hz to 2.1 × 10 −17 Hz respectively.

Constraints from temperature anisotropies
In the shortwave approximation, primordial gravitational waves behave like massless neutrinos and therefore contribute to the effective number of relativistic degrees of freedom, N eff .
Assuming adiabatic initial conditions, Ω gw h 2 can be calculated directly from N eff as where it has been assumed that the standard model value N eff = 3.046 holds in the absence of primordial gravitational waves. The constant in (2.4) comes from the definition of N eff , where ρ r is the energy density of all relativistic species and ρ γ is the energy density of photons. Consequently, the conversion factor between the density parameter and the number of gravitational wave degrees of freedom is The combination of Planck 2018 + BAO data give a constraint of N eff = 2.99 +0.34 −0.33 (2σ) [16], which, using the upper limit, corresponds to Ω gw h 2 < 1.6 × 10 −6 (although a proper analysis should use Ω gw h 2 with a prior Ω gw h 2 ≥ 0).
In this approximation, the perturbations of gravitational waves are identical to those of massless neutrinos, and are described by their density, δ gw , velocity, θ gw , shear, σ gw , and higher-order moments, with fluid equations given by Eqs. (4.7a). It is possible that these could have non-adiabatic initial conditions, depending on the source of PGWs. Adiabatic initial conditions would be the sensible choice if primordial gravitational waves were a thermalised particle species produced by the decay of the inflaton, however most known sources of a cosmological gravitational wave background, including quantum fluctuations during inflation, reheating and cosmic strings, produce an unperturbed background (see [3]  Consequently the second choice of gravitational wave initial conditions, homogeneous initial conditions, have no initial density perturbation (in the Newtonian gauge). In this case the gravitational wave perturbations evolve differently to the neutrino perturbations and consequently the degeneracy between Ω gw and N eff is broken.
[18] details how CMB temperature anisotropies can be used to constrain short wavelength gravitational waves for both adiabatic and homogeneous initial conditions, using observations from WMAP (first-year), SDSS and the Lyman-α forest. Tighter constraints are seen for homogeneous gravitational waves compared to adiabatic gravitational waves by a factor of ≈ 5 − 10. These constraints have been updated using WMAP seven-year data, finding Ω gw h 2 < 8.7 × 10 −6 for adiabatic and Ω gw h 2 < 1.0 × 10 −6 (both 2σ) for homogeneous gravitational waves [19]. The adiabatic result has been obtained for Planck 2015 data [20][21][22], finding Ω gw h 2 < 1.7 × 10 −6 , but the homogeneous result has not been recently updated. All of these CMB constraints are smaller but are of the same order of magnitude as BBN constraints, but are valid to lower frequencies owing to the larger horizon size by the time of recombination.

Low frequencies 10 −16 Hz
Rather than assume a power law spectrum (2.2), we reconstruct P T (k) in logarithmic frequency bins, between a minimum and maximum wavenumber of log 10 [k Mpc] = −3.5 and −0.3 respectively, with an interval of ∆ log 10 [k Mpc] = 0.2. Outside of this range, the tensor transfer functions have very little sensitivity, so we set P T (k) to the lower and upper bin values.
We use Planck 2018 data [16] in combination with baryon-acoustic oscillation (BAO) data from the Baryon Acoustic Oscillation Survey (BOSS) DR12 [23], 6dF Galaxy Survey (6dFGS) [24] and Sloan Digital Sky Survey 'main galaxy sample' (SDSS-MGS) [25]. This corresponds to the TT,TE,EE + lowE + lensing + BAO data-set used in [16]. The precise Planck likelihoods used are the TT, TE and EE spectra at l ≥ 30, the low-likelihood using the Commander component separation algorithm [26] and the low−l EE likelihood from the SimAll algorithm in combination with Planck 2018 lensing.
Although Planck measured the CMB polarization over the full sky, the sensitivity for intermediate angular scales can be improved by using results from the BICEP2/Keck Array, with bandpowers in the range 20 < < 330. We use the most recent analysis from [2], which includes new data from the Keck array at 220 GHz.
For the low-frequency reconstruction, we assume an otherwise standard LCDM model, with adiabatic scalar perturbations parameterised by a power-law spectrum with scalar amplitude A s and spectral index n s . We assume three neutrinos species, two of these massless and a single massive neutrino with mass 0.06 eV. The other model parameters are the baryon density ω b ≡ Ω b h 2 , the cold dark matter density ω c ≡ Ω c h 2 , the Hubble parameter H 0 , and the optical depth to reionization τ . We assume flat priors on these parameters, and marginalise over the standard nuisance parameters in the Planck and BICEP2/Keck likelihood codes.
We perform Metropolis-Hastings Markov-chain Monte Carlo (MCMC) using a modified version of the Cobaya and camb codes [27]. 2 We run four MCMC chains, stopping them when the Gelman and Rubin R − 1 statistic is < 0.05. The 2σ upper limits for each of the 16 bins of P T (k) are then converted to Ω gw (k)h 2 using eq. (2.1) and are shown in figure 1. There is maximal sensitivity for scales scales close to the horizon size at recombination, and due to the decay of modes once they enter the horizon, these limits become much weaker for f 10 −16 Hz. A similar result was recently found in [28]. Above these frequencies, tighter constraints come from the second order result described in section 5. For comparison, we also plot an inflationary model with the Planck and BICEP2/Keck upper bound of r 0.002 < 0.056 and n t = −0.007.

High frequencies 10 −15 Hz
The effective energy-momentum tensor for gravitational waves in the short wavelength limit is [12,13] where straight lines denote covariant derivatives with respect to the background metric and h µν is the tensor perturbation to the conformal FLRW background metric, defined by, where τ is the conformal time and a(τ ) is the scale factor. The angled brackets . . . denote averaging over many wavelengths. We now illustrate that PGWs in this limit act like a massless neutrino species. CMB -low frequency Slow roll Planck/BK limit CMB -second order CMB -shortwave (adiabatic) CMB -shortwave (homogeneous) BBN -17 -15 -13 -11 -9 log 10 (f/Hz) Figure 1. CMB constraints on the gravitational wave density, Ω gw h 2 , as a function of the wavenumber, k, and frequency, f . The filled green bars show constraints from a reconstruction of P T (k) in logarithmic frequency bins, and the dashed green line shows a slow-roll model with the Planck and BICEP2/Keck upper bound of r 0.002 = 0.056 and n t = −0.007. Shortwave constraints for adiabatic and homogeneous initial conditions are shown in red and orange (solid) respectively. The second-order back-reaction result, which we apply to scales between k = 0.1 Mpc −1 and k = 1 Mpc −1 is shown in blue (dot-dash). A constraint from BBN [29] (black dashed) is also shown for comparison. The second-order back-reaction and shortwave results are integrated constraints across the given frequency range.
The shortwave approximation (SWA) states that the perturbation is well inside the horizon and that it is oscillating much faster than the background time-scale, so we can assume that the background is approximately flat. Consequently, covariant derivatives become partial derivatives; Furthermore, the equation of motion becomes Considering plane waves along z the solutions of the equation of motion depend on the retarded time τ −z and consequently spatial and temporal derivatives are equivalent. Therefore the trace of the energy-momentum tensor which implies that the equation of state w gw = 1/3, as for massless neutrinos. The full solution is a sum of plane waves of positive and negative frequencies along all three spatial directions (such that the full solution is isotropic), but this argument applies to each of these spatial directions separately.

Initial conditions
The calculation of adiabatic initial conditions for a universe containing photons, neutrinos, cold dark matter and baryons with linear perturbations is detailed in [30] (section 7). The initial conditions for non-adiabatic modes were calculated in [31]. The general methodology is to solve the perturbed Einstein and conservation equations in series solutions for small kτ , where k is the wavenumber of the mode under consideration and τ is the conformal time.
Matching the coefficients of the expansion gives the initial conditions for the full solution of the set of differential equations. This was performed in the synchronous gauge for the above system with the addition of gravitational waves with density (δ gw ), velocity (θ gw ), and shear (σ gw ) perturbations. The gravitational wave perturbations are expanded as for the other species, for example the gravitational wave density perturbation is expanded as where the coefficient, a δ n , along with the coefficients for the other species and perturbations, are the quantities we want to determine to establish the early time behaviour. We have four equations from the Einstein equations and a set of fluid conservation equations for each species [30]. For SWA gravitational waves these fluid equations arė where h and η are the synchronous gauge metric perturbations and dots denote differentiation with respect to conformal time.
As mentioned previously the adiabatic mode has gravitational wave perturbations identical to the neutrino perturbations. For the homogeneous mode there is one free parameter in the set of coefficients which is fixed by transforming to the Newtonian gauge using the well-known transformation relations (see [30]) and enforcing the condition that the zeroth order gravitational wave density coefficient is zero,ã δ 0 = 0 (where the tilde denotes that this condition is imposed in the Newtonian gauge). The resulting adiabatic, homogeneous and gravitational wave isocurvature modes are shown in table 1. The fractional contribution to the radiation density R i = ρ i / j ρ j for i, j = γ, ν, gw. The perturbations have their standard definitions (and have a tilde in the Newtonian gauge), as do the metric perturbations. The behaviour of gravitational wave perturbations for the well-known adiabatic mode and neutrino density isocurvature mode are given along with the homogeneous gravitational wave mode and two new modes, the gravitational wave velocity isocurvature mode and the gravitational wave shear isocurvature mode.
The homogeneous mode calculated here differs from the one quoted in [32] and consequently the one used in [18]. This is for two main reasons; firstly, the gravitational wave density perturbations are not assumed to be sub-dominant to those for photons and neutrinos and secondly, and more importantly, because the homogeneous mode quoted in [32] is a linear sum of the homogeneous mode given here and the neutrino density isocurvature mode. In the homogeneous mode of [32] the photon and neutrino density perturbations are assumed equal. We make no such assumption but can obtain the same mode by combining the homogeneous and neutrino density isocurvature modes given in table 1. The combination of these two modes can be done in the initial condition correlation matrix [31] to give an equivalent result but it is the homogeneous mode given here that is the true independent mode for gravitational waves. Because of this, small differences between the results of [18] and this analysis should be expected for the homogeneous mode, even for the same data.
To calculate the CMB power spectrum we modified the camb code to include the gravitational wave equations of motion, by duplicating the massless neutrino equations, and setting the initial conditions according to table 1. The changes in the CMB power spectrum when including adiabatic or homogeneous gravitational waves are shown in figure 2. The contributions from the Sachs-Wolfe effect (SW) [33], integrated Sachs-Wolfe effect (ISW) and Doppler shift (DOP) are shown separately, along with the cross-correlations between them. The most noticeable difference between the adiabatic and homogeneous cases is that homogeneous gravitational waves decrease the total power whereas adiabatic gravitational waves increase the total power. This is because the homogeneous gravitational wave and photon perturbations are out of phase with each other when inside the horizon due to the initial conditions having opposite signs (see table 1).
Including adiabatic gravitational waves has a very small effect on the Doppler term. Most of the change in the first peak of the power spectrum is due to the SW effect and the SW-ISW cross-correlation. This is also true for the homogeneous mode with the addition of a large change in the low-part of the spectrum, driven by the SW, SW-ISW and Doppler terms. The enhancement of the first peak and the decrease at low-is a background effect (i.e. still observable when the gravitational wave perturbations are turned off).
It is also worth noting the behaviour of the new modes. For the gravitational wave shear isocurvature mode the neutrino and gravitational wave shear both have zero order initial conditions but they balance in such a way that the right-hand side of the relevant perturbed Einstein equation (∝ i R i σ i ) is zero. This is also true for the density and velocity perturbations such that there are no initial metric perturbations (in synchronous or Newtonian gauge). Because of this all perturbations other than those for neutrinos and gravitational waves stay zero for all times. Using the line-of-sight integral approach of [34] it is clear that there will be no contribution to the CMB power spectrum from the shear mode.
The gravitational wave velocity isocurvature mode is the analogue of the neutrino velocity isocurvature mode and consequently behaves very similarly. In the Newtonian gauge densities and potentials have terms that go as 1/(kτ ). This is a consequence of the Newtonian gauge being inadequate when there is a non-zero anisotropic stress and does not mean that the perturbations diverge as kτ → 0 [31]. Similarly, there is a gravitational wave density isocurvature mode not shown in table 1 that is a direct analogue of the neutrino density isocurvature mode.
Isocurvature modes are well constrained by current CMB observations such that we will only consider constraints to the adiabatic and homogeneous modes here [1]. The baryon and cold dark matter perturbations and modes are not included in table 1 for notational convenience and due to the matter perturbations being calculable from the perturbations of the relativistic species [31].

Adiabatic Homogeneous
Neut. Dens. IC GW Vel. IC GW. Shear IC  Table 1. Initial conditions on synchronous gauge (top) and Newtonian gauge quantities (bottomwith tildes) to second order in kτ (in the synchronous gauge) for the modes relevant to the gravitational wave (GW) analysis. The adiabatic and neutrino density isocurvature (IC) modes are well known and are extended here to include gravitational wave initial conditions. The homogeneous mode is the same as the homogeneous mode of [32], except for a decoupling of the neutrino density isocurvature mode, while the gravitational wave velocity IC and gravitational wave shear IC modes are new. There is a gravitational wave density IC mode that is a rescaling of the neutrino density IC mode and is not shown.

Parameter constraints
To obtain limits on the density of gravitational waves our modified version of camb [35] was integrated into the cosmological parameter estimation code CosmoMC [36] to perform an MCMC analysis.
We use the same data as in section 3, but do not include tensor modes and hence BICEP2/Keck data. The base cosmology is an otherwise standard LCDM model, with the addition of Ω gw h 2 . We obtain the following 95% upper limits on the gravitational wave density parameter; Ω gw h 2 < 1.7 × 10 −6 (Shortwave, adiabatic) , (4.8) Ω gw h 2 < 2.8 × 10 −7 (Shortwave, homogeneous) . (4.9) These constraints can be seen in figure 1, along with the Big Bang Nucleosynthesis (BBN) constraint from [29], the low-frequency constraint of section 3 and the intermediate frequency constraint of section 5. The CMB constraints extend to much lower frequencies than those from BBN -the exact range of validity for the shortwave approximation is discussed in section 5.2. Note that these results, in contrast to the direct reconstruction, are integrated constraints across the range of frequencies.
5 Intermediate frequencies 10 −15 Hz f 10 −16 Hz In order to consider gravitational waves without the restriction of the shortwave approximation (SWA), the work of [14,37] is closely followed. Here the effective density and pressure of PGWs are calculated using the second order back-reaction of the tensor fluctuations on the metric. The second order part changes the zeroth order (or background) Einstein equations, modifying the Friedmann and continuity equations. Expanding the Einstein equations to second order and averaging over all space (the spatial average of the linear terms are zero by definition), whereG µ ν andT µ ν are the background Einstein and energy-momentum tensors respectively, and δ (2) G µ ν and δ (2) T µ ν are the second order perturbations to the Einstein and energymomentum tensor respectively. This allows an effective energy-momentum tensor, τ µ ν to be defined, 2) whose components are τ µ ν = diag −ρ gw , p gw , p gw , p gw . The off-diagonal terms are zero under averaging, assuming an isotropic source of PGWs.
In general the choice of gauge is important when calculating the effective energymomentum tensor, but since the tensor perturbation h ij in the transverse-traceless gauge defined in eq. (4.2) is gauge-invariant, this will not be a problem (see [14] for details). Consequently, in vacuum, the evaluation of the effective energy-momentum tensor simplifies to evaluating the perturbed Einstein tensor, This can be done using eq. (35.58b) of [13]. eq. (5.3) is very similar to eq. (35.59b) of [13], where R (B) µν and R (2) µν (h) are the background and second-order Ricci tensors respectively. However, in this analysis the average is a spatial average, not over many wavelengths, so can be applied to super-horizon modes.
For the transverse-traceless perturbation defined in eq. (4.2) the components of the perturbed Einstein tensor are, where we have introduced H =ȧ/a, withȧ ≡ da/dτ . We note that these agree with Eqs. (8-11) of [37] except for the last two terms of eq. 5.7. It does however agree with eq. (2.30) of [38]. This could be due to an ambiguity in the notation, to be clear, byḣ km,i we mean ∂ i ∂ 0 h km .
Because the usage of the final expressions from this section are integral to the analysis of this section of this paper the following calculation is given in detail. From the expressions for the Einstein tensor the density and pressure are found by spatial averaging defined via [14,37], and using the gravitational wave equation of motion to find, where w (0) is the equation of state of the background spacetime, and we have also performed an average Q over the ensemble average of stochastic initial conditions. These satisfy the continuity equationρ The background equation of state appearing in eq. (5.10b) is related to an important assumption. The approach of [14] inherently assumes that the back-reaction is a small perturbation to the background spacetime. Consequently ρ gw /ρ crit 1 is required at all times.
We next Fourier transform the tensor metric perturbation, and decompose the Fourier components in terms of the polarisation tensor ij [39], where λ is the gravitational wave helicity andh is the gravitational wave amplitude. In Fourier space the spatial averaging of a product of two general functions, where we have rewritten the complex exponentials from the Fourier transforms as a Dirac delta function and used this to do one of the wavenumber integrals. Here we have also set V = 1 as it is only included to keep track of dimensions (see [40] chapter 3). We can do this integral in spherical polar coordinates to find, where we have evaluated the products of the polarisation tensor using [39], Using this for the density and pressure of gravitational waves, We can separate the initial condition from the time evolution of the gravitational wave amplitude as,h (k, τ ) = A k D(k, τ ) , (5.19) such that the primordial power spectrum, specifies the ensemble average of stochastic initial conditions.
The time evolution of the gravitational wave amplitude is now described by the function D(k, τ ) which obeys the gravitational wave equation of motion, which comes from the first order vacuum Einstein equation for gravitational waves R (1) µν (h) = 0. The helicity dependence ofh(k, τ ) was dropped earlier because this equation of motion is helicity independent.
We now have our final expressions for the density and pressure, and the primordial power spectrum is conventionally parameterised via, Consequently the methodology for calculating the density and pressure of gravitational waves is as follows. We solve the gravitational wave equation of motion (eq. (5.21)) for a given background cosmology and use the solution for D(k, τ ) to evaluate the k-space density and pressure using Eqs. (5.23). This is then integrated with the power spectrum of eq. (5.24) to get the total homogeneous density and pressure according to Eqs. (5.22).

Equation of state for single fluid backgrounds
The k-space density and pressure can be evaluated analytically for radiation, matter and de Sitter backgrounds. This can be used to consider the behaviour of the k-dependent equation of state before integrating over k to find the total gravitational wave density and pressure. For convenience we will define x = kτ . For a radiation background, D(x) = B sinc x, where B is the initial condition on the gravitational wave amplitude,ȧ/a = 1/τ and the equation of state of the background is 1/3. Consequently,ρ gw = B 2 4x 2 τ 2 − 7 + 2x 2 + 7 cos 2x + 6x sin 2x , (5.25) p gw = B 2 12x 2 τ 2 3(7 − 4x 2 ) cos 2x + 26x sin 2x + 2x 2 − 21 .
Finally, for a de Sitter background, a/a = −1/τ and the equation of state of the background is −1. Considering only the even part of D(x) such that C = 0 (though the conclusions are the same if C is included), For the above backgrounds we can calculate the equation of state parameter w gw (k, τ ) = p gw (k, τ )/ρ gw (k, τ ) in the super-Hubble regime (x 1) by expanding in x and in the sub-Hubble regime (x 1) by averaging trigonometric functions, e.g., sin 2x = cos 2x = 0. Doing this we find for radiation, matter and de Sitter backgrounds. The equation of state for a general ΛCDM background can be solved numerically, and is shown in figure 3 for Planck 2018 parameter values and k = 0.05 Mpc −1 . It starts at −1/3 when the mode is outside the horizon, then goes through a transition period where it goes through large negative and positive values before exhibiting stable oscillations about an average value of w gw = 1/3.

Shortwave validity
The SWA is valid on scales averaged over a "sufficient number of wavelengths". With a prescription how to evolve PGWs on super-horizon scales, we now quantify how many wavelengths are required to give accurate enough CMB constraints, and hence what frequency range the SWA can be applied on.
From figure 3 it is clear that there are two time-scales in the evolution of w gw : the high-frequency oscillatory behaviour, and the lower frequency transition from −1/3 to 1/3. Our task is to find the number of wavelengths required such that a constant w gw = 1/3 is a good approximation to calculate the CMB power spectrum.
A cosmic variance limited CMB experiment, up to a given , requires a precision of approximately [41] δC C = 3 .  in the power spectrum C . This corresponds to 0.1-0.2% for = 2000. To proceed we assume the PGW source is a δ-function for a given frequency, with energy density Ω gw h 2 . We then compute the fractional difference in C for the correct evolution, compared to the SWA with constant w gw = 1/3. This error decreases for higher frequency sources, as these enter the horizon and thus behave like a relativistic species earlier. We determine the minimum frequency such that δC /C < 0.2% for all . Of course, it is less important to have the correct evolution the smaller Ω gw h 2 is. We therefore choose Ω gw h 2 = 5.6 × 10 −6 , corresponding to ∆N eff = 1 in the SWA, such that the PGW source has an appreciable affect on the background evolution at matter-radiation equality.
We find that a source with k > 1 Mpc −1 is required for the SWA to satisfy δC /C < 0.2% for all . This corresponds to the mode undergoing 50 oscillations by the epoch of equality. This limit is indicated in figure 3 for a mode with a smaller k = 0.05 Mpc −1clearly by equality a constant w gw = 1/3 is not yet a good approximation. The SWA result in figure 1 therefore extends from k > 1 Mpc −1 . We note that, compared to figure 2 of [18], they use a factor of ∼ 20 wavelengths. Our analysis suggests that a slightly more conservative limit is required.

Behaviour of gravitational wave density and pressure
The gravitational wave density and equation of state exhibit a range of interesting physical behaviours. Figure 4 shows these for standard ΛCDM parameter values as a function of k and t. A smoothing has been applied to w gw to more clearly show the behaviour when the gravitational wave amplitude is highly oscillatory. The super-horizon and sub-horizon regimes can be seen clearly, along with the transition region between the two.
When super-horizon the gravitational wave equation of state is −1/3 as verified above, apart from at late times, during the matter to cosmological constant transition, when it goes below −1/3 (this can be seen by closely inspecting the top-left of the lower panel of figure  4). This is a phenomenon that has not previously been mentioned in the literature. We analytically verified the behaviour by solving the equation of motion in a matter-cosmological constant background for small k and matched solutions for different time regimes (see appendix for details). This showed that the equation of state of super-horizon modes dips below −1/3 during the matter to cosmological constant transition to values of ∼ −0.5 (dependent on various parameters) but returns back to −1/3 soon after the cosmological constant is dominating.
The integrated density of eq. (5.22) and the subsequent equation of state are shown in figure 5, for a representative PGW source with n t = 3, k min = 0.1 Mpc −1 and k max = 1 Mpc −1 . The lower cutoff is chosen to be compatible with the low-frequency constraint, and the spectral index must be relatively steep, n t 3, to also satisfy this constraint. The high-frequency cutoff is chosen as the SWA can be used for frequencies above this. The suband super-Hubble regimes are clear in both cases and the transition region between the two can also be seen.
We note that the energy density is negative for super-Hubble modes, as stated in section 4 of [37]. Since w gw = −1/3 it can be interpreted as additional positive curvature. This contribution can lead to a reduction in the expansion rate, depending on the integration limits and spectrum in eq. (5.22).

Perturbations
The effective energy-momentum tensor at linear order is given by where Π i j = τ i j − δ i j τ k k /3 is the anisotropic stress. Previously we have calculated the background energy density and pressure. The fluctuating part can be calculated by subtracting the average, ∆ µ ν = τ µ ν − τ µ ν . Numerically, however, these are much more challenging to calculate, as they cannot easily be written in terms of the initial spectrum of fluctuations. We therefore take a phenomenological approach to the PGW perturbations, treating them as an effective Parameterized Post-Friedmann (PPF) fluid.  The density has two regimes, one where it goes as a −2 with a negative density (dashed) and one where it goes as a −4 with a positive density (solid), with a transition in between. These regimes can be seen more clearly in the equation of state. There is some numerical noise as the equation of state approaches 1/3 but this has no observable consequences.
The PPF framework is usually used in 'smooth' dark energy models [42], but has several properties useful to model PGW perturbations. Firstly, it is able cross the w = −1 divide, which occurs for PGW oscillations after entering the horizon. Secondly, it is designed to conserve energy and momentum on large scales, where PGWs behave like positive curvature with w gw = −1/3. Finally, on small scales the PPF fluid is designed to be smooth compared to cold dark matter, which one would expect for PGWs due to the pressure support with w gw = 1/3. In our approach we use the default PPF parameters in camb, and leave a more detailed study of PGW perturbations to future work.
One point worth mentioning is that in the SWA, PGWs are modelled as an effective neutrino species, with a hierarchy of moments describing the perturbations. In the secondorder method no such hierarchy exists, and the fluid is described by its effective energy- momentum tensor. There are therefore two main observable differences expected compared to the shortwave treatment; (1) due to the background evolution, and (2) in the treatment of perturbations.

Observables
The effects of high-frequency gravitational waves (as in section 4) on cosmological observables can be compared to the intermediate gravitational wave analysis of this section. Figure 6 shows the CMB power spectrum for Ω gw h 2 = 5.6 × 10 −6 for SWA gravitational waves with adiabatic and homogeneous initial conditions. The intermediate gravitational waves are shown for the same density with representative parameters of n t = 3 and k min , k max = 0.1, 1 Mpc −1 . The intermediate analysis changes the temperature anisotropies in similar ways to the adiabatic SWA gravitational waves.
The fractional changes in the Hubble rate, H(z), and the scale of the sound horizon, r s , are shown in figure 7. As expected, due to these quantities only depending on the background and not the perturbations, the adiabatic and homogeneous high-frequency gravitational waves have identical affects on these parameters. The intermediate gravitational waves increase the Hubble rate similarly to the SWA result when dominated by high-frequency, w gw = 1/3 modes but decreases the Hubble rate at high redshift when dominated by w gw = −1/3 modes. The same affect is seen in the scale of the sound horizon but with opposite sign and a small move to lower redshift.
We note that the intermediate model shares some similarities with the axion-model that can potentially alleviate the Hubble tension [43,44]. In particular, there is an early dark energy (EDE) phase with w gw = −1/3 before a radiation phase with w gw = 1/3. Even though the energy density is negative when w gw = −1/3, the sound horizon can still be reduced at the time of recombination. We leave the study of whether the PGW model can reduce the Hubble tension to future work. This analysis assumes that the gravitational wave density is small enough that it can be calculated as a perturbation on a ΛCDM background. This was tested by iteratively recalculating the background including gravitational waves. Repeating this procedure until convergence shows an error of less than 0.01% in H(z) over all z, for the maximum value of Ω gw h 2 allowed by data. We conclude that this is a small enough error to use the approximation that gravitational wave back-reaction can be calculated on a standard ΛCDM background.

Parameter constraints
To obtain limits on Ω gw h 2 , a modified version of camb was integrated into cobaya to perform an MCMC analysis. We use the same data as in section 4, using an otherwise standard ΛCDM model. For PGWs, we choose k min = 0.1 Mpc −1 , as below this the low-frequency constraint dominates, and k min = 1 Mpc −1 , as above this the SWA can be used. We marginalise over the tensor spectral index, n t , in the prior range 3 to 5, where the lower limit is chosen to be compatible with the low-frequency constraint. The upper limit is chosen so as to include a range of short lasting early universe phenomena. Any production mechanism that produces gravitational waves in a short time frame will correspond to a large tilt. As examples, both cosmic strings and first order phase transitions can produce gravitational waves in the low to intermediate frequency regime [3,[45][46][47].
We obtain the following 95% upper limits on the gravitational wave density parameter; Ω gw h 2 < 1.0 × 10 −6 (Second-order) , (5.38) and plot this constraint in figure 1. We note this is similar to the shortwave adiabatic result.

Conclusions
In this paper we have presented constraints on primordial gravitational waves from the CMB for the entire range of observable frequencies. This includes updated constraints from B-mode polarisation at the lowest frequencies, the shortwave approximation at high frequencies, and a new intermediate constraint that bridges the region of applicability of the two. These constraints are compatible at their extremities and provide the tightest current constraints in particular frequency ranges.
The constraint from low polarisation shows that peak sensitivity occurs for scales close to the horizon size at recombination, corresponding to f ∼ 10 −17 Hz, with a gravitational wave density Ω gw h 2 ∼ 10 −16 . These limits become much weaker for f 10 −16 Hz, and at f ∼ 3×10 −16 Hz a stronger result comes from the second-order back-reaction of gravitational waves. This allows us to place a limit of Ω gw h 2 < 1.0 × 10 −6 (at 95% confidence), in a previously unconstrained frequency region of 10 −15 Hz f 3 × 10 −16 Hz. At higher frequencies, f 10 −15 Hz, we use the shortwave approximation (SWA) to update previous constraints and quantify the validity of the SWA, finding Ω gw h 2 < 1.7 × 10 −6 for adiabatic initials conditions and Ω gw h 2 < 2.8 × 10 −7 for homogeneous initial conditions (both at 95% confidence).
These constraints will be tightened by future ground and space based CMB observations from CMB-S4 and from polarisation via. LiteBIRD, CORE and PIXIE among others [48][49][50][51]. These will result in an order of magnitude improvement in the measurement of extra relativistic degrees of freedom and an even greater improvement in the tensor-to-scalar ratio. Combining these with other cosmological observables promises to further illuminate the early Universe.
There are several possibilities for future work. Due to the numerical challenges of calculating the fluctuations due to the second-order back-reaction, in this analysis we have treated them as an effective PPF fluid. In future work we plan to extend the line-of-sight CMB formalism to calculate these. It is worth noting though that, even in the shortwave limit, differences are expected compared to modelling them as an effective neutrino species with a hierarchy of moments. One further avenue might be investigating the possibility of PGWs alleviating the Hubble tension. The second-order model, with an appropriate source of PGWs, increases the relativistic degrees of freedom at recombination, thereby reducing the sound horizon, and having an early dark energy phase with w gw = −1/3. This would, however, require a non-standard source with a steep n t 3 spectrum peaking at f ∼ 10 −15 Hz. The general solution is D 0 (x) =D 0 − α coth x. Imposing that the gravitational wave amplitude is finite as x → 0, This constant is going to be set to 1 in most cases.

Perturbed solutions
The equation of motion to order κ 2 is, This can be solved in three separate regimes, low-x, intermediate-x and high-x and these solutions can be matched together at x a and x b .

Low-x solution
For small x the equation of motion is, with solution, where the initial condition is D 1 (0) = 0.

Intermediate solution
The where erf x is the error function, erfi x is the imaginary error function, p F q is the generalised hypergeometric function and the matching onto the low-x solution at x a determines the coefficients C 1 and C 2 .

High-x solution
For large x the equation of motion becomes,

GW equation of state parameter
The gravitational wave equation of state parameter for a ΛCDM background from the above analytics and from a numerical computation can be seen in figure 8. The analytic solutions were matched at x a = 0.35 and x b = 1.15 to get the best agreement with the numerics. They confirm the fact that the equation of state departs from −1/3 for super-horizon GWs during the matter-cosmological constant transition but return back to −1/3 when the cosmological constant comes to dominate. There are discontinuities due to imperfect matching of the solutions. Effectively the solutions are not of high enough order to fully encompass the behaviour in their specific regimes. This could be improved by using the intermediate solution twice and having four separate matched regimes but the analysis given here is sufficient to verify the numerical behaviour.