Stochastic backgrounds of relic gravitons: a theoretical appraisal

Stochastic backgrounds or relic gravitons, if ever detected, will constitute a prima facie evidence of physical processes taking place during the earliest stages of the evolution of the plasma. The essentials of the stochastic backgrounds of relic gravitons are hereby introduced and reviewed. The pivotal observables customarily employed to infer the properties of the relic gravitons are discussed both in the framework of the $\Lambda$CDM paradigm as well as in neighboring contexts. The complementarity between experiments measuring the polarization of the Cosmic Microwave Background (such as, for instance, WMAP, Capmap, Quad, Cbi, just to mention a few) and wide band interferometers (e.g. Virgo, Ligo, Geo, Tama) is emphasized. While the analysis of the microwave sky strongly constrains the low-frequency tail of the relic graviton spectrum, wide-band detectors are sensitive to much higher frequencies where the spectral energy density depends chiefly upon the (poorly known) rate of post-inflationary expansion.

1 The spectrum of the relic gravitons 1

.1 The frequencies and wavelengths of relic gravitons
Terrestrial and satellite observations, scrutinizing the properties of the electromagnetic spectrum, are unable to test directly the evolution of the background geometry prior to photon decoupling. The redshift probed by Cosmic Microwave Background (CMB in what follows) observations is of the order of z dec ≃ 1087 and it roughly corresponds to the peak of the visibility function, i.e. when most of the CMB photons last scattered free electrons (and protons). After decoupling the ionization fraction drops; the photons follow null geodesics whose slight inhomogeneities can be directly connected with the fluctuations of the spatial curvature present before matter-radiation decoupling.
The temperature of CMB photons is, today, of the order of 2.725 K. The same temperature at photon decoupling 2 must have been of the order of about 2962 K, i.e. 0.25 eV. The CMB temperature increases linearly with the redshift: this fact may be tested empirically by observing at high redshifts clouds of chemical compounds (like CN) whose excited levels may be populated thanks to the higher value of the CMB temperature [1,2].
The initial conditions for the processes leading to formation of CMB anisotropies are set well before matter radiation equality and right after neutrino decoupling (taking place for temperatures of the order of the MeV) whose associated redshift is around 10 10 . The present knowledge of particle interactions up to energy scales of the order of 200 GeV certainly provides important (but still indirect) clues on the composition of the plasma.
If ever detected, relic gravitons might provide direct informations on the evolution of the Hubble rate for much higher redshifts. In a rudimentary realization of the ΛCDM paradigm 3 , the inflationary phase can be modeled in terms of the expanding branch of de-Sitter space. Assuming that, right after inflation, the Universe evolves adiabatically and is dominated by radiation, the redshift associated with the end of inflation can be approximately computed as z + 1 ≃ 3 × 10 28 T γ0 2.725 K where g ρ denotes the weighted 4 number of relativistic degrees of freedom at the onset of the radiation dominated evolution and 106.75 corresponds to the value of the standard model of particle interactions. In Eq. (1.1) it has been also assumed, quite generally, that H ≃ 10 −5 M P as implied by the CMB observations. For more accurate estimates the quasi-de Sitter nature of the inflationary expansion must be taken into account. In the ΛCDM paradigm, the basic mechanism responsible for the production of relic gravitons is the parametric excitation of the (tensor) modes of the background geometry and it is controlled by the rate of variation of space-time curvature.
In the present article the ΛCDM paradigm will always be assumed as a starting point for any supplementary considerations. The reasons for this choice are also practical since the experimental results must always be stated and presented in terms of a given reference model. Having said this, 2 Natural unitsh = c = k B = 1 will be adopted. In this system, K = 8.617 × 10 −5 eV. 3 In the acronym ΛCDM, Λ qualifies the dark-energy component while CDM qualifies the dark matter component. 4 Fermions and bosons contribute with different factors to g ρ . By assuming that all the species of the standard model are in local thermodynamic equilibrium (for instance for temperature higher than the top quark mass), g ρ will be given by g ρ = 28 + (7/8)90 = 106.75 where 28 and 90 count, respectively, the bosonic and the fermionic contributions. most of the considerations presented here can also be translated (with the appropriate computational effort) to different models.
Given a specific scenario for the evolution of the Universe (like the ΛCDM model), the relic graviton spectra can be computed. The amplitude of the relic graviton spectrum over different frequencies depends upon the specific evolution of the Hubble rate. The theoretical error on the amplitude increases with the frequency: it is more uncertain (even within a specified scenario) at high frequencies rather than at small frequencies.
The experimental data, at the moment, do not allow either to rule in or to rule out the presence of a primordial spectrum of relic gravitons compatible with the ΛCDM scenario. The typical frequency probed by CMB experiments is of the order of 5 ν p = k p /(2π) ≃ 10 −18 Hz = 1 aHz where k p is the pivot frequency at which the tensor power spectra are assigned. CMB experiments will presumably set stronger bounds on the putative presence of a tensor background for frequencies O(aHz). This bound will be significant also for higher frequencies only if the whole post-inflationary thermal history is assumed to be known and specified. The typical frequency window of wide-band interferometers and of the electromagnetic spectrum (plot at the right). In both plots the common logarithms of the (comoving) frequency and of the (comoving) wavelengths are reported, respectively, on the horizontal and on the vertical axis.
(such as Ligo and Virgo) is located between few Hz and 10 kHz, i.e. roughly speaking, 20 orders of magnitude larger than the frequency probed by CMB experiments. The frequency range of wide-band interferometers will be conventionally denoted by ν LV . To compute the relic graviton spectrum over the latter range of frequencies, the evolution of our Universe should be known over a broad range of redshifts. We do have some plausible guesses on the evolution of the plasma from the epoch of neutrino decoupling down to the epoch of photon decoupling. The latter range of redshifts corresponds to an interval of comoving frequencies going from ν p ≃ 10 −18 Hz up to ν bbn ≃ 0.01 nHz ∼ 10 −11 Hz (at most).
In Fig. 1 (plot at the left) the frequency range of the relic graviton spectrum is illustrated, starting from the (comoving) frequency ν p whose associated (comoving) wavelength is of the order of of 10 26 m, i.e. roughly comparable with the Hubble radius at the present time. The low-frequency branch of the spectrum can be conventionally defined between ν p and ν bbn . The largest frequency of the relic graviton spectrum (i.e. ν max ) is of the order of 0.1 GHz in the ΛCDM scenario. Thus the highfrequency branch of the graviton spectrum can be conventionally defined for ν bbn < ν < ν max . In summary we can therefore say that • the range ν p < ν < ν bbn will be generically referred to as the low-frequency domain; in this range the spectrum of relic gravitons basically follows from the minimal ΛCDM paradigm; • the range ν bbn < ν < ν max will be generically referred to as the high-frequency domain; in this range the spectrum of relic gravitons is more uncertain.
The high-frequency branch of the relic graviton spectrum, overlapping with the frequency window of wide-band detectors (see shaded box in the left plot of Fig. 1), is rather sensitive to the thermodynamic history of the plasma after inflation as well as, for instance, to the specific features of the underlying gravity theory at small scales. This is why we said that the theoretical error in the calculation of the relevant observables increases, so to speak, with the frequency.
In Fig. 1 (plot at the right) the electromagnetic spectrum is reported in its salient features. It seems instructive to draw a simple minded parallel between the electromagnetic spectrum and the spectrum of relic gravitons. Consider first the spectrum of relic gravitons (see Fig. 1, plot at the left): between 10 −18 Hz (corresponding to ν p ) and 10 kHz (corresponding to ν LV ) there are, roughly, 22 decades in frequency. A similar frequency gap (see Fig. 1, plot at the right), if applied to the well known electromagnetic spectrum, would drive us from low-frequency radio waves up to x-rays or γ-rays. As the physics explored by radio waves is very different from the physics probed by γ rays, it can be argued that the informations carried by low and high frequency gravitons originate from two very different physical regimes of the theory. Low frequency gravitons are sensitive to the large scale features of the given cosmological model and of the underlying theory of gravity. High frequency gravitons are sensitive to the small scale features of a given cosmological model and of the underlying theory of gravity.
The interplay between long wavelength gravitons and CMB experiments will be specifically discussed in the subsection 1.2. The main message will be that, according to current CMB experiments, long wavelength gravitons have not been observed yet. The latter occurrence imposes a very important constraint on the low frequency branch of the relic graviton spectrum of the ΛCDM scenario whose salient predictions will be introduced in subsection 1.3. According to the minimal ΛCDM paradigm a very peculiar conclusion seems to pop up: the CMB constraints on the low-frequency tail of the graviton spectrum jeopardize the possibility of any detectable signal for frequencies comparable with the window explored by wide band interferometers (see subsection 1.4). The natural question arising at this point is rather simple: is it possible to have a quasi-flat low-frequency branch of the relic graviton spectrum and a sharply increasing spectral energy density at high-frequencies? This kind of signal is typical of a class of completions of the ΛCDM paradigm which have been recently dubbed TΛCDM (for tensor -ΛCDM). The main predictions of these models will be introduced in subsection 1.5. We shall conclude this introductory section with a discussion of two relevant constraints which should be applied to relic graviton backgrounds in general, i.e. the millisecond pulsar and the big-bang nucleosynthesis constraint (see subsection 1.6).

Long wavelength gravitons and CMB experiments
The bounds on the backgrounds of relic gravitons stemming from CMB experiments are phrased in terms of r T which is the ratio between the tensor and the scalar power spectra at the same conventional scale (often called pivot scale). While the use of r T is practical (see e. g. the 5-year WMAP data [3,4,5,6,7]), it assumes the ΛCDM scenario insofar as the curvature perturbations are adiabatic. Within the ΛCDM model, the tensor and scalar power spectra can be parametrized as A R = (2.41 ± 0.11) × 10 −9 , k p = 0.002 Mpc −1 .
where k p is the pivot wave-number and A R is the amplitude of the power spectrum of curvature perturbations 6 computed at k p ; n s and n T are, respectively, the scalar and the tensor spectral indices 7 . The value of k p is conventional and it corresponds to an effective harmonic ℓ eff ≃ 30. The figure for A R quoted in Eq. (1.3) corresponds to the value inferred from the WMAP 5-year data [3,4,5,6,7] in combination with the minimal ΛCDM model 8 . In the ΛCDM model the origin of A R stems from adiabatic curvature perturbations which are present after neutrino decoupling but before matter radiation equality (taking place at a redshift z eq = 3176 +151 −150 according to the WMAP 5-yr data [3,4,5]). The dominant component of curvature perturbations is adiabatic meaning that, over large scales, the fluctuations in the specific entropy are vanishing, at least in the minimal version of the model. The adiabatic nature of the fluctuations induces a simple relation between the first acoustic peak of the TT power spectra and the first anticorrelation peak of the TE power spectra [9]: this is, to date, the best evidence that curvature perturbations are, predominantly, adiabatic. It is useful to translate the comoving wave number k p into a comoving frequency ν p = k p 2π = 3.092 × 10 −18 Hz ≡ 3.092 aHz, (1.4) so, as anticipated, ν p is of the order of the aHz. The amplitude at the pivot scale 9 is controlled exactly by r T . The combined analysis of the CMB data, of the large-scale (LSS) structure data [13,14] and of the supernova (SN) data [15,16] can lead to quantitative upper limits on r T which are illustrated in Tabs. 1 and 2 as they are emerge from the combined analyses of different data sets.
The inferred values of the scalar spectral index (i.e. n s ), of the dark energy and dark matter fractions (i.e., respectively, Ω Λ and Ω M0 ), and of the typical wavenumber of equality k eq are reported in the remaining columns. While different analyses can be performed, it is clear, by looking at Tabs. 1 and 2, that the typical upper bounds on r T (k p ) range between, say, 0.2 and 0.4. More stringent limits can be obtained by adding supplementary assumptions. In Tab. 1 the quantity α S determines the frequency dependence of the scalar spectral index. In the simplest case α S = 0 and the spectral index Data r T n s α S Ω Λ Ω M0 k eq Mpc WMAP5 alone < 0. 58 Table 2: Same as in Tab. 1 but assuming no running in the (scalar) spectral index (i.e. α S = 0). is frequency-independent (i.e. n s does not run with the frequency). It can also happen, however, that α S = 0 which implies an effective frequency dependence of the spectral index. If the inflationary phase is driven by a single scalar degree of freedom (as contemplated in the minimal version of the ΛCDM scenario) and if the radiation dominance kicks in almost suddenly after inflation, the whole tensor contribution can be solely parametrized in terms of r T . The rationale for the latter statement is that r T not only determines the tensor amplitude but also, thanks to the algebra obeyed by the slow-roll parameters, the slope of the tensor power spectrum, customarily denoted by n T . To lowest order in the slow-roll expansion, therefore, the tensor spectral index is slightly red and it is related to r T (and to the slow-roll parameter) as where ǫ measures the rate of decrease of the Hubble parameter during the inflationary epoch 10 . Within the established set of conventions the scalar spectral index n s is given by n s = (1 − 6ǫ + 2η) and it depends not only upon ǫ but also upon the second slow-roll parameter η = M 2 P V ,ϕϕ /V (where V is the inflaton potential, V ,ϕϕ denotes the second derivative of the potential with respect to the inflaton field and M P = 1/ √ 8πG). It is sometimes assumed that also n T is not constant but it is rather a function of the wavenumber, i.e.
where α T now measures the running of the tensor spectral index.
As already mentioned, among the CMB experiments a central role is played by WMAP [3,4,5,6,7] (see also [8,9,10] for first year data release and [11,12] for the third year data release. In connection with [3,4,5,6,7], the WMAP 5-year data have been also combined with observations of the Acbar 11 experiment [17,18,19,20] . The TT, TE and, partially EE angular power spectra 12 have been measured by the WMAP experiment. Other (i.e. non space-borne) experiments are now measuring polarization observables, in particular there are • the Dasi (degree angular scale interferometer) experiment [21,22,23] operating at south pole ; • the Capmap (cosmic anisotropy polarization mapper) experiment [24,25] ; • the Cbi (cosmic background imager) experiment [26,27]; • the Quad experiment [28,29,30]; as well as various other experiments at different stages of development 13 . In the near future the Planck explorer satellite [31] might be able to set more direct limits on r T by measuring (hopefully) the BB angular power spectra.

The relic graviton spectrum in the ΛCDM model
Having defined the frequency range of the spectrum of relic gravitons, it is now appropriate to illustrate the possible signal which is expected within the ΛCDM scenario. In Fig. 2 the spectrum of the relic gravitons is reported in the case of the minimal ΛCDM scenario for different values of r T . On the 11 The Arcminute Cosmology Bolometer Array Receiver (ACBAR) operates in three frequencies, i.e. 150, 219 and 274 GHz. 12 Following the custom the TT correlations will simply denote the angular power spectra of the temperature autocorrelations. The TE and the EE power spectra denote, respectively, the cross power spectrum between temperature and polarization and the polarization autocorrelations. 13 Other planned experiments have, as specific target, the polarization of the CMB. In particular it is worth quoting here the recent projects Clover [32], Brain [33], Quiet [34] and Spider [35] just to mention a few.
horizontal axis the common logarithm of the comoving frequency is reported. The spectral energy density per logarithmic interval of frequency and in critical units is illustrated on the vertical axis. More quantitatively Ω GW (ν, τ 0 ) is defined as 14 Ω GW (ν, τ 0 ) = 1 ρ crit dρ GW d ln ν , (1.7) where ρ crit = 3H 2 0 M 2 P is the critical energy density. Since ρ crit depends upon H 2 0 (i.e. the present value of the Hubble rate), it is practical to plot directly h 2 0 Ω GW (ν, τ 0 ) at the present (conformal) time τ 0 . The proper definition of Ω GW (ν, τ 0 ) in terms of the energy-momentum pseudo-tensor in curved space-time is postponed to section 5. The salient features of the relic graviton spectra arising in the context of the ΛCDM scenario can be appreciated by looking carefully at Fig. 2.
The infra-red branch of the relic graviton spectrum (see also Fig. 1) extends, approximately, from ν p up to a new frequency scale which can be numerically determined by integrating the evolution equations of the tensor modes and of the background geometry across the matter-radiation transition. A semi-analytic estimate of this frequency is given by Hz. (1.8) At intermediate frequencies Fig. 2 exhibits a further suppression which is due to the coupling of the tensor modes with the anisotropic stress provided by the collisionless species which are present prior to matter-radiation equality. This aspect has been recently emphasized in Ref. [36] (see also [37,38,39]). Figure 2 assumes that the only colisionless species are provided by massless neutrinos, as the ΛCDM model stipulates and this corresponds, as indicated, to R ν = 0.405. The quantity R ν measures the contribution of N ν families of massless neutrinos to the radiation plasma: The frequency range of the suppression due to neutrino free-streaming extends from ν eq up to ν bbn which is given, approximately, by ν bbn = 2.252 × 10 −11 g ρ 10.75 Hz ≃ 0.01 nHz. (1.10) Both in Eqs. (1.8) and (1.10) Ω M0 and Ω R0 denote, respectively, the present critical fraction of matter and radiation with typical values drawn from the best fit to the WMAP 5-yr data alone and within the ΛCDM paradigm. In Eq. (1.10) g ρ denotes the effective number of relativistic degrees of freedom entering the total energy density of the plasma. While ν eq is still close to the aHz, ν bbn is rather in the nHz range. In Fig. 2 (plot at the left) the spectral index n T is frequency independent; in the plot at the right, always in Fig. 2, the spectral index does depend on the wavenumber. These two possibilities correspond, respectively, to α T = 0 and α T = 0 in Eq. (1.6). In the regime ν < ν eq a numerical calculation of the transfer function is mandatory for a correct evaluation of the spectral slope. In the approximation of a sudden transition between the radiation and matter-dominated regimes the spectral energy density goes, approximately, as ν −2+n T . The spectra illustrated have been computed within the approach developed in [40,41] and include also other two effects which can suppress the amplitude of the quasi-flat plateau. These effects are related to the contribution of the dark energy and to the evolution of the effective number of relativistic species. From a quantitative point of view both effects are, however, less relevant than neutrino free streaming.
Apart from the modification induced by the neutrino free-streaming the slope of the spectral energy density for ν > ν eq is quasi flat and it is determined by the wavelengths which reentered the Hubble radius during the radiation-dominated stage of expansion. The suggestion that relic gravitons can be produced in isotropic Friedmann-Robertson-Walker models is due to Ref. [42] (see also [43]) and was formulated before the inflationary paradigm. After the formulation of the inflationary scenario the focus has been to compute reliably the low frequency branch of the relic graviton spectrum. In [44,45,46] the low-frequency branch of the spectrum has been computed with slightly different analytic approaches but always assuming an exact de Sitter stage of expansion prior to the radiation-dominated phase. The analytical calculation (whose details will be described in section 6) shows that in the range ν p < ν < ν eq , the spectral energy density of the relic gravitons (see Eq. (1.7)) should approximately go as Ω GW (ν, τ 0 ) ≃ ν −2 . Within the same approximation, for ν > ν eq the spectral energy density is exactly flat (i.e. Ω GW (ν, τ 0 ) ≃ ν 0 ). This result, obtainable by means of analytic calculations (see also [47,48,49,50]), is a bit crude in the light of more recent developments. To assess the accurately spectral energy density it is necessary to take into account that the infrared branch is gradually passing from a quasi-flat slope (for ν > ν eq ) to the slope ν −2 which is the one computed within the sudden approximation [47,48,49,50]. It is useful to quote some of the previous reviews which covered, in a more dedicated perspective, the subject of the stochastic backgrounds of relic gravitons. The review article by Thorne [51] does not deal solely with relic graviton backgrounds while the reviews of Refs. [52,53,54] are more topical.
The flat plateau of the spectral energy density extends, approximately, between ν eq and a certain ν max . Also the maximal amplified frequency can be computed once the model of smooth transition between inflation and radiation is known. The smoothness of the transition determines specifically the precise amount of exponential suppression for ν > ν max . A simple estimate of ν max is given by GHz, (1.11) where, as in Eqs. (1.2) and (1.3), A R denotes the amplitude of the power spectrum of curvature perturbations evaluated at the pivot wavenumber ν p . It is worth noticing that between ν bbn and ν max there are approximately 20 orders of magnitude in frequency. In the ΛCDM scenario the spectrum has, in this range, always the same slope (i.e. n T is frequency-independent in Eq. (1.2)).
Some details of the calculations leading to the spectral energy densities illustrated in Fig. 2 can be found in sections 5 and 6. Without dwelling on the details it is however clear, as anticipated, that the constraints on the long wavelength gravitons make it difficult (if not impossible) to have a detectable spectral energy density at the scale of wide-band interferometers. The latter statement, valid in the minimal ΛCDM scenario, will be sharpened in the following subsection.

Short wavelength gravitons and wide-band interferometers
In the ΛCDM scenario the spectral energy density of the relic gravitons has its larger amplitude in the low-frequency branch. As the frequency increases the spectral energy density diminishes so that it is plausible to expect a rather small amplitude over the frequencies corresponding to wide-band interferometers (see, for instance, Fig. 2 for ν ≃ ν LV = 100 Hz).
Wide-band interferometers operate in a window ranging from few Hz up to 10 kHz (see also Fig. 1). The available interferometers are Ligo [55], Virgo [56], Tama [57] and Geo [58]. In loose terms these instruments are Michelson interferometers with two important differences: the mirrors are suspended and Fabry-Pérot cavities are used to increase the optical path of the photons. It would be too pretentious to describe in detail, in the present script, also the experimental apparatus and we therefore suggest Ref. [59] where the basics of wide-band interferometers are introduced in a self-contined perspective.
The sensitivity of a given pair of wide-band detectors to a stochastic background of relic gravitons depends upon the relative orientation of the instruments. The wideness of the band (important for the correlation among different instruments) is not as large as 10 kHz but typically narrower and, in an optimistic perspective, it could range up to 100 Hz. The putative frequency of wide-band detectors will therefore be indicated as ν LV , i.e. in loose terms, the Ligo/Virgo frequency. There are daring projects of wide-band detectors in space like the Lisa [60], the Bbo [61] and the Decigo [62] projects. The common feature of these three projects is that they are all space-borne missions and that they are all sensitive to frequencies smaller than the mHz (1 mHz = 10 −3 Hz). While wide-band interferometers are now operating and might even reach their advanced sensitivities during the incoming decade, the wished sensitivities of space-borne interferometers are still on the edge of the achievable technologies.
Since ν bbn < ν LV < ν max , wide-band interferometers represent an ideal instrument to investigate the relic graviton spectrum at high-frequencies. The spectral energy density of the relic gravitons produced within the ΛCDM model is quite minute and it is undetectable by interferometers even in their advanced version where the sensitivity is expected to improve by 5 or even 6 orders of magnitude in comparison with the present performances [63,64,65] (see also [66] and [67]). In Fig. 3 the spectral energy density is reported for ν = ν LV and always in the case of the prediction stemming from the minimal ΛCDM scenario. In Fig. 3, the common logarithm of the spectral energy density is illustrated  Figure 3: The spectral energy density of the relic gravitons in the context of the ΛCDM model evaluated at the Ligo/Virgo frequency as a function of the tensor-to-scalar ratio. To guide the intuition consider that (advanced) wide-band interferometers will achieve, optimistically, a sensitivity of 10 −10 in Ω GW (ν, τ 0 ). Taking into account the current bounds on r T (see Tabs. 1 and 2) it is clear that the ΛCDM scenario is definitively too small to provide a serious target for wide-band detectors.
as a function of the common logarithm of r T .
In Ref. [64] (see also [63,65]) the current limits on the presence of an isotropic background of relic gravitons have been assessed. According to the Ligo collaboration (see Eq. (19) of Ref. [64]) the spectral energy density of a putative (isotropic) background of relic gravitons can be parametrized as 15 : (1.12) The parametrization of Eq. (1.12) fits very well with Fig. 3 where the pivot frequency ν LV = 100Hz coincides with the pivot frequency appearing in the parametrization (1.12). For the scale-invariant case (i.e. β = −3 in eq. (1.12)) the Ligo collaboration sets a 90% upper limit of 1.20 × 10 −4 on the amplitude appearing in Eq. (1.12), i.e. Ω GW,−3 . Using different sets of data (see [63,65]) the Ligo collaboration manages to improve the bound even by a factor 2 getting down to 6.5 × 10 −5 . Thus Fig.  3 together with the upper limit of Eq. (1.12) shows that the current Ligo sensitivity is still too small to detect the relic graviton background arising within the ΛCDM paradigm.

Beyond the ΛCDM paradigm and high-frequency gravitons
In the case of an exactly scale invariant spectrum the correlation of the two (coaligned) LIGO detectors with central corner stations in Livingston (Lousiana) and in Hanford (Washington) might reach a sensitivity to a flat spectrum which is [68,69,70] where T denotes the observation time and SNR is the signal to noise ratio. Equation (1.13) is in close agreement with the sensitivity of the advanced Ligo apparatus [55] to an exactly scale-invariant spectral energy density [71,72,73,74]. Equation (1.13) together with the plots of Fig. 3 suggest that the relic graviton background predicted by the ΛCDM paradigm is not directly observable by wide-band interferometers in their advanced version.
CMB observations probe the aHz region of the spectral energy density of Fig. 2. Wide-band interferometers probe a frequency range between few Hz and 10 kHz. In both ranges, the signal of the ΛCDM scenario might be too small to be directly detectable. In Fig. 4 the spectral energy density is computed in an extension of the ΛCDM paradigm which has been dubbed tensor-ΛCDM (TΛCDM for short) [40,41]. In the TΛCDM scenario the transition from the inflationary phase to the radiation-dominated epoch is mediated by a rather long stiff phase. By stiff phase we mean a phase where the total sound speed of the plasma is larger than the sound speed of a radiation-dominated plasma (i.e. 1/ √ 3 in natural units). In the simplest realization of the scenario the barotropic index w t is constant during the stiff phase. For instance, in Fig. 4 the cases w t = 1 and w t = 0.6 have been illustrated. Since w t denotes the ratio between the total pressure and the total energy density, it is rather plausible to demand that w t ≤ 1. The latter requirement implies that the sound speed is always smaller than the speed of light. As suggested in [75] (see also [76,77]) the presence of a stiff phase can have the effect of increasing the spectral energy density at high frequencies. The increase In both plots the other cosmological parameters have been chosen to coincide with the best fit to the WMAP 5-yr data alone.
takes place for frequencies larger than the mHz and is typically maximal in the GHz region. The spectral energy densities illustrated in Fig. 4 suggest that it is not impossible to imagine situations where the spectral energy density of the relic gravitons satisfies all the constraints demanded by CMB physics but, at the same time, it is sufficiently large to be observed by wide-band interferometers. The results reported in Fig. 4 refer to the minimal TΛCDM model where only one post-inflationary phase stiffer than radiation is contemplated. The barotropic index could have, however, a more complicated dependence. Already in the examples of Fig. 4 the numerical integration implies that the barotropic index does depend, effectively, upon the scale factor (see, e. g., discussions in section 6 on the transfer function for the spectral energy density).
The comparison of Fig. 2 and 4 suggests, in short, the following subjects of reflection: • the theoretical error in the estimate of the spectral energy density increases with the frequency; • departures from the standard post-inflationary thermal history can be directly imprinted in the primordial spectrum of the relic gravitons; • in the incoming decade the observations of wide-band interferometers could be analyzed in conjunction with more standard data sets (i.e. CMB data supplemented by large-scale structure data and by the observations of type Ia supernovae) to constrain the spectral energy density of the relic gravitons both at small and at high frequencies.
The presence of post-inflationary phases stiffer than radiation is, after all, rather natural and this was the original spirit of [75]. We do not know which was the rate of the post-inflationary expansion and since guesses cannot substitute experiments it would be productive to use the TΛCDM paradigm as reference model for a unified analysis of the low-frequency data stemming from CMB and of the high-frequency data provided by wide-band interferometers. Already in [75] (see also [76,77]) a rather special candidate for a post-inflationary phase stiffer than radiation was the case when the sound speed equals the speed of light, i.e. the case when the energy density of the sources driving the geometry is dominated by the kinetic term of a (minimally coupled) scalar field. This particular case was also prompted by various classes of quintessence models. A specific example of this dynamics was provided in [78].
A more detailed account of the techniques leading to Fig. 4 will be swiftly presented in section 6 and can be found in [40,41]. Without going through the details it is however important to stress that the calculations should be accurate enough not only in the high-frequency region but also in the low-frequency part of the spectrum. Indeed, as stressed above, one of the purposes of the TΛCDM scenario is to convey the idea that low-frequency and high-frequency measurements of the relic graviton background can be analyzed in a single theoretical framework.

The millisecond pulsar bound and the nucleosynthesis constraint
The spectral energy density of the relic gravitons must be compatible not only with the CMB constraints (bounding, from above, the value of r T ) but also with the pulsar timing bound [79,80] and the big-bang nucleosynthesis constraints [81,82,83]. The pulsar timing bound demands Ω(ν pulsar , τ 0 ) < 1.9 × 10 −8 , ν pulsar ≃ 10 nHz, (1.14) where ν pulsar roughly corresponds to the inverse of the observation time during which the pulsars timing has been monitored. The spectral energy densities illustrated in Figs. 2 and 4 satisfy the pulsar timing bound.
The most constraining bound for the high-frequency branch of the relic graviton spectrum is represented by big-bang nucleosynthesis. Gravitons, being relativistic, can potentially increase the expansion rate at the BBN epoch. The increase in the expansion rate will affect, in particular, the synthesis of 4 He. To avoid the overproduction of 4 He the expansion rate the number of relativistic species must be bounded from above.
The BBN bound is customarily expressed in terms of (equivalent) extra fermionic species. During the radiation-dominated era, the energy density of the plasma can be written as ρ t = g ρ (π 2 /30)T 4 where T denotes here the common (thermodynamic) temperature of the various species. An (ultra)relativistic fermion species with two internal degrees of freedom and in thermal equilibrium contributes 2 · 7/8 = 7/4 = 1.75 to g ρ . Before neutrino decoupling the contributing relativistic particles are photons, electrons, positrons, and N ν = 3 species of neutrinos, giving The neutrinos have decoupled before electron-positron annihilation so that they do not contribute to the entropy released in the annihilation. While they are relativistic, the neutrinos still retain an equilibrium energy distribution, but after the annihilation their (kinetic) temperature is lower, T ν = (4/11) 1/3 T . Thus after electron-positron annihilation. By now assuming that there are some additional relativistic degrees of freedom, which also have decoupled by the time of electron-positron annihilation, or just some additional component ρ X to the energy density with a radiation-like equation of state (i.e. p X = ρ X /3), the effect on the expansion rate will be the same as that of having some (perhaps a fractional number of) additional neutrino species. Thus its contribution can be represented by replacing N ν with N ν + ∆N ν in the above. Before electron-positron annihilation we have ρ X = (7/8)∆N ν ρ γ and after electron-positron annihilation we have ρ X = (7/8)(4/11) 4/3 ∆N ν ρ γ ≃ 0.227 ∆N ν ρ γ . The critical fraction of CMB photons can be directly computed from the value of the CMB temperature and it is notoriously given by h 2 0 Ω γ ≡ ρ γ /ρ crit = 2.47 × 10 −5 . If the extra energy density component has stayed radiation-like until today, its ratio to the critical density, Ω X , is given by (1.17) If the additional species are relic gravitons, then [81,82,83]:  [82] and therefore the bounds on ∆N ν can be translated into bounds on the energy density of the relic gravitons.
A review of the constraints on ∆N ν can be found in [82]. Depending on the combined data sets (i.e. various light elements abundances and different combinations of CMB observations), the standard BBN scenario implies that the bounds on ∆N ν range from ∆N ν ≤ 0.2 to ∆N ν ≤ 1. Similar figures, depending on the priors of the analysis, have been obtained in a more recent analysis [83]. All the relativistic species present inside the Hubble radius at the BBN contribute to the potential increase in the expansion rate and this explains why the integral in Eq. (1.18) must be performed from ν bbn to ν max (see also [76] where this point was stressed in the framework of a specific model). The existence of the exponential suppression for ν > ν max (see Figs. 4) guarantees the convergence of the integral also in the case when the integration is performed up to ν → ∞. The constraint of Eq. (1.18) can be relaxed in some non-standard nucleosynthesis scenarios [82], but, in what follows, the validity of Eq. (1.18) will be enforced by adopting ∆N ν ≃ 1 which implies, effectively (1.19) The spectral energy densities illustrated in Figs. 2 and 4 are both compatible with the big-bang nucleosynthesis bound. Thus the big-bang nucleosyntheis constraint does not forbid a potentially detectable signal in the high-frequency branch of the relic graviton spectrum. Potential deviations of the thermal history of the plasma must anyway occur before big-bang nucleosynthesis.
2 The polarization of relic gravitons and of relic photons

Basic notations
As discussed in the introduction, in the ΛCDM paradigm the background line element can be written where, in the spatially flat case, γ ij will coincide with δ ij and the Friedmann-Lemaître equations can be written as where H = a ′ /a; the prime denotes a derivation with respect to the conformal time coordinate τ . The Hubble rate is customarily defined in the synchronous frame where the time coordinate (conventionally denoted by t) obeys dt = a(τ )dτ . Denoting with a dot a derivation with respect to the cosmic time t, H =ȧ/a, and, by definition, H = aH. In Eqs. (2.2)-(2.4) ρ t and p t are, respectively, the total energy density and the total pressure of the plasma, i.e.
The total matter fraction of the critical energy density, i.e. Ω M0 = ρ M0 /ρ crit consists of baryons and (i.e. ρ b ) and cold dark matter particles (i.e. ρ c ). In Tabs. 1 and 2 the values of Ω M0 are given as they are inferred within the ΛCDM scenario. In similar terms Ω Λ = ρ Λ /ρ crit denotes the critical fraction of dark energy. In what follows, if not otherwise stated, the cosmological parameters will be fixed to the best fit of the WMAP-5yr data alone, i.e.
At the beginning of the previous section we started by stressing analogies and differences between relic gravitons and relic photons. The most important one is that both gravitons and photons carry two polarizations. This observation is important for a quantitative understanding of the present endevours aimed at measuring the E-mode and the B-mode polarization of the CMB. In the present section the description of the polarization of the gravitons will be developed by stressing, when possible, the analogy with polarization observables of the electromagnetic field.

Linear and circular tensor polarizations
Recalling that i, j, k, ... are indices defined on the three-dimensional Euclidean sub-manifold, the tensor fluctuations of the geometry are parametrized in terms of the rank-two tensor h ij where ∇ i is the covariant derivative with respect to γ ij ; if γ ij = δ ij , ∇ i = ∂ i . In Eq. (2.8) the subscript refers to the tensor nature of the fluctuation while the superscript denotes the perturbative order. The tensor fluctuation h ij ( x, τ ) can be decomposed in terms of the two linear polarizations, i.e.
where λ = ⊕, ⊗ denote the two polarizations and where Ifk coincides with the radial direction, the unit vectorsk,â andb can be chosen, in spherical coordinates, as:k = (sin ϑ cos ϕ, sin ϑ sin ϕ, cos ϑ), Since ǫ As in the case of electromagnetic waves, it is often desirable to pass from the linear to the circular polarizations: The transformation properties of the circular polarization under a rotation in the plane orthogonal to the direction of propagation are closely analog to the transformation properties, under the same rotation, of the polarization of the electromagnetic field. This analogy will now be exploited to introduce the E-mode and B-mode polarization.
Before proceeding with the discussion it is appropriate to recall a very basic aspect of rotations which can have, however, some confusing impact of the polarization analysis especially in the case of the tensor modes. Consider, for simplicity, a coordinate system characterized by two basis vectors, i.e.ê x =r cos ϑ andê y =r sin ϑ. If we now perform a clockwise (i.e. right-handed) rotation of the axesê x andê y , the rotated basis (ê ′ x ,ê ′ y ) will be given as in Eqs. (2.22) and (2.23) by replacing x ,ê ′ y ) and (â,b) → (ê x ,ê y ). Some authors, for different reasons, instead of rotating the coordinate system prefer to rotate the polarization vector. If angles are in the right-handed sense for the rotation of the axes, they are in the left-handed sense for the rotation of the vectors.

Polarization of the CMB radiation field
The radiation field can be described by the polarization tensor, i.e.
where E i and E j are the electric components of the radiation field. Assuming, for sake of simplicity, that the radiation field propagate along theẑ axis, then the various entries of ρ ij can be written in a matrix form where x and y denote the components of the electric field orthogonal to the direction of propagation coinciding, in this set-up, with the third Cartesian axis. A full description of the radiation field can be achieved by studying the four Stokes parameters [84] conventionally named I, Q, U and V : It is immediate, from the definitions (2.29) and (2.30), to write the intensities of the radiation field along the different Cartesian axis as a function of the Stokes parameters, i.e.
Equations (2.31) and (2.32) can be inserted back into Eq. (2.28) with the result that where From Eq. (2.33) it also follows that where 1 is the 2 × 2 unit matrix and

E-and B-modes
The fluctuations of the geometry induce fluctuations of the Stokes parameters whose spectral properties are, ultimately, the aim of CMB polarization experiments. In general terms the fluctuation of each Stokes parameter can be written as In Eqs. (2.38)-(2.40) the superscript reminds that the various fluctuations of the Stokes parameters are induced, respectively, by the tensor, scalar and vector modes of the geometry. While some of the results of the present section will be generally valid, the focus, in what follows, will be on the tensor contribution. Defining the two linear combinations ∆ ± (n, τ ) = ∆ Q (n, τ ) ± i∆ U (n, τ ), (2.41) and denoting with a tilde the transformed quantities, Eq. (2.37) implies that ∆ ± (n, τ ) transform as In more general terms, consider a generic function ofn (be it f (n)). Under a rotation of an angle α on the plane orthogonal ton, f (n) is said to transform as a function of (integral) spin-weight ±s providedf (n) = e ∓isα f (n). (2.43) In other words ∆ + (n) and ∆ − (n) transform, respectively, as functions of spin weight +2 and −2.
The circular polarizations of the gravitons introduced in Eqs. (2.24)-(2.25) transform, respectively, as functions of spin weight +2 and −2. The brightness perturbations for the intensity of the radiation field (i.e. ∆ I (n)) transform, on the contrary, as quantities of spin weight 0. The fluctuations in the intensity of the radiation field, being a spin-0 quantity, can be expanded in ordinary spherical harmonics as The spin-s quantity will naturally be expanded in terms of a generalization of the ordinary spherical harmonics which are called spin-s spherical harmonics or also spin-weighted spherical harmonics. Owing to this observation, ∆ ± (n, τ ) can be expanded in terms of spin-±2 spherical harmonics ±2 Y ℓ m (n), i.e.
Given a quantity of spin-weight s it is possible to construct quantities of spin-weight 0 by the repeated use of appropriate differential operators which can either raise or lower the spin-weight of a given function (see subsection 2.5 for a specific discussion). Consequently, from ∆ + (n, τ ) and ∆ − (n, τ ) it is possible to construct two fluctuations of spin 0 which can be eventually expanded in ordinary spherical harmonics. By demanding that the two fluctuations of spin-weight 0 are eigenstates of parity the E-and B-modes are defined as ℓ m the angular power spectra can be defined. In particular the EE, BB, TT and TE angular power spectra are given by: where ... denotes the ensemble average. Two further power spectra can be defined and they are: (2.50) Overall, the existence of linear polarization allows for 6 different power spectra.
In the minimal version of the ΛCDM paradigm the adiabatic fluctuations of the scalar curvature lead to a polarization which is characterized exactly by the condition a 2, ℓm = a −2, ℓm , i.e. a (B) ℓm = 0. This observation implies that, in the ΛCDM scenario, the non-vanishing angular power spectra are given by the TT, EE and TE correlations. In the TΛCDM scenario the TT, EE and TE angular power spectra are supplemented by a specific prediction for the B-mode autocorrelation (see section 7).

Spin-2 spherical harmonics
Spherical harmonics of higher spin appear in matrix elements calculations in nuclear physics (see e.g. the classic treatise of Blatt and Weisskopf [85], and, in a similar perspective the book of Edmonds [86]). The comprehensive treatments of Biedenharn and Louk [87] and of Varshalovich et al. [88] can also be usefully consulted.
The spin-s harmonics have been introduced, in their present form, by Newman and Penrose [89] and their group theoretical interpretation has been discussed in [90]. The spin-s spherical harmonics have been applied to the discussion of CMB polarization induced by relic gravitons in a number of papers [91,92,93]. They are rather crucial in the formulation of the so-called total angular momentum approach. Discussions of the spin-weighted spherical harmonics in a cosmological context can also be found in [94,95]. The spin weighted spherical harmonics will now be introduced by following the spirit of Ref. [90] which has been also used, with different conventions, in [91]. In subsection 2.6 the (equivalent) approach of [92,93] will be more specifically outlined.
Th functions ±2 Y ℓ m (n) appearing in Eq. (2.45) are the spin-2 spherical harmonics [90]. Consider the representations of the operator specifying three-dimensional rotations, i.e.R; this problem is usually approached within the matrix element, i.e. D where α, β and γ (set to zero in the above definition) are the Euler angles defined as in [96]. If s = 0, D 0, m (α, β, 0) = (2ℓ + 1)/4πY ℓ m (α, β) where Y ℓ m (α, β) are the ordinary spherical harmonics. The spin-s spherical harmonics can be obtained from the spin-0 spherical harmonics by using repeatedly the differential operators: (2.53) The notation spelled out in Eqs. (2.52) and (2.53) (which is not usual) will be employed to emphasize the interpretation of K (s) ± as ladder operators (see [90]). The operator K (s) + raises the spin weight of a function by one unit. Consider, therefore, the ordinary spherical harmonics Y ℓ m (n) defined as where P ℓ (µ) are the Legendre polynomials and P m ℓ (µ) the associated Legendre functions. It is appropriate to mention here that the factor (−1) m (i.e. Condon-Shortley phase) can either be included in the normalization factor N m ℓ or (as it has been done) in the definition of the associated Legendre functions appearing in Eq. (2.55). When using the recurrence relations of the associated Legendre functions the Condon-Shortley phase introduces a sign difference every time m is odd. The conventions expressed by Eqs. (2.54) and (2.55) will be followed throughout the present discussion and, in particular, in section 7 where the correlation functions of the E-modes and of the B-modes will be specifically computed with different techniques.
According to Eq. (2.43), Y ℓ m transform with s = 0, i.e. they have spin weight 0. By applying once K (0) + (n) to Y ℓ m (n) we do get a function of spin weight 1, i.e.
We can then apply once more 16 K (1) The result of this simple manipulation will be The right hand side of Eq. (2.57) is +2 Y ℓ m (ϑ, ϕ) (up to an overall normalization). Including the appropriate normalization factor, ±2 Y ℓ m (ϑ, ϕ), i.e. the spherical harmonics of spin weight s = ±2 are given by: The spin weights s = ±2 are both needed since the transformation of the polarization involve both spin weights (see Eq. (2.45)). In fact, since ±2 Y ℓ m (ϑ, ϕ) form a complete and orthogonal basis on the sphere, i.e.
The coefficients off the expansion will be given by Integrating by parts in Eqs. (2.61) and (2.62) allows for a different form of the expansion coefficients a ±2, ℓ m : where, as already mentioned, N ℓ = (ℓ − 2)!/(ℓ + 2)!. In Eqs. (2.63) and (2.64) there appear only ordinary (i.e. spin-weight 0) spherical harmonics. This occurrence suggests a complementary approach to the problem: instead of expanding ∆ ± (n, τ ) in terms of spin-2 spherical harmonics, fluctuations of spin-weight 0 can be directly constructed (in real space) from ∆ ± (n, τ ) by repeated application of the ladder operators defined in Eqs. (2.52) and (2.53).
The E-mode and B-mode polarization in real space are then, in explicit terms: The quantities ∆ E (n, τ ) and ∆ B (n, τ ) can be expanded in terms of ordinary spherical harmonics, as already suggested in Eq. (2.46): The "electric" and "magnetic" components of polarization are eigenstates of parity and may be defined, from a ±, ℓ m as already mentioned in Eq. (2.47): Under parity the components appearing in Eqs. A parity transformation (i.e. a space inversion) implies, in spherical coordinates, that The transformation (2.71) implies that the two basis vectors defined in Eq. (2.70) transform asê 1 →ê 1 andê 2 → −ê 2 , i.e. whileê 1 does not changeê 2 flips its sign under space inversion. It follows that space-inversion does not flip the sign of ∆ E (n) but it does flip the sign of ∆ B (n), i.e. under the The contribution of long wavelength gravitons to Eqs. (2.72) and (2.73) will be discussed in section 7. It is often useful to observe that the differential operators appearing in the definition of the spinweighted spherical harmonics (see, e.g. Eq. (2.58)) can be expressed in terms of the usual differential operators arising in the theory of the orbital angular momentum in non-relativistic quantum mechanics (see, e. g. [96]). Indeed, recalling that it can be easily deduced that Equations (2.76)-(2.78) allow often to express combinations of spin-2 spherical harmonics in terms of ordinary (i.e. spin-weight 0) spherical harmonics using the properties of the ladder operators associated to the (orbital) angular momentum, i.e. L ± : where L ± and L z obey the well known commutation relations [L ± , L z ] = ∓L ± and [L + , L − ] = 2L z .
Looking at Eq. (2.79) it is tempting draw a parallel between the (orbital) ladder operators and the ladder operators raising (or lowering) the spin weight of a given function (see Eqs. (2.52) and (2.53). This problem has been discussed and solved in [90]. It is possible to formulate the parallel in terms of a putative O(4) group. Half of the generators will be connected with the orbital angular momentum operators, while the other half will allow to increase (or decrease) the spin weight of a given function. The two sets of generators commute. The operators K (s) ± are not directly, though, the ladder operators stemming from the second set of generators. This has to do with the fact that in Eq. (2.51) the third Euler angle (i.e. γ) has been fixed to zero. The K

Polarization on the 2-sphere
In a more geometric perspective, the spin-2 harmonics are introduced by describing the polarization tensor on the 2-sphere which represents the microwave sky. In Eq. (2.34) the tensor P ij describes the properties of the radiation field and it is symmetric and trace-free (i.e. P ij = P ji and P i i = 0). Equation (2.34) holds in flat space-time. On the 2-sphere the line element can be written as The polarization matrix P ij will now be generalized as satisfying P ab = P ba , and g ab P ab = 0, wheren is a unit vector in the direction (ϑ, ϕ). The sign of the off-diagonal entries in Eq. (2.81) is opposite with respect to the one obtained in Eq. (2.34). This is just because we want to match with the conventions adopted, for instance, in [93,94,95]. To avoid possible confusions, furthermore, the Latin indices a, b, c, d, .... run over the two-dimensional space.
As already mentioned, for scalar functions defined on the 2-sphere, such as the temperature anisotropies, the spherical harmonic functions Y ℓm (n) are the complete orthonormal basis. For the 2 × 2 tensors defined on the 2-sphere, such as P ab in Eq.(2.81), the complete orthonormal set of tensor spherical harmonics can be written as [93,94,95]: where ∇, in this subsection, denotes a covariant derivation on the 2-sphere, e.g.
In Eq. (2.83) the normalization factor is given by 18 is the Levi-Civita symbol on the 2-sphere. The differential operators acting in Eqs. (2.82) and (2.83) are interpreted as a generalized gradient and curl operators, i.e.
The explicit form of the various components of Y G (ℓm)ab and Y C ab can be computed. For instance using Eqs. (2.82), (2.84) and (2.86), the explicit components of Y G (ℓm)ab : The same exercise can be conducted for the various components of Y (ℓm)ab and Y (C) (ℓm)ab can be written in the form of 2 × 2 matrices: and as In terms of the spin-2 harmonics ±2 Y (lm) (n) which is Eq. (2.58). Using the orthonormality of the spherical harmonics Y ℓm (n) it is easy to prove the orthonormality conditions, i.e.
dn Y where dn = sin ϑdϑdϕ denotes, as usual, the integration over the solid angle. Since Y G (ℓm)ab and Y C (ℓm)ab form an appropriate orthonormal basis, the polarization can be expanded as where expansion coefficients a (G) ℓm and a (C) ℓm represent the electric and magnetic type components of the polarization, respectively. Note that the sum starts from ℓ = 2, since relic gravitons generate only perturbations of multipoles from the quadrupoles up. The expansion coefficients are obviously ℓm already introduced in Eq. (2.68). The relation between the two sets of expansion coefficients is simply: The two approaches to the spin weighted spherical harmonics described in the present section are equivalent and can be used interchangeably depending upon the specific problem.
3 The action of the relic gravitons

Second-order fluctuations of the Einstein-Hilbert action
By perturbing the Einstein-Hilbert action, to second-order in the amplitude of the tensor fluctuations we have, formally, that: where R denotes the background Ricci scalar; δ (1) R and δ (2) R denote respectively, the first and secondorder fluctuations of R = g µν R µν . In Eq. (3.1) the possible coupling to the anisotropic stress has been neglected. This is customary during the early evolution of the geometry since, in the context of the ΛCDM paradigm, during the early inflationary phase the sources of anisotropic stress can be safely ignored unless the number of effective e-folds is close to minimal. Later on the anisotropic stress of the fluid plays a role and cannot be neglected at least if we aim at a reasonable quantitative discussion of the relic graviton spectrum (see also section 1 and Fig. 3). By introducing the first-order fluctuations of the background geometry g µν we have that (3.5) The first-and second-order fluctuations of the Christoffel connections are: where the prime denotes a derivation with respect to the conformal time coordinate. Using the result of Eq. (3.6) the first-and second-order fluctuations of the Ricci tensor can be written in explicit terms: The Ricci scalar is zero to first order in the tensor fluctuations, i.e. δ t R = 0. This is due to the traceless nature of these fluctuations. To second-order, however, δ (2) t R = 0 and its form is: Using the results of Eqs. (3.7)-(3.10) into Eq. (3.1) the second-order action for the tensor modes reads, up total derivatives, where, as already mentioned in section 1,

Lagrangian densities
The action (3.11) can be written in various ways which differ by the addition (or subtraction) of a total conformal time derivative. Recalling the standard notations the Langrangian density L (1) ( x, τ ) can be recast in the form where the canonical amplitude h has been introduced By now introducing the canonical amplitude as ah = µ, Eq. (3.13) can be transformed as If a total τ derivative is dropped, an equivalent form of the Lagrangian density can be obtained All the three Lagrangian densities of Eqs. (3.14), (3.16) and (3.17) lead to the same Euler-Lagrange equations.

Hamiltonian densities
In Eq. (3.14) the canonical field is h and the canonical momentum is Π = ah ′ . Conversely, in Eq. (3.16) the canonical field is µ and the associated canonical momentum isπ = µ ′ − Hµ. Finally, according to Eq. (3.17) the canonical momentum is π = µ ′ while the canonical field is always µ.
By taking the functional derivative of G (1)→ (2) [µ, Π] with respect to Π and with respect to µ we get, up to a sign, the connection between the new and old pivot variables, namely: Since the generating functional depends explicitly upon time, the new Hamiltonian will be related to the old one through a partial time derivative of the generating functional, i.e.
depending upon the old coordinates (i.e. µ) and upon the new momenta (i.e. π). The relations between the new and old variables are given bỹ stipulating that, in this case, the canonical momentum gets shifted by Hµ while the canonical field is left invariant. Since the generating functional depends explicitly upon the conformal time coordinate, we will simply have that

Evolution equations in different regimes
From Eq. (3.11) the evolution equations of h j i will be given by The canonical field h (see Eq. (3.15)) will also obey Eq. (3.27). The Hamilton equations derived from Eq. (3.18) read: which has exactly the same content as Eq. (3.27). In similar terms the Hamilton's equations can be derived from Eq. (3.19) and the result is Bearing in mind that ah = µ, Eqs. . It is practical, for some applications, to change the time parametrization. For instance, in terms of the rescaled time coordinate σ we will have that the evolution for the canonical amplitude h obeys the simple equation Before concluding this section it should be pointed out that Eq. (3.27) is accurate as long as the sources of anisotropic stress are totally absent. This approximation is, strictly speaking, unrealistic. Indeed we do know that there are sources of anisotropic stress. Typically, after neutrino decoupling, the neutrinos free stream and the effective energy-momentum tensor acquires, to first-order in the amplitude of the plasma fluctuations, an anisotropic stress, i.e.
where Π j i is the contribution of the anisotropic stress, satisfying ∇ i Π i j = 0 and Π i i = 0. The presence of the anisotropic stress clearly affects the evolution of the tensor modes. To obtain the wanted equation we perturb the Einstein equations to first-order and we get: This form of the evolution equation for the tensor modes is the one required to compute the effects related to the finite value of the anisotropic stress.

Quantization of the tensor modes
There are analogies between the quantum state of relic gravitons and the quantum treatment of visible light. Quantum effects are not crucial to treat first-order interference of the radiation field (i.e. Young interferometry) [97]. First-order interference in quantum optics correspond to the calculation of the two-point function of the relic gravitons. Quantum effects arise, in optics, from second-order interference, i.e. when computing (and measuring) the interference between the intensities of the radiation field. Second-order interference effects are associated with the possibilities of counting photons and have been pioneered by Hanbury-Brown and Twiss in the early fifties [97,98]. Hanbury-Brown-Twiss interferometry is based on photon counting statistics.
Having said that we are not even close (experimentally) to study graviton counting statistics (as we do it with the photons), second order interference effects would allow, in principle, to assess the coherence properties of relic graviton backgrounds. The quantum state of the relic gravitons can be described in terms of a generalized coherent state usually called squeezed state. Squeezed states can be described in terms of quadrature operators where one of the modes of the radiation field is always broadened by the time evolution, while the other one is squeezed.

Heisenberg description
The quantization of the canonical Hamiltonian of Eq. (3.20) is performed by promoting the normal modes of the action to field operators in the Heisenberg description and by imposing (canonical) equal-time commutation relations: The operator corresponding to the Hamiltonian (3.20) becomes: In Fourier space the quantum fieldsμ andπ can be expanded aŝ Demanding the validity of the canonical commutation relations of Eq. (4.1), the Fourier components must obey: (4.5) 19 In order to derive the following equation, the relationsμ † − k The evolution ofμ andπ is therefore dictated, in the Heisenberg representation, by: where, as usual, unitsh = 1 are assumed. Using now the mode expansion (4.3) and the Hamiltonian in the form (4.5) the evolution for the Fourier components of the operators iŝ It is not a surprise that the evolution equations of the field operators, in the Heisenberg description, reproduces, forμ k the classical evolution equation derived before. The general solution of the system is thenμ where the mode functions obey: If the form of the Hamiltonian is different by a time-dependent canonical transformation, also the canonical momenta will differ and, consequently, the relation of g to f may be different. For instance, in the case of the Hamiltonian of Eq. (3.19) we will have, instead, Since, by construction, the Hamiltonians of Eqs. (3.19) and (3.20) are related by canonical transformations, the mode functions of Eqs. (4.11) and (4.12) will have both to obey Eq. (4.13). In different terms, the commutation relations between field operators should be preserved by the time evolution and this is equivalent to the Wronskian normalization condition of Eq. (4.13).

Generalized coherent states of relic gravitons
Consider the Hamiltonian given in Eq. (3.19) in the spatially flat case: dropping, for simplicity, the tilde from the momenta. Defining the creation and annihilation operatorŝ and recalling thatμ † − k = µ k andπ † − k = π k , Eq. (4.15) implŷ , inserting Eq. (4.16) into Eq. (4.14),Ĥ gw can be written asĤ The evolution ofâ k (τ ) andâ † − k (τ ) obeys: The solution of Eq. (4.18) is:â where τ 0 is the initial integration time. The unitarity of the time evolution demands that |u k (τ )| 2 − |v k (τ )| 2 = 1. A useful parametrization of u k (τ ) and v k (τ ) is given in terms of a real amplitude and two phases as: Equation (4.18) determine the evolution equations for u k (τ ) and v k (τ ). Using then Eq. (4.21) the evolution equations for r k , ϕ k and ϑ k can be obtained: There is a slight difference in the normalizations adopted between Eqs. (4.9)-(4.10) and Eqs. (4.25)- (4.26). This difference is due to the fact that, in Eqs. (4.9)-(4.10) the mode functions f k are normalized, asymptotically, in such a way that f k → 1/ √ 2k. In Eqs. (4.15)-(4.16) the factors √ 2k and k/2 have been included in the definition creation and annihilation operators.
After simple calculations the two-point functions of the field operators and of their related canonical momenta becomes : (4.28) where r = | x − y|. Again, as already remarked, the non-strandard pre-factors apperaing in the Fourier amplitudes of Eqs. (4.27) and (4.28) are a consequence of the normalizations of Eq. (4.15). In the limit | x − y| → 0 (and making use of the definitions of Eq. (4.21)), Eqs. (4.27) and (4.28) lead to Equations (4.29) and (4.30) show that the canonical field is broadened while the conjugate momentum gets squeezed by keeping constant the product of their respective root mean squares. The latter behaviour is evident as soon as the relevant wavelengths are larger than the Hubble radius. To demonstrate the two previous statements consider, indeed, Eqs. (4.22) and (4.23). Their solution for wavelengths larger than the Hubble radius (i.e. to leading order in kτ < 1) is: The rationale for the behaviour exhibited by Eq. (4.31) can be understood also from a slightly different perspective. The relation between u k (τ ) and v k (τ ) and the tensor mode functions f k (τ ) and g k (τ ) is where the evolution of f k (τ ) and g k (τ ) is obtained by solving Eq. (4.12). Indeed, by solving Eq. (4.12) in the limit From Eq. (4.35) the first derivative of f k with respect to τ is nothing but By computing g k (τ ) from Eq. (4.35) it is clear that, in the limit kτ ≪ 1 which has exactly the same physical content of Eqs. (4.32) and (4.33). When the Universe expands, g k (τ ) decreases and that the solution associated with A 2 (k) becomes progressively subleading. However, this observation does not imply that g k (τ ) disappears since the evolution must be unitary. This feature of squeezed quantum state suggests the possibility of associating an effective entropy to the process of graviton production [104,105,106].
In the Schrödinger description quantum state of the relic gravitons is closely related to the squeezed states of the radiation field [99] (see also [100,101]). Defining, for practical reasons, the wavefunction of the ground state will be, for a given mode k, where, as previously discussed, ϕ k → 0 has been assumed. In Eq. (4.39) Σ k (r k , ϕ k ) is the so-called two-mode squeezing operator [102,103] which will now be written in its generic form: where the creation and destruction operators are the ones computed in τ 0 , i.e.â k ≡â k (τ 0 ) by definition of Schrödinger description. The state |0 is annihilated both byâ k and byâ − k . These two-modes appear simultaneously since gravitons are produced from the vacuum whose total momentum vanish. The x andx have been dubbed, in the literature, as superfluctuant operators (see, e. g., [104,105,106]).
The statistical properties of squeezed states can be addressed by employing a useful analogy with quantum optics. The idea is to pretend to resolve single gravitons and to study the statistical properties of the second-order interference effects. In the case of the gravitons this problem is addressed by defining the (normalized) Glauber correlation function not for the photons (as customarily done) but for the gravitons. For simplicity let us consider a single mode of the field. In this case the normalized Glauber intensity correlation can be written as the colons denote normal ordering andÎ denotes the operator corresponding to the intensity of the radiation field. The normal ordering is related to the fact that, in the optical domain, most measurements of the electromagnetic field are based on the absorption of photons via the photoelectric effect 20 . In the case of a single mode of the field Eq. (4.41) can also be written as In Eq. (4.42) the normalized two point function is written for coincident spatial points (i.e. ∆t = 0). The Hanbury-Brown and Twiss experiment, in some sense, probes directly the properties of g (2) (0). In the case of coherent states, g 2 (0) = 1: this is the case when the photoelectric counts obey a Poisson statistics. In the case of chaotic light, the joint detection probability greater than that for two independent events. This can be verified from Eq. (4.42) by assuming that the state of the photons/gravitons is given by a thermal mixture. Equation (4.42) can also be recast in the form where (∆n) 2 = n 2 − n 2 is the photon number variance andn =â †â . From Eq. (4.43) it is also customary to define the Mandel's Q-parameter, i.e.
which vanishes exactly in the case of a coherent state. The latter statement can be easily appreciated since, for a coherent state,â|α = α|α . Thus, (∆n) 2 = |α| 2 = n . In the case of chaotic (thermal) light it turns out that g (2) (0) = 2. This result can be easily drived by using the so-called Glauber-Sudarshan P -representation of the density matrix, i.e.
where n is the Bose-Einstein occupation number. In the P -representation we have that where |α| 2 = d 2 α|α| 2 P (α). By performing the required integrations it is easy to show that g (2) (0)− 1 = 1, i.e. g (2) (0) = 2. So far it has been shown that while a purely coherent state implies g (2) = 1 (i.e. Poissonian statistic) a thermal state implies that g (2) (0) = 2. In the case of the squeezed states it can be shown that g (2) (0) = 3 + 1 n , Q = 2 n + 1, (4.47) where n = sinh 2 r k is the multiplicity. The coherent state leads to a radiation field with Poissonian statistics. Thermal states (as well as squeezed states) have a statistics which is, according to the quantum optical terminology, superpoissonian. The latter statement is often dubbed by saying that if g (2) (0) > 1 photons are bunched while, in the opposite case (i.e. g (2) (0) < 1) the photons are said to be anti-bunched. The quantum optical language is much more effective for a mathematical description of the semi-classical limit than the usual considerations related to the limith → 0. Squeezed states are genuine quantum states with many particles. They are, in some sense, like coherent states with the crucial difference that their statistics is super-Poissonian. The possibility of scrutinizing the statistical properties of many-gravitons systems would rely on our ability of resolving single gravitons which is not even close to the present technological capabilities.

Relic graviton backgrounds: observables
In the literature relic graviton backgrounds are characterized in terms of different quantities and, in particular, the most common ones are 21 : • the power spectrum P T (k, τ ); • the spectral energy density of the relic gravitons Ω GW (k, τ ); • the spectral amplitude S h (ν, τ ).
The three listed variables can be related in different regimes. For instance the power spectrum has a simple relation to the spectral energy density when the relevant wavelengths are inside the Hubble radius. In section 6 it will be argued that, for numerical applications, the transfer function for the spectral energy density is more practical to compute than the transfer function for the power spectrum or for the spectral amplitude itself. The power spectrum is actually a strongly oscillating function of the conformal time coordinate for wavelengths shorter than the Hubble radius (i.e. kτ > 1); in the same limit the spectral energy density is asymptotically constant.

The tensor power spectrum and the spectral amplitude
The two-point function of the tensor modes of the geometry is defined as: where the state |0 is annihilated byâ k,λ for λ = ⊗, ⊕. Recalling Eqs. (2.9) and (4.3),μ ij = aĥ ij , the expansion ofĥ ij will then be: where ǫ ij (k) has been defined in Eqs. (2.10)-(2.11) and where In Eq. (5.2) F k,λ (and the associated G k,λ ) are the (complex) tensor mode functions obeying It is immediate to realize that where f k,λ and g k,λ have been introduced, for a single graviton polarization, in Eqs. (4.10)-(4.11).
After computing the expectation value we get: which, thanks to spatial isotropy, can also be written as: where P T (k, τ ) is, by definition, the tensor power spectrum: In Eq. (5.7), P T (k, τ ) is the tensor power spectrum which has been implicitly introduced in Eq. (1. 2) when talking about the customary parametrization of tensor power spectra in the context of the ΛCDM model. It is often useful, for practical applications, to consider h ij ( x, τ ) as a classical field characterized by a Fourier amplitude obeying a specific stochastic average. Then we can write where the Fourier amplitudes obey: Following the suggestions of [91], it is useful to introduce, for some applications the stochastic fields where P T (k) denotes the tensor power spectrum and where the factor 2 in front of the averages arises as a consequence of the √ 2 appearing in Eq. (5.11). In Eqs. (5.11), (5.12) and (5.13) the conformal time coordinate is absent. In Eq. (5.10) the conformal time appears explicitly. Indeed, Eqs. (5.11), (5.12) and (5.13) tacitly assume that h ⊕ ( k, τ ) = e ⊕ ( k)T e (k, τ ) and that h ⊗ ( k, τ ) = e ⊗ ( k)T e (k, τ ). This factorization is related to the concept of transfer function for the amplitude which will be discussed in section 6. The decomposition of Eqs. (5.11), (5.12) and (5.13) is useful when all the polarization have to be treated simultaneously typically in problems involving long wavelength gravitons (see Eqs. (7.85)-(7.86) and discussion therein). Furthermore the decomposition of Eqs. (5.11), (5.12) and (5.13) allows to factorize the dependence upon the initial spectrum which is useful for numerical applications.
It is finally common to characterize the stochastic backgrounds of relic gravitons in terms of the socalled spectral amplitude. The definition of the spectral amplitude can be read-off from the definition of P T (k, τ ). More precisely, taking the limit x → y in Eq. (5.6) we will have that where the second equivalence defines the spectral amplitude S h (ν) by recalling, once more, that the comoving wavenumber is related to the comoving frequency as k = 2πν.

Energy-momentum tensors for the relic gravitons
According to Eq. (3.11), each polarization of the graviton obeys the action of a minimally coupled scalar field. This observation was discussed, in particular, in [109,110] by Ford and Parker (see also, e.g. [48]). Following then Refs. [109,110] and within our set of notations the energy-momentum tensor of the relic gravitons becomes which can be obtained from the action of the gravitons by taking the functional derivative with respect to g µν . By making explicit the sum over the polarizations, Eq. (5.15) becomes: i.e., even more explicitly, where ∂ τ denotes a derivation with respect to the conformal time coordinate τ . Since, by definition, the energy density and pressure of the relic gravitons can then be written as The superscript in the energy density and pressure (i.e. ρ GW and p GW ) is convenient since different prescriptions for assigning the energy-momentum pseudo-tensor will be compared in a moment.
The other possible definition of the energy-momentum pseudo-tensor stems from the generalization, to curved space-times, of the usual flat space-time approach [111] (see also [112,113]). The nonlinear corrections to the Einstein tensor, will consist, to lowest order, of quadratic combinations of h ij that can be formally expressed as 25) where the superscript at the right hand side denotes the second-order fluctuation of the corresponding quantity while the subscript refers to the tensor nature of the fluctuations. This procedure is essentially the one described in [112,113] and has been re-explored, in a cosmological context, in [114,115] (see also [116]). Recalling the form of the Einstein tensor, we obtain ℓ 2 P T 00 = Hh ′ kℓ h kℓ + where D 00 is a total derivative, i.e.
From Eqs. (3.9) and (3.10) it is also possible to write: where D ij and D R are further total derivative Therefore, up to total derivatives, the following result holds: and To pass from doubly covariant indices to mixed ones, it is useful to recall that, to second order, By looking at the form of the specific terms arising in the previous equation it is clear that T 0 0 = g 00 T 00 and that T j i = g jk T ki . The expressions for T 0 0 and T j i are with Σ i i = 0. (see also Eqs. (5.36) and (5.37)). These expressions coincide with the ones obtained, for instance, in [114,115] and are also consistent with the ones of [112,113]. From Eqs. (5.34) and (5.36) the components of the energy and pressure density can be easily obtained since, by definition, ρ (2) gw = T 0 0 and p (2) gw = −T /3. As discussed in Eqs. (5.16)-(5.23) it is rather useful to derive the explicit form of ρ (2) gw and p (2) gw in terms of the normalized canonical amplitude defined in Eq. (3.15). The result of this calculation is simply: GW seems to be superficially different. As it will be shown in a moment the equivalence of the two approaches is clear as soon as the relevant wavelengths are larger than the Hubble radius at a given time.
Before proceeding with this step, it is relevant to remark that the components of the energymomentum pseudo-tensor given in Eqs. (5.34) and (5.35) are not covariantly conserved. However, since the Bianchi identity ∇ µ G µ ν = 0 should be valid to all orders, we will also have that: whose explicit form is Recalling now the components of the energy-momentum pseudo-tensor and the results for the fluctuations of the Christoffel symbols we have that can also be written as ∂ρ gw ∂τ + 3H(ρ gw + P gw ) = 0 (5.43) where

The energy density of the relic gravitons
By choosing the prescription of Eq. (5.19) the energy density of the relic gravitons will be given by: where the functions f k (τ ) and g k (τ ) are the mode functions obeying: If the energy momentum pseudo-tensor has the form derived in Eq. (5.34), then the energy density will have, apparently, a slightly different form: , Ω The two relevant physical limits are for modes inside the Hubble radius (i.e. kτ > 1) and for modes outside the Hubble radius (i.e. kτ < 1). If k/H > 1, then f k (τ ) will be, in the first approximation, plane waves and g k (τ ) ≃ ±ikf k (τ ). In this limit: In the limit kτ > 1 we will have H 2 ≪ k 2 . Thus Eqs. (5.50) and (5.51) coincide (up to corrections O(H 2 /k 2 )). In this limit it is also possible to express Ω GW (k, τ ) solely in terms of the power spectrum. Indeed, since we will also have (5.54) In the opposite limit, i.e. when the given wavelengths are smaller than the Hubble radius, the same analysis implies Ω GW (k, τ ). In short the idea is that, for kτ ≪ 1, g k = Hf k ; indeed, for kτ ≪ 1, Eq. (5.4) implies where the second term (going as B k /a(τ )) in Eq. (5.55) is actually negligible for large times. Therefore, in the limit kτ ≪ 1 it can be easily shown that Ω GW (k, τ ) = Ω GW (k, τ ). It is curious to notice that, up to a factor of 2, the first energy density (i.e. Eq. (5.46)) leads to the same value obtained for modes inside the Hubble radius. Therefore, for modes which are inside the Hubble radius: Consequently, we can define the critical fraction of relic gravitons at a given time as where, in the third equality it has been used that H = aH. Recalling Eq. (5.14) we can also express the critical fraction of relic gravitons in terms of the spectral amplitude. Indeed, according to Eq. (5.14) and, therefore, As long as the relevant wavelengths are shorter than the Hubble radius at a given time, different prescriptions for assigning the energy-momentum pseudo-tensor lead to the same result (see also the discussion in section 6). In the opposite limit different choices may exhibit quantitative differences. The limit of short wavelengths in comparison with the Hubble radius is the relevant one when discussing wide-band interferometers. Conversely, the initial conditions for the CMB anisotropies are set when the relevant wavelengths are larger than the Hubble radius before equality.
6 Relic gravitons from the ΛCDM scenario 6.1 Inflationary power spectra During the inflationary phase, the tensor power spectrum can be easily computed by solving Eq. (5.4): .
is the Hankel function of first kind [117,118] and where ǫ = −Ḣ/H 2 . To pass from Eq. (5.4) to Eq. (6.1) it is useful to bear in mind the following pair of identities .
The second equality in Eq. (6.2) follows (after integration by parts) from the relation between cosmic and conformal times, i.e. a(τ )dτ = dt. By substituting Eq. (6.1) into Eq. (5.7) the standard expression of the tensor power spectrum can be obtained. In particular, when the relevant modes exited the Hubble radius during inflation the power spectrum becomes where the small argument limit of the Hankel functions has been taken. In the slow-roll approximation, M 2 P H 2 ≃ V /3; then Eq. (6.3) implies that 22 The spectral index defined from Eq. (6.4) is, by definition, where the second equality follows from the identities obeyed by the slow-roll parameters [119]. The tensor and scalar power spectra are customarily assigned at a reference scale (usually dubbed pivot scale): where, by definition, A T is the amplitude of the tensor power spectrum evaluated at the pivot scale k p . In the case of single-field inflationary models, the power spectrum computed at the pivot scale k p (i.e. A T ) and the spectral index n T can be related. Bearing in mind that the power spectrum of curvature perturbations is given, in single field inflationary models, as [119] the ratio between the tensor and the scalar power spectra is given by Equation (6.8) implies, recalling Eq. (6.5), that r T = −8n T . Since there is a direct relation of the tensor spectral index to r T , the number of the parameters can be reduced from two to one. In Tabs. 1 and 2 the values of r T have been reported as they can be estimated in few different analyses of the cosmological data sets. 22 Within the present notations, as already established in Eq. (3.12), ℓ P = √ 8πG = 1/M P = √ 8πM P .
6.2 Transfer functions for inflationary power spectra Equation (6.6) correctly parametrizes the spectrum only when the relevant wavelengths are larger than the Hubble radius before matter-radiation equality (i.e. kτ ≃ k/H < 1 for τ < τ eq ). To transfer the spectrum inside the Hubble radius the procedure is to integrate numerically Eq. (5.4) (as well as Eqs. (2.2)-(2.4)) across the relevant transitions of the background geometry. While the geometry passes from inflation to radiation, Eq. (6.6) implies that the tensor mode function is constant if the wavelength associated with the given Fourier mode is larger than the Hubble radius at the corresponding epoch: The term proportional to B k in Eq. (6.9) leads to a decaying mode and F k (τ ) is therefore determined, for |kτ | ≪ 1, by the first term whose squared modulus coincides with the spectrum computed in Eq. (6.4) and parametrized as in Eq. (6.6). The latter statement is true if the inflationary phase is suddenly followed by the radiation dominated epoch since, during radiation, a(τ ) ≃ τ . The situation can be different if, after inflation, a stiff age takes over. The stiffer case still compatible with causality (i.e. w t = 1) leads approximately to a scale factor a(τ ) ≃ √ τ . In the latter situation the second term in Eq. (6.9) grows logarithmically and cannot be neglected in comparison with the constant contribution. The evolution of the background (i.e. Eqs. (2.2)-(2.4)) and of the tensor mode functions (i.e. Eq. (5.4)) should therefore be solved across the radiation matter transition and the usual approach is to compute the transfer function for the amplitude [120] i.e.
In Eq. (6.10), F k (τ ) denotes the approximate form of the mode function (holding during the matterdominated phase); F k (τ ) denotes, instead, the solution obtained by fully numerical methods. The averages appearing in Eq. (6.10) refer to the average over the oscillations: as the wavelengths are inside the Hubble radius, the solutions are all oscillating 23 . The calculation of T h (k) requires a careful matching over the phases between the numerical and the approximate (but analytic) solution. After matter-radiation equality, the scale factor is going, approximately, as a(τ ) ≃ τ 2 and, therefore, the (approximate) solution of Eq. (5.4) is given by which is constant for kτ < 1. The result of the numerical integration for the amplitude transfer function T 2 h (k/k eq ) has been recently revisited in [40,41] and it turns out to be where The latter result agrees with the findings of [120] who obtain c 1 = 1.34 and b 1 = 2.50. The value of k eq can be obtained directly from the experimental data (see, for instance, last column of Tabs. 1 and 2 implying k eq ≃ O(0.009) Mpc −1 ). The WMAP 5-yr data combined with the supernova data and with the large-scale structure data would give k eq = 0.00999 +0.00028 −0.00027 Mpc −1 . A rather good analytic estimate for k eq is k eq = 0.0082879 where the typical value selected for h 2 0 Ω R0 is given by the sum of the photon component (i.e. h 2 0 Ω γ0 = 2.47 × 10 −5 ) and of the neutrino component (i.e. h 2 0 Ω γ0 = 1.68 × 10 −5 ): the neutrinos, consistently with the ΛCDM paradigm, are taken to be massless and their (present) kinetic temperature is just a factor (4/11) 1/3 smaller than the (present) photon temperature. Equation (6.14) stems from the observation that the exact solution of Eqs. (2.2)-(2.4) for the matter-radiation transition can be given as a(τ ) = a eq [y 2 + 2y] where y = τ /τ 1 . The time-scale τ 1 = τ eq ( √ 2 + 1) is related to the equality time τ eq which can be estimated as Mpc. (6.15) In the case of the WMAP 5-yr data combined with the supernova and large-scale structure data h 2 0 Ω M0 = 0.1368 0.0038 −0.0037 . Consequently, Eqs. (6.10), (6.11) and (6.12) imply that the spectrum of the tensor modes is given, at the present time, as Within the standard approach, Eq. (6.16) is customarily connected to the spectral energy density of the relic gravitons. In [40,41] it has been observed that it is simpler and more accurate to compute directly the transfer function for the spectral energy density. In the following subsection this procedure will be illustrated in two different cases.

Transfer function for the spectral energy
Instead of computing the transfer function for the power spectrum it is more direct (and more accurate) to compute directly the spectral energy density and its related transfer function. As previously discussed there can be ambiguities in assigning the energy density because of possible different forms of the energy-momentum pseudo-tensor. In particular in Eqs. (5.46) and (5.48) two different expressions have been proposed on the basis of slightly different physical considerations. The result of the numerical calculation are reported in Fig. 5 in terms of ∆ (1) ρ (κ, x) (derived from Eq. (5.46)) and in terms of the transfer function of the spectral energy density (denoted by T ρ (κ)). The quantities ∆ (1) ρ (κ, x) (and, analogously, ∆ (2) ρ (κ, x)) are defined as As a function of x = kτ and κ = k/k eq , both ∆ (1) ρ (κ, x) and ∆ (2) ρ (κ, x) go to a constant value when the relevant modes are evaluated deep inside the Hubble radius. This occurrence allows to introduce the energy transfer function which is defined as: The specific form of the energy-momentum tensor is immaterial for the determination of T 2 ρ (κ): different forms of the energy-momentum tensor of the relic gravitons will lead to the same result. This aspect can be appreciated by looking at Fig. 6 where ∆ (1,2) ρ (k, τ ) has been reported for κ = 10 −2 (plot at the left) and for κ = 10 −4 (plot at the right). The dashed and the dot-dashed curves (in both plots) correspond, respectively, to ∆ (1) ρ (κ, x) and to ∆ (2) ρ (κ, x). The full line, in both plots, corresponds to the combination k 2 |f k (τ )| 2 + |g k (τ )| 2 = k(|c + (k)| 2 + |c − (k)| 2 ), (6.21) where c ± (k) are the mixing coefficients which parametrize the solution for the tensor mode functions when the relevant wavelengths are, asymptotically, inside the Hubble radius, i.e.
In Eq. (6.22) f k (τ ) and g k (τ ) are the solutions to leading order in the limit kτ ≫ 1. From Eq. (6.22), c ± (k) are given by Using Eqs. (6.22)-(6.23), Eqs. (6.17)-(6.18) can be directly calculated in the limit x = kτ ≫ 1 with the result that which proofs that the oscillating contributions are suppressed as x −1 f for x f ≫ 1. To get to the results illustrated in Figs. 5 and 6 the evolution equations of the mode functions have been integrated by setting initial conditions deep outside the Hubble radius (i.e. x = kτ ≪ 1), by following the corresponding quantities through the Hubble crossing (i.e. x ≃ 1) and then, finally, deep inside the Hubble radius (i.e. x ≫ 1). The integration of the mode functions is most easily performed in terms off κ (x) = kf k (τ ). Usingf κ (x), ∆ ρ (κ, x) ≡ ∆ (1) ρ (κ, x) can be written as In the case of the conventional ΛCDM scenario, the Universe is suddenly dominated by radiation as soon as inflation ends. Equation (5.4) implies that F k is constant when kτ ≪ 1. The initial conditions are fixed by requiring that, at the initial time of the numerical integration, In Figs. 5 and 5, x i = 10 −5 even if, for practical reasons, the scale on the horizontal axis has been narrowed. Bearing in mind Eq. (6.15), we can also write, for x i ≪ 1, To avoid unnecessary complications, the initial condition of the integrations illustrated in Figs. 5 and 6 have been set asf κ (x i ) = x i , i.e. the initial spectrum has been rescaled. The transfer function, by definition, must always depend only on the dynamics of the transition and not upon the features (e.g slope, amplitude) of the initial power spectrum.
In the plot at the right of Fig. 5, the fit to the energy transfer function is reported with the full (thin) line on top of the diamonds defining the numerical points. The analytical form of the fit can then be written as: T ρ (k/k eq ) = 1 + c 2 k eq k + b 2 k eq k 2 , c 2 = 0.5238, b 2 = 0.3537. (6.28) Equation (6.28) permits the accurate evaluation of the spectral energy density of relic gravitons, for instance, in the minimal version of the ΛCDM paradigm.
Yet another relevant physical situation for the present considerations is the one where the background geometry, after inflation, transits from a stiff epoch to the ordinary radiation-dominated epoch. In the primeval plasma, stiff phases can arise for various reasons. Zeldovich [122] (see also [123]) suggested this possibility in connection with the entropy problem. In [75,76,77,68] it has been suggested that the stiff phase could take place after the inflationary phase with the main purpose of identifying a potential source of high-frequency gravitons. This possibility was also prompted by a possible post-inflationary dominance of a quintessence field.
The simplest consideration leading to the possibility of a post-inflationary phase stiffer than radiation is connected with our extreme ignorance of the thermal history of the Universe after inflation. In a model-independent approach, it is plausible to think that the onset of the radiation-dominance could be delayed. This may happen, in concrete models, for various reasons. One possibility could be that the inflaton field does not decay but rather changes its dynamical nature by acting as quintessence field [78] (see also [126]). In this kind of situations we are the geometry passes from a stiff phase where to a radiation-dominated phase where c st = 1/ √ 3. Note that, according to Eqs. (6.29) and (6.30), c 2 st = w t iff the (total) barotropic index is constant in time. In the limiting case w t = 1 = c 2 st and the speed of sound coincides with the speed of light. As argued in [124], barotropic indices w t > 1 would not be compatible with causality (see, however, [125]). The presence of a suitable stiff phase has been also discussed recently as an effective way of suppressing entropic fluctuations [127] which are observationally constrained by the WMAP 5-yr data.
As in the case of the matter-radiation transition the transfer function only depends upon κ which is defined, this time, as κ = k/k s , where k s = τ −1 s and τ s parametrizes the transition time. A simple analytical form of the transition regime is given by a(y) = a s y 2 + 2y, where, by definition, ρ si = ρ s (τ i ) and ρ Ri = ρ R (τ i ). Equation (6.31) is a solution of Eqs. (2.2)-(2.4) when the radiation is present together with a stiff component which has, in the case of Eq. (6.31) a sound speed which equals the speed of light. In the limit y → 0 the scale factor expands as a(y) = √ 2y while, in the opposite limit, a(y) ≃ y. In Fig. 7 (plot at the right) ∆ ρ (κ, x) is illustrated for different values of κ. We shall not dwell here (again) about the possible different forms of the energy momentum pseudo-tensor: provided the energy density is evaluated deep inside the Hubble radius the different approaches to the energy density of the relic gravitons give the same result. The transfer function for the spectral energy density is numerically illustrated always in Fig. 7 (plot at the left). The semi-analytical form of the transfer function becomes, this time, where k s = τ −1 s . The value of k s can be computed in an explicit model 24 but it can also be left as a free parameter. Taking into account that the energy density of the inflaton will be exactly ρ si ≃ H 2 i M 2 P , the value of k s (as well as the duration of the stiff phase) will be determined, grossly speaking, by H i /M P . In Fig. (7) (plot at the left) the full line superimposed to the numerical points (illustrated

Analytic results for the mixing coefficients
The analytic results for the mixing coefficients are rather useful to obtain the final expression of the various transfer functions. Indeed, defining as k * the typical wavenumber of the transition (e.g. k * = k eq in the case of the radiation-mattter transition), the slope of the transfer function of the spectral energy density can analytically obtained in the limit κ ≫ 1. This observation helps when we have, for instance, to fit the numerical data points with an analytical expression which will however reproduce the data not only for κ > 1 (as Figs. 5 and 7 clearly show).
For illustration of the method it is practical to consider the transition from a generic accelerated phase to a decelerated stage of expansion. In this situation, by naming the transition point −τ 1 , the continuous and differentiable form of the scale factors can be written as: where the scale factors are continuous and differentiable at the transition point which has been generically indicated as τ 1 . The pump fields of the tensor mode functions turn out to be: The mode functions can then be written as: The continuity of the tensor mode functions at the transition point implies that the mixing coefficients are given by: where, according to the notations previously established, y 1 = y(−τ 1 ) = (α/β)x 1 . The case α = β = 1 corresponds to a transition from the inflationary phase to a radiation-dominated phase. In this case we do know which are the mixing coefficients. The previous expressions give us: which clearly agree with previous results [44,48,50]. In the case of Eq. (6.40) |c + (k)| 2 − |c − (k)| 2 = 1 and k 4 |c − (k)| 2 is exactly scale-invariant. Another interesting situation is the one of the transition from inflation to stiff, i.e. β = 1, α = 1/2, y 1 = x 1 /2 which leads to a logarithmic enhancement at small wavenumbers [75,76]. In this situation the mixing coefficients can be written as: 0 (x 1 /2) 0 (x 1 /2)] , (6.41) The above result can be expanded in for x 1 ≪ 1 and the result is:

Exponential damping of the mixing coefficients
Consider the case of the ΛCDM paradigm where the inflationary epoch is almost suddenly followed by the radiation-dominated phase. By denoting the transition time as τ i , it is plausible to think that all the modes of the field such that k > a i H i ≃ τ −1 i are exponentially suppressed [129,130]. For the modes kτ i > 1, the pumping action of the gravitational field is practically absent. There will be a given k, be it k max , for which k ≃ τ −1 i . The latter wavenumber is, in practice, the maximal k to be amplified and it can be estimated as : By taking as typical values of the curvature perturbations at the pivot scale the one endorsed, for instance, by the WMAP 5-yr data alone we will have we will have that Eqs. (6.45) and (6.46) can be written as: where the typical values of the slow-roll parameter have been derived by taking into account that, in the absence of running of the tensor spectral index, r T = 16ǫ; since, according to the WMAP 5-yr data alone, r T < 0.43, ǫ ≤ 0.01. For phenomenological purposes it can be also interesting to know what kind of exponential suppression we can expect. From the analysis of various transitions it emerges that the mixing coefficients for k > k max (or ν > ν max ) will satisfy |c + (k)| 2 − |c − | 2 = 1, |c + (k)| 2 + |c − (k)| 2 = e −2β k kmax + 1.   Typically, however, β > 2 for sufficiently smooth transitions. To justify this statement it is interesting to consider the following toy model where the scale factor interpolates between a quasi-de Sitter phase and a radiation-dominated phase: For τ → −∞ (i.e. τ ≪ −τ i ) , a(τ ) ≃ −a i /τ and the quasi de-Sitter dynamics is recovered. In the opposite limit (i. e. τ ≫ +τ i ), a(τ ) ≃ a i τ and the radiation dominance is recovered. In Fig. 6 (plot at the left) the exponential damping of the mixing coefficients is numerically illustrated. The curve at the top (full line) illustrates the case κ = 1. The cases κ = 2 and κ = 3 are barely distinguishable at the bottom of the plot. Notice, always in the right plot, the rather narrow range of times which are reported in a linear scale. In the plot at the right the asymptotic values of the mixing coefficients are reported for different values of κ = k/k max . By fitting the numerical data with with an equation of the form given in Eq. (6.48), the value of β = 6. 33. Different examples can be presented on the same line of the one discussed in Fig. 6. While it is pretty clear that the decay is indeed exponential, the value of β may well vary. This can be summarized, for instance, in a rescaling of k max , i.e. by positing, for instance that k max →k max /β. Thus, the dynamics of the transition can slightly shift the numerical value of the upper cut-off by a numerical factor which depends upon the width of the transition regime.

Nearly scale-invariant spectra
By using the transfer function for the tensor amplitude, the spectral energy density for frequencies ν ≫ ν eq can be simply given by: where d A (z * ) is the (comoving) angular diameter distance to decoupling. In Eqs. (6.50)-(6.51) (as well as in the program used for the numerical calculations) there are two complementary options. The first one is to use the angular diameter distance to decoupling which is directly inferred from the CMB data. For instance, in the case of the 5-yr WMAP data alone, d A (z * ) = 14115 Mpc +188 −191 . This is the strategy also adopted in other studies [132] (in connection, obviously, with earlier releases of WMAP data). At the same time it is also possible to take the best fit value of the total matter fraction (i.e. Ω M0 = 0.258 for the case of the WMAP 5-yr data alone) and compute the comoving angular diameter distance according to the well know expression for spatially flat Universes: The latter strategy has been used, for instance, in [131]. The two strategies are compatible and, moreover, this explains why, in Eq. (6.51) the dependence upon Ω M0 does not cancel. In Eq. (6.50) n T denotes, as usual, the tensor spectral index which can be also written as If α T = 0, the standard case is recovered. It should be finally mentioned that, in the limit ν ≫ ν eq , the oscillating terms have to be appropriately averaged: this is done by setting the terms going as cos 2 (2πντ 0 ) to 1/2. The latter procedure has been employed, for instance, in the analyses of [132,38]. In Fig. 9 the transfer function for the spectral energy density (introduced in Eqs. (6.20) and (6.28)) has been consistently employed to estimate the spectral energy density itself and the spectral amplitude. According to Eq. (5.14) the spectral energy density can be written as where we used that k = 2πν and that Ha = H. By making explicit the numerical factors in Eq. where, we recall, H 0 = 3.24078 × 10 −18 h 0 Hz. In Fig. 9 the spectral energy of the relic gravitons as well as S h (ν, τ 0 ) are reported for different values of r T . The spectra of Fig. 9 have been obtained from the direct integration of the mode functions but can be parametrized, according to Eq. (6.28) as (6.57) By comparing Eqs. (6.50)-(6.51) to Eqs. (6.56)-(6.57), the amplitude for ν ≫ ν eq differs, roughly, by a factor 2. This coincidence is not surprising since Eqs. (6.50)-(6.51) have been obtained by averaging over the oscillations (i.e. by replacing cosine squared with 1/2) and by imposing that |g k | = k|f k |. These manipulations are certainly less accurate than the procedure used to derive the transfer function for the spectral energy density.
So far the evolution of the tensor modes has been treated in the absence of anisotropic stress. This approximation is, strictly speaking, unrealistic. Indeed we do know that there are sources of anisotropic stress. After neutrino decoupling, the neutrinos free stream and the effective energy-momentum tensor acquires, to first-order in the amplitude of the plasma fluctuations, an anisotropic stress, i.e.
The presence of the anisotropic stress clearly affects the evolution of the tensor modes. To obtain the wanted equation we perturb the Einstein equations to first-order and we get: Equation (6.59) reduces to an integro-differential equation which has been analyzed in [36] (see also [37,38,39]). The overall effect of collisionless particles is a reduction of the spectral energy density of the relic gravitons. Assuming that the only collisionless species in the thermal history of the Universe are the neutrinos, the amount of suppression can be parametrized by the function where R ν is the fraction of neutrinos in the radiation plasma. In the case R ν = 0 (i.e. in the absence of collisionless patrticles) there is no suppression. If, on the contrary, R ν = 0 the suppression can even reach one order of magnitude. In the case N ν = 3, R ν = 0.405 and the suppression of the spectral energy density is proportional to F 2 (0.405) = 0.645. This suppression will be effective for relatively small frequencies which are larger than ν eq and smaller than the frequency corresponding to the Hubble radius at the time of big-bang nucleosynthesis, i.e. Hz. (6.61) The effect of neutrino free-streaming (as well as all the other late time effects addressed in this subsection) has been taken into account in Fig. 2 (see section 1) but it is absent in Fig. 9. The second effect which has been taken into account in Fig. 2 is the damping associated with the (present) dominance of the dark energy component. Indeed the redshift of Λ-dominance is simply defined by Consider now the mode which will be denoted as k Λ , i.e. the mode reentering the Hubble radius at τ Λ . By definition k Λ = H Λ a Λ must hold. But for τ > τ Λ is constant, i.e. H Λ ≡ H 0 where H 0 is the present value of the Hubble rate. Using now Eq. (6.62), it can be easily shown that k Λ = (Ω M 0 /Ω Λ ) 1/3 k H where k H = a 0 H 0 . The frequency interval between ν H and ν Λ is rather tiny. Indeed, it turns out that Hz. (6.64) For the same choice of parameters of Eq. (6.64), ν H = 3.708 × 10 −19 Hz which is not so different than ν Λ = 2.607 × 10 −19 Hz. It would therefore seem that this branch of the spectrum could be easily neglected. However, it turns out that the adiabatic damping of the mode function across τ Λ reduces the amplitude of the spectral energy density by a factor (Ω M0 /Ω Λ ) 2 . For the typical choice of parameters of Eqs. (6.63) and (6.64) we have that the suppression is of the order of 0.12. Again this is a tiny number which is, anyway, comparable with the suppression due, for instance, to the neutrino free streaming. This class of effects has been repeatedly in a number of recent papers [137] (see also [145]). The essence of the effect is captured by the following observation. Consider a mode k which reenters before τ Λ . The present value of the amplitude F k (τ ) = f k (τ )/a(τ ) will be adiabatically suppressed since, as repeatedly stressed, in this regime f k (τ ) will simply be plane waves. Consequently, defining asF k * the amplitude at k * = H * a * when the given mode crosses the Hubble radius, we will also have that where the subscripts (in the first equality) denote the time range over which the corresponding redshift is computed, i.e. either matter-dominated or Λ-dominated stages. The second equality follows from the first one by appreciating that a(k * ) ≃ τ 2 * ≃ k −2 and by using Eq. (6.62). Equation (6.65) implies, immediately, that the spectral energy density of relic gravitons is corrected in two different fashions. For ν < ν H the frequency dependence will be different and will be proportional to Ω GW (ν, τ 0 ) ∝ (ν/ν H ) n T −2 (Ω M0 /Ω Λ ) 2 . Vice versa, in the range ν > ν H the frequency dependence will be exactly the one already computed but, overall, the amplitude will be smaller by a factor (Ω M0 /Ω Λ ) 2 . The modification of the frequency dependence is only effective between 0.36 aHz and 0.26 aHz: this effect is therefore unimportant and customarily ignored (see, for instance, [132,137]) for phenomenological purposes. On the other hand, the overall suppression going as (Ω M0 /Ω Λ ) 2 must be taken properly into account on the same footing of other sources of suppression of the spectral energy density. There is, in principle, a third effect which may arise and it has to do with the variation of the effective number of relativistic species. Recall, indeed, that the total energy and the total entropy densities of the plasma can be written as For temperatures much larger than the top quark mass, all the known species of the minimal standard model of particle interactions are in local thermal equilibrium, then g ρ = g s = 106.75. Below, T ≃ 175 GeV the various species start decoupling, the notion of thermal equilibrium is replaced by the notion of kinetic equilibrium. The time evolution of the number of relativistic degrees of freedom effectively changes the evolution of the Hubble rate. In principle if a given mode k reenters the Hubble radius at a temperature T k the spectral energy density of the relic gravitons is (kinematically) suppressed by a factor which can be estimated as (see, for instance, [137]) At the present time g ρ0 = 3.36 and g s0 = 3.90. In general terms the effect parametrized by Eq. (6.67) will cause a frequency-dependent suppression, i.e. a further modulation of the spectral energy density Ω GW (ν, τ 0 ). The maximal suppression one can expect can be obtained by inserting into Eq. (6.67) the largest g s and g ρ . So, in the case of the minimal standard model this would imply that the suppression (on Ω GW (ν, τ 0 )) will be of the order of 0.38. In popular supersymmetric extensions of the minimal standard models g ρ and g s can be as high as, approximately, 230. This will bring down the figure given above to 0.29.
All the three effects estimated in the last part of the present section (i.e. free streaming, dark energy, evolution of relativistic degrees of freedom) have common features. Both in the case of the neutrinos and in the case of the evolution of the relativistic degrees of freedom the potential impact of the effect could be more pronounced. For instance, suppose that, in the early Universe, the particle model has many more degrees of freedom and many more particles which can free stream, at some epoch. At the same time we can say that all the aforementioned effects decrease rather than increasing the spectral energy density. Taken singularly, each of the effects will decrease Ω GW by less than one order of magnitude. The net result of the combined effects will then be, roughly, a suppression of Ω GW (ν, τ 0 ) which is of the order of 3 × 10 −2 (for 10 −16 Hz < ν < 10 −11 Hz) and of the order of 4 × 10 −2 for ν > 10 −11 Hz. These figures are comparable with the possible inaccuracies stemming from the calculation of the transfer function and, therefore, this is a further motivation, to use the transfer function of the spectral energy density. Finally the late time effects reduce a quantity which is already pretty small, i.e., as computed, h 2 0 Ω GW (ν, τ 0 ) ≃ 10 −15 for ν ≫ ν eq .

B-modes induced by long wavelength gravitons
In the minimal realization of the ΛCDM scenario the scalar fluctuations of the geometry induce an E-mode polarization which has been observed and which is now subjected to closer scrutiny [3,4,5,6,7] The tensor modes of the geometry not only induce an E-mode polarization but also a Bmode polarization. The detected angular power spectra due to the presence of a putative (adiabatic) curvature perturbation are the temperature autocorrelation (TT angular power spectrum) the E-mode autocorrelation (EE angular power spectrum) and their cross correlation (i.e. the TE angular power spectrum). The various angular power spectra of the temperature and polarization observables have been already defined in section 2 (see, in particular, Eqs. (2.48)-(2.49) and discussions therein). Long wavelength gravitons contribute not only to the TT, EE and TE angular power spectra but also to the B-mode autocorrelations, i.e. the BB angular power spectra. The effect of long wavelength gravitons on the temperature and polarization observables can be studied by deriving the evolution equations of the brightness perturbations which are related, in loose terms, to the fluctuations of the Stokes parameters. The tensor nature of the fluctuation defined in Eq. (2.8) plays, in this respect, a decisive role. In particular the following two points should be borne in mind: • in the case of the scalar modes of the geometry the heat transfer equations have an azimuthal symmetry; • in the case of the tensor modes the fluctuations of the brightness do depend, both, upon µ = cos ϑ as well as upon ϕ; this is ultimately, the rationale for the existence of a B-mode polarization.
The heat transfer equations can be schematically written as where ǫ ′ = x eñe σ T a is the differential optical depth 25 . In the expression for the differential optical depth σ T = (8/3)πr 2 0 where r 0 = e 2 /m e is the classical radius of the electron. The (differential) cross section for Compton scattering of polarized photons can be written in terms of σ T and it takes the usual Klein-Nishina form: whereê andê ′ are, respectively, the outgoing and ingoing polarizations of the photons; E f and E i are, respectively, the energies of the outgoing and of the ingoing photons. In the limit E f ≃ E i the differential cross section becomes, as anticipated 3) The right left side of Eq. (7.1) constitutes the collisionless term while the right hand side is the collisional contribution. At the right-hand side of Eq. (7.3) M(Ω, Ω ′ ) is, in general, a matrix whose dimensionality depends upon the specific problem. As it will be shown M(Ω, Ω ′ ) can be easily computed from Eq. (7.3). In similar terms f (Ω) should be understood as a column matrix whose components are the various Stokes parameters.

Collisionless Boltzmann equation for the tensor modes
The collisionless part of the Boltzmann equation can be written as: where q i is the comoving three-momentum of the photon. If the tensor fluctuation of the geometry is parametrized as in Eq. (2.8) δ (1) t g ij = −a 2 (τ )h ij ( x, τ ), (7.5) then, the right hand side of Eq. (7.4) can be written as where the prime denotes a derivation with respect to τ and where the following identities have been used In Eq. (7.7)n i denotes the direction of the photon momentum. The third identity appearing in Eq. (7.7) stems directly from the definition of comoving three-momentum, i.e.
where P i is a generic spatial component of the canonical momentum obeying Using Eq. (2.8) the full metric g ij appearing in Eq. (7.8) becomes Consequently the comoving three-momentum q i and its conformal time derivative can be expressed as: To obtain an explicit expression for dq i /dτ the conformal time derivative of the canonical momentum should be made explicit by using the geodesic equation, i.e.
To first order in the tensor fluctuations of the geometry the Christoffel connection will then lead to the following expression: (7.14) Using Eq. (7.14) into Eqs. (7.12) and (7.13), the third identity reported in Eq. (7.7) is quickly recovered. Depending upon the problem at hand the collisionless part of the Boltzmann equation can be written in different (but equivalent) forms. Equation (2.9) allows to write h ′ ij in terms of the two polarizations: Using Eqs. (7.7) and (7.15), Eq. (7.6) can be written, in Fourier space, as where C coll denotes, for completeness, the collisional part of the Boltzmann equation; µ =k ·n is the projection of the photon momentum on the direction of the Fourier mode. If the direction of propagation of the gravitational wave coincides with the third (Cartesian) component, i.e.k =ẑ. It is then easy to see that the unit vectorsâ andb must be directed, respectively, along thex andŷ Cartesian directions. This particular choice of the coordinate system implies then

(7.39)
M(Ω, Ω ′ ) will be a 3 × 3 matrix whose entries can be computed in full analogy with what has been already done in Eqs. (7.33)-(7.36). In general terms, the ingoing electric field can be expressed as implying that the outgoing Stokes parameters can be written as ]. (7.43) where, for the moment, the collisionless part of the Boltzmann equation has been neglected. Inserting Eq. (7.40) into Eqs. (7.41)-(7.43) the explicit form of M(Ω, Ω ′ ) since the Stokes parameters of the ingoing radiation field are nothing but The terms M 11 , M 12 , M 21 , M 22 will be exactly the ones already evaluated in Eqs. (7.33)-(7.36). The remaining entries are found to be  • for the same transformation, the remaining entries flip their respective sign.

Different parametrizations of the full equation
The phase space distribution obeying the collisionless part of the Boltzmann equation (see, e. g., Eq. (7.16)) does depend upon k = | k| (i.e. the Fourier mode). In the previous subsection the dependence upon k has been suppressed (for sake of simplicity) since all the aforementioned considerations could be separately repeated for each Fourier mode. Similarly, from Eq. (7.16), it is clear that the phasespace distribution does depend not only uponn (i.e. the direction of the photon) but also upon its comoving three-momentum. The phase space distribution consists of 3 independent functions whose dependence can be written, in Fourier space, as: f (k, τ, µ, ϕ, q) = [I x (k, τ, µ, ϕ, q), I y (k, τ, µ, ϕ, q), U(k, τ, µ, ϕ, q)]. (7.50) As already discussed it is more handy to use directly µ = cos ϑ as independent variable. The full Boltzmann equation can then be solved by iteration around the equilibrium configuration provided by the Bose-Einstein distribution f 0 (q); this means that the 3 components of f (τ, µ, ϕ) can be written as: I y (k, τ, µ, ϕ, q) = f 0 (q)[1 + I y (k, τ, µ, ϕ)], (7.52) U(k, τ, µ, ϕ, q) = f 0 (q)U(k, τ, µ, ϕ).
The final result of the previous pair of manipulations can be explicitly written as 26 where Ψ(k, τ ) is given by Concerning the result obtained in Eqs. (7.77) and (7.78) few comments are in order: • in Eqs. (7.77)-(7.78) h ′ λ denotes indifferently either h ′ ⊕ or h ′ ⊗ : the derivation reported in the case of h ⊕ can be repeated in the case of h ⊗ bearing in mind the differences in the angular structure (i.e. Eqs. (7.64)-(7.66) should be used instead of Eqs. (7.61)-(7.63)); • by changing ϕ → −ϕ and ϕ ′ → −ϕ ′ , Eqs. (7.61)-(7.63) and Eqs. (7.64)-(7.66) will be different because if the various sines appearing the various expressions; therefore the angular structure of the Stokes parameters will change but Eqs. (7.77) and (7.78) will keep their form.
The rationale for the latter statement is that while ϕ → −ϕ and ϕ ′ → −ϕ ′ changes the angular structure of the Stokes parameters, also some of the entries of M(Ω, Ω ′ ) will change. The net result, as already mentioned, will be that Eqs. (7.77) and (7.78) will still be valid. The result expressed by Eqs. (7.77) and (7.78) has been firstly obtained by Polnarev [133] and then exploited for different 26 As in the previous sections the h ′ ⊕ and h ′ ⊗ denote the time derivatives of the polarizations with respect to the conformal time coordinate τ . This notation has been avoided in the previous equations of the present section since it could have been confused with the angular variables describing the polarizations of the outgoing photons. From now on this possible clash of notations does not arise.
semi-analytical discussions of the problem (see, e. g. [134,135,136]). The polarization decomposition leading to Eqs. (7.77) and (7.78) can be related to slightly different treatments (see, for instance, [91]) using the brightness perturbations rather than β and ζ. The connection between the different formalisms will be explored later in this section. By summing up Eqs. (7.77) and (7.78) the collision terms cancel: Defining ξ = β + ζ the collision term of Eq. (7.76) can also be written as: The collision term can be further simplified by expanding in multipoles the unknowns β(k, µ, τ ) and ζ(k, µ, τ ) where P ℓ (µ) denote, as usual, the Legendre polynomials. By recalling that the first three Legendre polynomials of even order are the orthogonality of the Legendre polynomials imply that the integral over x in Eq. (7.80) leads to So far the discussion has been conducted in terms of one single polarization at the time. For both polarizations the evolution equations of ζ and β turn out to be, of course, the same. However the azimuthal structure of the relevant Stokes parameters will be different for different poolarizations. In view of specific applications it is useful to introduce a unified notations allowing for the simultaneous treatment of the two different cases: I( k, τ, µ, ϕ) = I x ( k, τ, µ, ϕ) + I y ( k, τ, µ, ϕ) where, following Eqs. (5.11)-(5.13) (see also [91]) the stochastic variables terms of an appropriate transfer function for the amplitude. By linearly combining Eqs. (7.85) and (7.86) it is also easy to obtain Q( k, τ, µ, ϕ) + iU ( k, τ, µ, ϕ) = (1 + µ) 2 e 2iϕ e 1 ( k) + (1 − µ) 2 e −2iϕ e 2 ( k) β(k, µ, τ ) √ 2 , (7.88) The dependence of Eqs. (7.78) and (7.79) on the derivative of f 0 (q) with respect to the comoving three-momentum q can be further simplified by rephrasing the Boltzmann equations in terms of the appropriate brightness perturbations. Recalling that, by definition of brightness perturbations, the intensity I introduced in Eq. (7.84) can be written as I = −∆ I (∂ ln f 0 /∂ ln q). Consequently, Eqs. (7.84) and (7.88)-(7.89) can also be written as: where the relation of ∆ T (k, τ, µ) and ∆ P (k, τ, µ) to β(k, τ, µ) and ζ(k, τ, µ) is Using Eq. (7.94) inside Eqs. (7.78)-(7.79), the evolution equations obeyed by ∆ P (k, τ, µ) and ∆ T (k, τ, µ) can be obtained and it is: In Eq. (7.96) h is the canonical amplitude of the graviton. Recalling the results of section 3 it has been defined that h ⊕ = √ 2ℓ P h and h ⊗ = √ 2ℓ P h. The √ 2 factor simplifies once the brightness perturbations of Eqs. (7.94) are used. Using the common practice Eqs. (7.95)-(7.97) have been written in units ℓ P = 1.
After Eqs. (7.19) and (7.20) it has been stressed that by flipping the sign of ϕ and ϕ ′ the outgoing and ingoing polarizations change but the choice of orthnormal frame is still valid. The discussion conducted so far can be repeated by assuming the polarizations of Eqs. (7.19) and (7.20) but with ϕ → −ϕ and ϕ ′ → −ϕ ′ . The same sign flip can be followed throughout the derivation and it can be appreciated that the results of Eqs. (7.95), (7.96) and (7.97) are left unchanged (see also discussion after Eqs. (7.48)-(7.49)). The flip in the sign of ϕ and ϕ ′ affects, however, the brightness perturbations. The intensity of the radiation field, as expected, is left invariant when ϕ → −ϕ (see Eq. (7.91)).

B-modes from relic gravitons
The analytic and numerical solutions of the evolution equations for the brightness perturbations allow for an explicit evaluation of the temperature and polarization observables. The results obtained so far imply that long wavelength gravitons will produce both temperature and polarization anisotropies. More specifically, following the terminology of section 2 (see, in particular, Eqs. (2.48)-(2.50)) the relevant angular power spectra induced by the relic gravitons will be the TT, EE, TE and BB angular power spectra.
Bearing in mind the connection between spherical harmonics and the associated Legendre functions, the Y ℓ m (n) will be essentially given by the product of P m ℓ (µ) and of e imϕ , up to a well know normalization coefficient. In Eq. (7.108) the integration over ϕ can be carried on in explicit terms and since there is a sin 2ϕ in the integrand the whole integral will be proportional to (δ m, 2 − δ m −2 ); in more precise terms the integration over ϕ will bring Eq. (7.108) in the form: (7.109) Equation (7.109) can be further simplified by using the following three relations: Equations (7.111) and (7.112) are well known recurrence relations for the for the associated Legendre functions (see p. 24 of Ref. [86] for a derivation). Equation (7.110) can be written in slightly different ways. Some authors do not include, at the right hand side of Eq. (7.110), the factor (−i) j inside the sum (see for instance [93]). The analog of Eq. (7.109) can also be obtained within the approach of Ref. [93] whose conventions (for the ⊕ polarization) are For × polarization, the angular structure of Eq. (7.114) will be different and can be obtained from the results for the ⊕ polarization by replacing cos 2ϕ → sin 2ϕ and sin 2ϕ → − cos 2ϕ. Eqs. (7.113) and (7.114) are nothing but the explicit expression of the polarization tensor on the sphere introduced in Eq. (2.81) with the difference that the indices (in Eq. (2.81)) are both covariant while in Eq. (7.113) they are both contravariant. The conventions of [93] imply also a difference of a factor of 1/8 in comparison with the conventions adopted here. Equations (7.111) and (7.112) can be used inside Eq. (7.109). As an example the term (1 − µ 2 )µ∆ ′′ P , after using Eq. (7.110) will produce a factor µ(1 − µ 2 )P ′′ j , i.e. (1 − µ 2 )∆ ′ P will produce a factor (1 − µ 2 )P ′ j , i.e.
The final result for the B-mode will then be: ]. (7.118) If the factor (−i) j would be absent from the expansion of Eq. (7.110) there would not be (−i) ℓ at the right-hand side of Eq. (7.118); for the same reason, the relative sign of the two terms inside the squared bracket would be plus instead of minus. Finally, the factor (1/8) appearing at the right-hand side of Eq. (7.113) will also modify slightly the prefactor of Eq. (7.118) which will be ℓ m . In the case of the E-mode, always working with the ⊕ polarization, Eqs. (7.104) and (7.107) will give, after integration over ϕ, (7.120) where Eqs. (2.54) and (2.55) have been used. To complete the calculation it is necessary to perform the integration over µ. This can be done, as previously shown, by using wisely the recursion relations of the associate Legendre functions. From Eq. (7.120) we have that, using the expansion of Eq. (7.110), As in the case of the B-mode, the main idea of the calculation is to express the relevant part of the integrand appearing at the right hand side of Eq. (7.121) in terms of products of associated Legendre functions with the same m but different values of j so that the orthogonality relation (7.117) can be exploited. Using the definition of the associated Legendre functions of Eq. (2.55) it is easy to show that (7.122) where the minus sign in the second term at the right hand side of Eq. (7.122) arises since the Condon-Shortley phase has been included in the definition of the associated Legendre functions (see Eqs. (2.54)-(2.55) and discussion therein). The right hand side of Eq. (7.122) now be simplified by using the following steps: • Eq. (7.111) can be used to simplify the term µ 2 P (2) j (µ); • Eq. (7.112) can be used to simplify the term 8µ • the term (1 − µ 2 )P j (µ) can be simplified by using the equation of the Legendre polynomials and the recursion relations (7.111) and (7.112).
Using these three steps, a rather lengthy but straightforward algebra leads to the following expression: j−2 (µ). (7.123) To get to Eq. (7.123) it is useful to recall, on a side, that , (7.124) which can be derived by recalling the equation obeyed by the Legendre polynomials P j (µ) as well as the recursion relations reported in Eqs. (7.111)-(7.112). Using Eq. (7.124) inside Eqs. (7.120)- (7.121) it is easy to obtain, after simple algebra, that As in the case of the B-mode the presence of the factor (−i) j in Eq. (7.110) leads to a sign difference in the terms ∆ P ℓ±2 appearing in Eq. (7.125). The absolute values of the coefficients of the three terms appearing (inside the the square bracket) in Eq. (7.125) are independent on the conventions related to Eq. (7.110). The coefficient of the term ∆ P ℓ appears to be different from the homologue term reported in Eq. (4.41) of [93]. After many checks and cross-checks of the derivations leading to Eq. (7.125) it has been concluded that the difference (i.e. 6ℓ(ℓ + 1) instead of 6(ℓ − 1)(ℓ + 2) in the numerator of the coefficient of ∆ Pℓ ) is just the result of a typo and the correct expression is 6(ℓ − 1)(ℓ + 2) as in Eq. (7.125).

Angular power spectra induced by long wavelength gravitons
The tensor modes of the geometry produce both, temperature and polarization anisotropies. Consider, for simplicity, the ΛCDM parameters determined from the best fit of the WMAP 5-yr data alone. These parameters have been reported in Eq. (2.7). In Fig. 10 the TT angular power spectra stemming from the standard adiabatic mode and the temperature autocorrelations induced by the tensor modes are illustrated. In both plots of Fig. 10 the full line denotes the TT angular power spectra corresponding  Fig. 10 (plot at the right) the dashed line denotes the tensor contribution in the case r T = 1 while the dot-dashed line denotes the total TT angular power spectrum. The temperature and polarization observables can be obtained numerically (see, e.g. [91,93]) but useful analytic results do exist in the literature [133,134,135,136,137]. An essential step for both approaches is to derive the angular power spectra in more explicit terms and as a function of the solution of the heat transfer equations. The solution of Eqs. (7.95), (7.96) and (7.97) can be formally written by using the integration along the line of sight which is also common to the scalar case (see, e.g. [119], for an introduction): where x = k(τ 0 − τ ). Note that S T (k, τ ) and S P (k, τ ) are related to the source terms of, respectively, Eqs. (7.95) and (7.96), i.e.
S P (k, τ ) = −K(τ ) S(k, τ ), S T (k, τ ) = −h ′ e −ǫ(τ,τ 0 ) + K(τ ) S(k, τ ), (7.128) where S(k, τ ) is given by Eq. (7.97) and where is the visibility function expressed in terms of the differential optical depth ǫ ′ and in terms of the optical depth ǫ(τ, τ 0 ). The visibility function K(τ ) gives that probability that a photon is last scattered between τ and τ + dτ . In loose terms the visibility function is peaked around the redshift of recombination. In analytic discussions the visibility function is often approximated by means of a Gaussian profile (see, e. g., [120,137,119] and also [138,139,140,141]): As indicated in Eq. (7.130), N (σ rec ) is determined by requiring that the integral of K(τ ) over τ is normalized to 1. The WMAP data suggest a thickness (in redshift space) ∆z rec ≃ 195 ± 2 [8] which would imply that σ rec , in units of the (comoving) angular diameter distance to recombination, can be estimated as σ rec /τ 0 ≃ 1.43 × 10 −3 . When τ 0 ≫ τ rec and τ 0 ≫ σ rec the normalization appearing in Eq. (7.130) can be estimated as N (σ rec ) → σ −1 rec 2/π. In [137] it has been suggested that a better approximation, for the case of the tensor modes, is to assume that K(τ ) has two different widths, respectively, for τ < τ rec and for τ > τ rec . Recalling now Eq.  (7.131) where the second equality follows from the Fourier transform of ∆ I (n) and by using Eq. (7.98) in the obtained expression. The angular integrations can now be performed with the approach already exploited in previous subsection for the polarization observables. It is however useful to solve Eq.
In I (T) ± the integration over ϕ is trivial since the (ordinary) spherical harmonics depend upon ϕ as e imϕ (see, e. g., Eqs. (2.54) and (2.55)). Thus, as expected on the basis of the results of the previous subsection, I (T) ± ∝ 2πδ m, ±2 . The latter results allows for a simplification so that, for instance, − (ℓ, m, x) can be performed exactly in the same way. By making explicit the ensemble averages and by using Eq. (2.49), the angular power spectrum illustrated in Fig. 10 can be obtained and it is given by: where, as in Eq. (2.46), N ℓ = (ℓ − 2)!/(ℓ + 2)!. Having discussed the temperature autocorrelation induced by the long wavelength gravitons it is now the moment of discussing the polarization observables. In Fig. 11 (plot at the left) the EE angular power spectra (full line) are compared with the B-mode autocorrelations (dashed line) induced by the tensor modes in the case r T = 0.1 while the other cosmological parameters are fixed to the best-fit values of the WMAp 5-yr data alone. Always in Fig. 11 (plot at the right) the full line illustrates the B-mode autocorrelation in the case r T = 0.1 while the dashed line illustrated the EE angular power spectrum stemming from the standard adiabatic mode. In Fig. 12 (plot at the left) the BB angular power spectra are illustrated for different values of r T compatible with current bounds on r T (see Tabs. 1 and 2).
The E-mode and the B-mode angular power spectra can be obtained from Eqs. (7.99) and (7.100) by carefully following all the steps leading to the temperature autocorrelation of Eq. (7.135). The crucial difference will be, of course, that a • the evolution for ∆ P (k, µ, τ ) can then be solved in Fourier space as in Eq. (7.126).
Now Eqs. (7.142) and (7.143) can be used into Eqs. (2.48) and the angular power spectra become: where The tensor modes also induce a TE angular power spectrum (i.e. a cross-correlation between the temperature and the E-mode polarization). In Fig. 12 (plot at the right) the absolute value of the TE (tensor) angular power spectrum is compared with the corresponding angular power spectrum induced by the standard adiabatic mode. The cross-correlation between temperature and polarization can be obtained inserting Eqs. (7.133) and (7.141) into Eq. (2.49) and the result is where, ∆ Tℓ (k, τ 0 ) and ∆ Eℓ (k, τ 0 ) are given, respectively, by Eqs. (7.135)   8 High-frequency spikes in the relic graviton background In the ΛCDM paradigm long wavelength gravitons can affect the CMB polarization. As the frequency increases towards the region accessible to wide band interferometers, the ΛCDM signal can only decrease for number of independent reasons: • in the ΛCDM paradigm the spectral energy density is only approximately scale-invariant (see, e.g. Fig. 2) but the scaling violations always tend to make the spectral energy density smaller at high frequencies; • there are various secondary effects associated, for instance, with the variation of the effective number of relativistic degrees of freedom, with the neutrino anisotropic stress and with the transition to a phase dominated by the dark-energy contribution: all these features reduce the spectral energy density in different frequency regions; • the addition of supplementary scalar fields driving the inflationary phase does not change the two previous statements.
The previous conclusions are all based (either directly or indirectly) on the assumption that the thermal history of the Universe is minimal in the sense that, after inflation, the Universe soon becomes dominated by relativistic particles so that the sound speed of the plasma soon reaches values c st = 1/ √ 3. The post-inflationary thermal history might not be minimal. For instance it could happen that the transition to the radiation-dominated regime is not instantaneous [75]. More specifically it can happen that, after inflation, the sound speed of the plasma is such that c st > 1/ √ 3. When the sound speed is larger than 1/ √ 3, the fluid is said to be stiff. In the system of units used in the present paper the speed of light is such that c = 1. A natural upper limit for the sound speed is exactly 1 which is the maximally stiff fluid compatible with causality [124] (see, however, also [125]). If the thermal history of the Universe contemplates a post-inflationary phase stiffer than radiation, a spike in the relic graviton spectrum is expected at high frequencies [75] (see also [76,77] as well as Eqs. (6.29)-(6.30) and discussion therein). The possibility of having a post-inflationary phase stiffer than radiation has been also investigated in different contexts such as in [146,147,148].
In the early Universe, the dominant energy condition might be violated and this observation will also produce scaling violations in the spectral energy density [149]. If we assume the validity of the ΛCDM paradigm, a violation of the dominant energy condition implies that, during an early stage of the life of the Universe, the effective enthalpy density of the sources driving the geometry was negative and this may happen in the presence of bulk viscous stresses [149] (see also [150,151] for interesting reprises of this idea). In what follows the focus will be on the more mundane possibility that the thermal history of the plasma includes a phase where the speed of sound was close to the speed of light.
Absent any indirect tests on the thermal history of the Universe prior to the formation of light nuclear elements, it is legitimate to investigate situations where, before nucleosyntheis, the sound speed of the plasma was larger than 1/ √ 3, at most equalling the speed of light. In this plausible extension of the current cosmological paradigm, hereby dubbed Tensor-ΛCDM (i.e. TΛCDM) scenario, highfrequency gravitons are copiously produced [40,41]. Without conflicting with the bounds on the tensor to scalar ratio stemming from the combined analysis of the three standard cosmological data sets (i.e. cosmic microwave background anisotropies, large-scale structure data and observations of type Ia supenovae), the spectral energy density of the relic gravitons in the TΛCDM scenario can be potentially observable by wide-band interferometers (in their advanced version) operating in a frequency window which ranges between few Hz and few kHz.
The presence of a stiff phase increases the spectral energy density for frequencies larger than a pivotal frequency ν s which is related to the total duration of the stiff phase. If the stiff phase takes place before BBN, then ν s > 10 −2 nHz. If the stiff phase takes place for equivalent temperatures larger than 100 GeV, then ν s ≥ µHz. If the stiff phase takes place for T ≥ 100 TeV, then ν s > mHz. The frequency ν s marks the beginning of a branch of the spectrum where the tilt of the spectral energy density is blue (i.e. increasing in slope) rather than nearly scale invariant or slightly red (as it is the case in the conventional scenario).

Scaling violations
In the ΛCDM paradigm, for frequencies larger than ν eq the spectral energy density of the relic gravitons is, approximately, scale invariant (see Fig. 2). In the context of the TΛCDM scenario the approximate scale-invariance of the flat plateau is violated. This situation is illustrated, for instance, in Fig. 4. In the present subsection a more detailed account of the typical frequencies of the problem will be presented. If there is some delay between the end of inflation and the onset of radiation the maximal wavenumber of the spectrum will be given by: where α = 2/[3(w t + 1)] (with w t > 1/3) is related to the specific kind of stiff dynamics. Equation (8.1) can also be written as In the case Σ = O(1) (as it happens in the case α = 1/3 if the initial radiation is in the form of quantum fluctuations) ν max = k max /(2π) ≃ 100 GHz, more precisely: Hz. (8.4) On top of the standard parameters of the ΛCDM scenario (see, e. g., Eq. (2.7)) the minimal TΛCDM scenario demands two supplementary parameters • the frequency ν s defining the region of the spectrum at which the scaling violations take place; • the slope of the spectrum arising during the stiff phase.
The frequency ν s can be dynamically related to the frequency of the maximum and, consequently, the first parameter can be trated for Σ. As it will be shown in a moment, typical values for Σ range between 0.01 and 0.5. The slope of the spectrum during the stiff phase depends upon the total barotropic index and can therefore be traded for w t . The curvature scale H r determines k s (or ν s ), i.e. the frequency at which the spectral energy density starts increasing. Supposing that from the end of inflation there is a single stiff phase (as it is natural to assume in a minimalistic persepective) the value of k s is Using the relation of H r to Σ, Eq. (8.5) can also be written as .
and as Hz. (8.7) The quantity Σ which parametrizes the location, in the spectrum, where the scaling violations appear must be smaller than 1 or, at most, of order 1. This is what happens within specific models. For instance, if the radiation present at the end of inflation comes from amplified quantum fluctuations (i.e. Gibbons-Hawking radiation), quite generically, at the end of inflation ρ r ≃ H 4 . More specifically In Eq. (8.8) N eff is the number of species contributing to the quantum fluctuations during the quasi-de Sitter stage of expansion. In [128] (see also [76,77,78]) it has been argued that this quantity could be evaluated using a perturbative expansion valid in the limit of quasi-conformal coupling. It should be clear that N eff is conceptually different from the number of relativistic degrees of freedom g ρ . Given H and N eff the length of the stiff phase is fixed, in this case, by [78]  where we used the fact that α = 2/[3(w + 1)] and where we defined λ = N eff /(480π 2 ). Equation (8.9) implies that a i a r = λ Recalling that H r = H(a i /a r ) 1/α , we also have that Eq. (8.10) implies Using Eq. (8.11) and the definition of Σ (see Eq. (8.3)) it turns out that Σ = λ 1/4 , which is always smaller than 1 and, at most, O(1).
Instead of endorsing an explicit model by pretending to know the whole thermal history of the Universe in reasonable detail, it is more productive to keep Σ as a free parameter and to require that the scaling violations in the spectral energy density will take place before BBN. Consequently the variation of Σ, w and r T can be simultaneously bounded [40,41]. If the stiff dynamics takes place before big-bang nucleosynthesis, then ν s > ν bbn (see also Eq. (6.61)). This requirement guarantees that the stiff dynamics will be over by the time light nuclei start being formed. In a complementary  Hz. (8.14) The latter requirement would imply that the stiff age is already finished by the time the Universe has a temperature of the order of 100 TeV when, presumably, the number of relativistic degrees of freedom was much larger than in the minimal standard model (in Eq. (8.14) the typical value of g ρ is the one arising in the minimal supersymmetric extension of the standard model).
In Fig. 13 (plot at the left) the constraints on Σ, w t are illustrated. The value of Σ controls the position of the frequency at which the nearly scale-invariant slope of the spectrum will be violated. The barotropic index w t is taken to be always larger than 1/3 and with a maximal value of 1. The shaded region corresponds to the region excluded in the most constraining case, i.e. the one demanding ν s > ν bbn . Also the values of r T are bounded by the same kind of considerations. Indeed, ν s depends upon ǫ which is related to r T (see, e.g., Eq. (8.7)). Therefore, a lower bound on ν s also implies a bound in the (r T , w t ) plane.  Fig. 4. These examples will now be discussed in greater detail. The spectral energy density has been computed by using the numerical approach presented in section 6.

Spectral energy density in the minimal TΛCDM scenario
In the two examples of Fig. 13 (plot at the right) the ΛCDM parameters are fixed to the values reported in Eq. (2.7) (see [3,4,5,6,7]). The spectral index has been allowed to run, i.e. α T = 0 (see Eqs. (1.6) and (6.53)). The two supplementary parameters should be identified with the sound speed during the stiff phase (i.e. c st ) and the threshold frequency (i.e. ν s ). Besides c st and ν s , there will also be r T which controls, at once, the normalization and the slope of the low-frequency branch of the spectral energy density. At the moment wide band interferometers have sensitivities which are insufficient for cutting through the phenomenologically interesting region [63,64,65]. In the near future, however, there is the hope of a dramatic improvement of the sensitivity: even 5 or 6 orders of magnitude at least heeding the original design (see e.g. [55]) together with the recent proposals for an advanced Ligo program.
As specifically discussed in Eqs. (8.5)-(8.6) the frequency of the elbow, i.e. ν s , is fully determined by Σ (see Eq. (8.3) and discussion therein). The two supplementary parameters ν s and c st can be traded for Σ and w t as already done in Fig. 13 (plot at the left). In doing so there is also a potential advantage since, according to Eq. (8.4), Σ shifts the maximal frequency of the spectrum.
As soon as the frequency increases from the aHz up to the nHz (and even larger) the spectral energy density increases sharply in comparison with the nearly scale-invariant case where the spectral energy density was, for ν > nHz, at most O(10 −16 ). In the case of Fig. 13 the spectral energy density is clearly much larger. The accuracy in the determination of the infra-red branch of the spectrum is a condition for the correctness of the estimate of the spectral energy density of the high-frequency branch. The plots of Fig. 13 (see also Fig. 4) demonstrate that the low-frequency bounds on r T do not forbid a larger signal at higher frequencies.
A decrease of r T implies a suppression of the nearly scale-invariant plateau in the region ν eq < ν < ν s . At the same time the amplitude of the spectral energy density still increases for frequencies larger than the frequency of the elbow (i.e. ν s ). The latter trend can be simply understood since, at high frequency, the transfer function for the spectral energy density grows faster than the power spectrum of inflationary origin. For instance, in the case w t = 1 and neglecting logarithmic corrections, Ω GW (ν, τ 0 ) ∝ ν n T +1 for ν ≫ ν s . Now, recall that n T is given by Eq. (1.4). If r T → 0, the combination (n T + 1) will be much closer to 1 than in the case when, say, r T ≃ 0.3. This aspect can be observed in Fig. 4 where different values of r T have been reported. By decreasing the w t from 1 to, say, 0.6 the extension of the nearly flat plateau gets narrower. This is also a general effect which is particularly evident by comparing the two curves of Fig. 13 (plot at the right). The slope of the high-frequency branch of the graviton energy spectrum can be easily deduced with analytic methods and it turns out to to be d ln Ω GW d ln ν = 6w t − 2 3w t + 1 , ν > ν s , (8.15) up to logarithmic corrections. The result of Eq. (8.15) stems from the simultaneous integration of the background evolution equations and of the tensor mode functions according to the techniques described in section 3. The semi-analytic estimate of the slope (see [75]) agrees with the results obtained by means of the transfer function of the spectral energy density.
To conclude the discussion it is appropriate to elaborate on the interplay between the stiff spectra and the phenomenological bounds mentioned in subsection 1.6. Let us start with millisecond pulsar bound of Eq. (1.14).
Assuming the maximal growth of the spectral energy density and the minimal value of ν s , i.e. ν bbn The detectability constraint (full line in the plot at the right) stemming from the putative sensitivities of wide-band interferometers in their advanced version. The points corresponding to the spectral energy density should lie above the full lines to be potentially interesting for those instruments.
As noticed in the past [121,75], the most significant constraint on the stiff spectra stems from BBN. The models illustrated in Fig. 13 (plot at the right) are on the verge of saturating the bounds of Eqs. (1.18)- (1.19). This conclusion stems directly from the form of spectral energy density: the broad spike dominates the (total) energy density of relic gravitons which are inside the Hubble radius at the time of big bang nucleosynthesis. For consistency with the low-frequency determinations of the tensor power spectrum r T must be bounded from above according to the values reported, for instance, in Tabs. 1 and 2. Once the the value of r T has been selected, the constraints of Eqs. (1.18)-(1.19) can be imposed. From Fig. 13 a large signal is expected for ν LV ≃ 100 Hz for 0.35 < w t < 0.6. This range turns out to be compatible with the bounds of Eqs. (1.18)-(1.19). In the opposite limit (e. g. w t ≃ 1) the spike becomes narrower, the elbow frequency augments and the signal at the interferometer scale diminishes.
In Fig. 14 19)). As the scalar spectral index diminishes, the constraints are better satisfied since n s controls α T and, consequently, the frequency dependence of the tensor spectral index n T (see Eq. (1.4)) in the case α T = 0.

Detectability prospects
By lowering w t , h 2 0 Ω GW (ν, τ 0 ) increases for ν = ν LV ≃ 0.1 kHz. This trend can be inferred from Fig. 14 (plot at the right) where the spectral energy density is evaluated exactly for ν = ν LV . To be detectable by wide band interferometers the parameters of the TΛCDM must lie above the full lines. The region of low barotropic indices emerging neatly from Fig. 14, leads to spectral energy densities which are progressively flattening as w t diminishes towards 1/3. Low values of w t bring the frequency of the elbow, i.e. ν s below 10 −10 Hz which is unacceptable since it would mean that, during nucleosynthesis, the Universe was dominated by the stiff fluid. In Fig. 13 (plot at the left) the region above the full line corresponds to a range of parameters for which ν s > ν bbn : in such a range a decrease of w t demands an increase of Σ. The latter occurrence is illustrated in Fig. 15 where, at the left, w t = 0.5. The full, dashed and dot-dashed curves illustrated in Fig. 15 (plot at the left) are incompatible with the phenomenological constraints since the frequency of the elbow is systematically smaller than ν bbn . Once more, this choice of parameters would contradict the bounds of Fig. 13 and would imply that the stiff phase is not yet finished at the BBN time. In the left plot of Fig. 15 the diamonds denote a model which is compatible with BBN considerations but whose signal at the frequency of interferometers is rather small (always three orders of magnitude larger than in the case of conventional inflationary models).
The compatibility with the phenomenological constraints demands that the parameters of the TΛCDM paradigm must lie above the full lines of Fig. 13 (plot at the left). The requirements of Fig. 13 suggest, therefore, that Σ should be raised a bit. In this case the frequency of the elbow gets shifted to the right but, at the same time, the overall amplitude of the spike diminishes. The putative amplitude remains still much larger than the conventional inflationary signal reported, for instance, in Fig. 3. In Fig. 16 the spectral energy density of the relic gravitons is illustrated as a function of r T for a choice of parameters which is compatible with all the bounds applicable to the stochastic backgrounds of the relic gravitons. The three curves refer to three different frequencies, i.e. 0.1 kHz,  Figure 16: The graviton energy spectrum is illustrated, in the TΛCDM scenario, for ν = ν LV and as a function of r T (plot at the left) in the case α T = 0. In the plot at the right α T = 0. 1 kHz and 10 kHz. Indeed, if the spectrum is nearly scale-invariant (as in the case o Fig. 3) we can compare the potential signal with the central frequency of the window. If the signal increases with frequency it is interesting to plot the same curve for some significant frequencies inside the window of wide-band interferometers. Even if the frequency window extends from few Hz to 10 kHz the maximal sensitivity is in the central region and depends upon various important factors which will now be briefly discussed.
To illustrate more quantitatively this point we remind the expression of the signal-to-noise ratio (SNR) in the context of optimal processing required for the detection of stochastic backgrounds:  (F depends upon the geometry of the two detectors and in the case of the correlation between two interferometers F = 2/5; T is the observation time). In Eq. (8.17), S (k) n (f ) is the (one-sided) noise power spectrum (NPS) of the k-th (k = 1, 2) detector. The NPS contains the important informations concerning the noise sources (in broad terms seismic, thermal and shot noises) while γ(ν) is the overlap reduction function which is determined by the relative locations and orientations of the two detectors. In [68] Eq. (8.17) has been used to assess the detectability prospects of gravitons coming from a specific model of stiff evolution with w t = 1. At that time the various suppressions of the lowfrequency amplitude as well as the free-streaming effects were not taken into account. Furthermore, the evaluation of the energy transfer function was obtained, in [70], not numerically but by matching of the relevant solutions. We do know, by direct comparison, that such a procedure is justified but intrinsically less accurate than the one proposed here. It would be interesting to apply Eq. (8.17) for the (more accurate) assessment of the sensitivities of different instruments to a potential signal stemming from the stiff age 28 . Equation (8.17) assumes that the intrinsic noises of the detectors are stationary, Gaussian, uncorrelated, much larger in amplitude than the gravitational strain, and statistically independent on the strain itself [71,72,73,74]. The integral appearing in Eq. (8.17) extends over all the frequencies. However, the noise power spectra of the detectors are defined in a frequency interval ranging from few Hz to 10 kHz. In the latter window, for very small frequencies the seismic disturbances are the dominant source of noise. For intermediate and high frequencies the dominant sources of noise are, respectively, thermal and electronic (i.e. shot) noises. The wideness of the band is very important when cross-correlating two detectors: typically the minimal detectable h 2 0 Ω GW will become smaller (i.e. the sensitivity will increase) by a factor 1/ √ ∆νT where ∆ν is the bandwidth and T , as already mentioned, is the observation time. Naively, if the minimal detectable signal (by one detector ) is h 2 0 Ω GW ≃ 10 −5 , then the cross-correlation of two identical detector with overlap reduction γ(ν) = 1 will detect h 2 0 Ω GW ≃ 10 −10 provided ∆ν ≃ 100 Hz and T ≃ O(1yr) (recall that 1yr = 3.15 × 10 7 Hz −1 ). The achievable sensitivity of a pair of wide band interferometers crucially depends upon the spectral slope of the theoretical energy spectrum in the operating window of the detectors. So, a flat spectrum will lead to an experimental sensitivity which might not be similar to the sensitivity achievable in the case of a blue or violet spectra. Previous calculations [68,69,70] showed that, however, to get a reasonable idea of the potential signal it is sufficient to compare the signal with the sensitivity to flat spectrum which has been reported in Eq. (1.13). Of course any experimental improvement in comparison with the values of Eq. (1.13) will widen the detectability region by making the prospects of the whole discussion more rosy.
In the TΛCDM paradigm the maximal signal occurs in a frequency region between the MHz and the GHz. This intriguing aspect led to the suggestion [68,69] that microwave cavities [152] can be used as GW detectors precisely in the mentioned frequency range. Prototypes of these detectors [153] have been described and the possibility of further improvements in their sensitivity received recently attention [154,155,156,157,158,159]. Different groups are now concerned with high-frequency gravitons. In [155] the ideas put forward in [152,153,154] have been developed by using electromagnetic cavities (i.e. static electromagnetic fields). In [156,157,158] dynamical electromagnetic fields (i.e. wave guides) have been studied always for the purpose of detecting relic gravitons. Yet a different approach to the problem has been described in [159]. In [158] an interesting prototype detector was described with frequency of operation of the order of 100 MHz (see also [160]). It is not clear if, in the near future, the improvements in the terrestrial technologies will allow the detection of relic gravitons for frequencies, say, larger than the MHz. It is, however, a rather intriguing possibility.
at most as ν and since there is a ν −6 in the denominator, the main contribution to the integral should occur for ν < 0.1 kHz. This argument can be explicit verified in the case of the calculations carried on in [68] and it would be interesting to check it also in our improved framework.

Final remarks
The incoming score year might witness direct experimental evidences of relic gravitons either from small frequency experiments or from high-frequency experiments. By low-frequency experiments we mean, as in the bulk of this review, the CMB experiments possibly analyzed together with the two remaining cosmological data sets (i.e. large-scale structure determinations of the matter power spectrum and type Ia supernova observations). By high-frequency experiments we mean the appropriately advanced versions of wide band interferometers such as Ligo, Virgo, Geo and Tama.
The main observables related to relic gravitons have been reviewed in a self-contained manner and in the framework of the ΛCDM paradigm with specific attention to the complementarity between the low-frequency and the high-frequency branches of the relic graviton spectrum. Any model claiming a signal coming from high-frequency gravitons should be compatible with the ΛCDM paradigm in the low-frequency branch of the relic graviton spectrum.
It is instructive to go back to the comparison drawn, in Fig. 1, between the electromagnetic spectrum and the relic graviton spectrum. The gap between the graviton frequencies explored by CMB experiments and the graviton frequencies probed by wide-band interferometers is of the order of the frequency gap between radio waves and γ-rays.
A detection of long wavelength gravitons in CMB experiments can be direct (i.e. from the B-mode polarization) or indirect (i.e. from some global fit of CMB observables including the tensor contribution). In the context of the ΛCDM paradigm the CMB detection of long wavelength gravitons will fix the overall normalization of the spectral energy density. Even in the absence of such a direct detection, the current upper limits on the contribution of long wavelength gravitons to CMB observables implies a minute signal at higher frequencies. The (hoped) sensitivities achievable by the advanced wide-band interferometers are still insufficient for a direct detection of high-frequency gravitons.
The latter statement summarizes the standard lore of the problem which may well be realized. At the same time it might be unwise to assume (or presume) that the current success of the ΛCDM paradigm also fixes the whole thermal history of the Universe for temperatures larger than the MeV. The general ideas conveyed in the present review suggest that the high-frequency branch of the relic graviton spectrum is rather sensitive to the whole post-inflationary thermal history of the Universe. If the post-inflationary evolution is dominated by stiff sources, for instance, it is not impossible, as explicitly shown, to have positive detection of relic gravitons both at small and high frequencies even enforcing the current bounds on the tensor contribution to CMB observables.
If we assume (or strongly believe) the standard lore (i.e. that relic gravitons will probably not be seen by wide-band detectors) it is useless to demand theoretical accuracy. For instance it is useless to ask what would be the effect of changing Ω Λ on a signal which is anyway 6 or even 7 orders of magnitude smaller than the most optimistic sensitivities. A positive detection of relic graviton backgrounds at high-frequencies would demand, however, more accurate estimates of the theoretical signal in different models. Absent a direct detection of relic gravitons by wide-band inteferometers, accurate theoretical calculations can be used to set bounds on possible deviations of the post-inflationary thermal history.