Post-inflationary thermal histories and the refractive index of relic gravitons

We investigate the impact of the post-inflationary thermal histories on the cosmic graviton spectrum caused by the inflationary variation of their refractive index. Depending on the frequency band, the spectral energy distribution can be mildly red, blue or even violet. Wide portions of the parameter space lead to potentially relevant signals both in the audio range (probed by the advanced generation of terrestrial interferometers) and in the mHz band (where space-borne detectors could be operational within the incoming score year). The description of the refractive index in conformally related frames is clarified.


Introduction
Stochastic backgrounds of cosmological origin have been suggested more than forty years ago [1,2,3] as a genuine general relativistic effect in curved space-times. Since the evolution of the tensor modes of the geometry is not Weyl-invariant [1], the corresponding classical and quantum fluctuations can be amplified not only in anisotropic metric but also in conformally flat background geometries [2,3] (see also [4]). For this reason backgrounds of relic gravitons are expected, with rather different properties, in a variety of cosmological scenarios and, in particular, during an isotropic phase of quasi-de Sitter expansion [5]. The backgrounds of cosmic gravitons are analyzed in terms of the spectral energy distribution in critical units, conventionally denoted by Ω gw (ν, τ 0 ) where τ 0 is the present value of the conformal time coordinate and ν is the comoving frequency whose numerical value coincides (at τ 0 ) with the value of the physical frequency 2 . The transition from the radiation-dominated to the matter stage of expansion leads to the infrared branch of the spectrum ranging between the aHz and 100 aHz [6,7,8]. The standard prefixes shall be used throughout, i.e. 1 aHz = 10 −18 , Hz, 1 mHz = 10 −3 Hz, 1 MHz = 10 6 Hz and so on and so forth.
Between few aHz and 100 aHz the low frequency branch of the spectrum is universal and it is caused by the tensor modes of the geometry reentering after matter-radiation equality. For higher frequencies the spectral energy distribution bears the mark of the evolution of the Hubble rate prior to the radiation-dominated epoch. The simplest possibility (so far consistent with observational data) is that a quasi-de Sitter phase of expansion is followed by a radiationdominated stage: in this case the spectral energy density is quasi-flat [5,9,10,11] between 100 aHz and 100 MHz. Neglecting all possible complications (damping of the tensor modes due to neutrinos [12,13], evolution of relativistic species [14,15], late time dominance of the dark energy [15]) we can estimate 3 the typical amplitude of the spectral energy distribution in critical units which is h 2 0 Ω gw (ν, τ 0 ) ≤ O(10 −16.5 ) for frequencies ranging between the mHz and 10 kHz. This minute result follows from the absolute normalization of the spectral energy distribution fixed by the upper limit on the tensor to scalar ratio r T (ν p ) at the pivot frequency 4 ν p = k p /(2π) = 3.092 aHz.
In spite of the fact that the combination of different Cosmic Microwave Background (CMB in what follows) observations imply various sets of upper bounds on r T (ν p ), here we shall be enforcing the limit r T < 0.06 for the tensor spectral index. This limit follows from a joint analysis of Planck and BICEP2/Keck array data [16] (see also [17]). However, in a less conservative perspective we could require that r T (ν p ) < 0.17, as demanded by the WMAP9 results [18,19]. This particular figure holds if the WMAP9 data are combined with the baryon acoustic oscillation data [20], with the South Pole Telescope data [21] and with the Atacama Cosmology Telescope data [22]. The WMAP9 data (combined with further data sets) lead to bounds on r T (ν p ) that are grossly similar, for the present ends, to the Planck Explorer data suggesting r T (ν p ) < 0.11 [23].
One of the tacit (but key) assumptions of the concordance scenario is that radiation dominates almost suddenly after the end of inflation. This assumption is used, among other things, to assess the maximal number of inflationary efolds today accessible by CMB observations. It is however not unreasonable to presume that in its early stages the Universe passed through different rates of expansion deviating from the radiation-dominated evolution. The slowest possible rate of expansion occurs when the sound speed of the medium coincides with the speed of light [24] (see also [25]). Expansion rates even slower than the ones of the stiff phase can only be realized when the sound speed exceeds the speed of light. This possibility is however not compatible with the standard notion of causality. The plausible range for the existence of such a phase is between the end of inflation and the formation of the light nuclei [26,27,28]. If the dominance of radiation is to take place already by the time of formation of the baryon asymmetry, then the onset of radiation dominance increases from few MeV to the TeV range. In particular, if the post-inflationary plasma is dominated by a stiff source (i.e. characterized by a barotropic index w = p/ρ larger than 1/3) the corresponding spectral energy density inherits a blue (or even violet) slope for typical frequencies larger than the mHz and anyway smaller than 100 GHz. In this case, depending on the parameters characterizing the stiff evolution, the spectral energy distribution can be of the order of 10 −10 in the audio band while in the mHz band is at most 10 −15 .
In quintessence scenarios the present dominance of a cosmological term is translated into the late-time dominance of the potential of a scalar degree of freedom that is called quintessence (see e.g. [29]). If we also demand the existence of an early inflationary phase accounting for the existence of large-scale inhomogeneities, the inflaton potential must dominate at early times while the quintessence potential should be relevant much later (see second and third papers in Ref. [26] and [27]). In between the scalar kinetic term of inflaton/quintessence field dominates the background. When the inflaton and the quintessence field are identified the existence of this phase is explicitly realized [27] even if a similar phenomenon may take place also in slightly different scenarios.
Gravitational waves might acquire an effective index of refraction when they travel in curved space-times [30,31] and this possibility has been recently revisited by studying the parametric amplification of the tensor modes of the geometry during a quasi-de Sitter stage of expansion [32]: when the refractive index mildly increases during inflation the corresponding speed of propagation of the waves diminishes and the power spectra of the relic gravitons are then blue, i.e. tilted towards high frequencies. The purpose of this paper is to compute the spectral energy distribution of the relic gravitons produced by a dynamical refractive index without assuming a standard post-inflationary thermal history.
Even though the current upper limits on stochastic backgrounds of relic gravitons are still far from their final targets [33,34], the advanced Ligo/Virgo projects are described in [35,36].
We shall then suppose, according to Refs. [35,36] that the terrestrial interferometers (in their advanced version) will be one day able to probe chirp amplitudes O(10 −25 ) corresponding to spectral amplitudes h 2 0 Ω gw = O(10 −11 ). In the foreseeable future there should be at least one supplementary interferometer operational in the audio band namely the Japanese Kamioka Gravitational Wave Detector (for short Kagra) [48,49] which is, in some sense, the prosecution and the completion of the Tama-300 experiment [45]. In the class of wide-band detectors we should also mention the GEO-600 experiment [46] (which is now included in the Ligo/Virgo consortium [47]) and the Einstein telescope [50] whose sensitivities should definitively improve on the advanced Ligo/Virgo targets.
The target sensitivity to detect the stochastic background of inflationary origin should correspond to a chirp amplitude h c = O(10 −29 ) (or smaller) and to a spectral energy distribution in critical units h 2 0 Ω gw = O(10 −16 ) (or smaller). These orders of magnitude estimates directly come from the amplitude of the quasi-flat plateau produced in the context of single-field inflationary scenarios; in this case the plateau encompasses the mHz and the audio bands with basically the same amplitude. Even though these sensitivities are beyond reach for the current interferometers, a number of ambitious projects will be supposedly operational in the future. The space-borne interferometers, such as (e)Lisa (Laser Interferometer Space Antenna) [37], Bbo (Big Bang Observer) [38], and Decigo (Deci-hertz Interferometer Gravitational Wave Observatory) [39,40], might operate between few mHz and the Hz hopefully within the following score year. While the sensitivities of these instruments are still very hypothetical, we can suppose (with a certain dose of optimism) that they could even range between h 2 0 Ω gw = O(10 −12 ) and h 2 0 Ω gw = O(10 −15 ). The layout of the paper is the following. In section 2 the basic action of the problem shall be analyzed in its different forms. For the sake of completeness the relation between different parametrizations will also be discussed with the purpose of arguing that the physical description does not change. In section 3 we shall analyze the evolution of the effective horizon and the the amplification of the relic gravitons. Section 4 will be focussed on the analytic (though approximate) estimates of the graviton spectra while section 5 contains the discussion of the detectability prospects. Some concluding remarks are collected in section 6.

Gauge-invariance and frame-invariance
Gravitational waves might acquire a refractive index when they evolve in curved space-times [30,31] and the impact of this idea on a quasi-de Sitter stage of expansion has been explored in [32] where the presence of a (time dependent) refractive index has been introduced for the first time. For standard dispersion relations the propagating speed of the tensor modes of the geometry in natural units coincides with the inverse of the refractive index (i.e. c gw (τ ) = 1/n(τ )) and the basic action can be written, in a covariant language, as: where P = √ 8πG = 1/M P and P µν (E) is the spatial projector tensor orthogonal to u (E) µ : In the case n → 1 the action of Eq. (2.1) reproduces the original Ford and Paker result [3]; in comoving coordinates 5 and in a conformally flat metric Eq. (2.1) assumes the following form: 3) The inverse of the refractive index multiplies each spatial derivative of the tensor amplitude (see Eq. (2.3)). This is the parametrization employed in Refs. [30,31,32] and it is physically motivated. It is however possible to adopt a somehow contrived viewpoint and to describe the dynamics of the refractive index with an apparently different action namely: To get from Eq. (2.3) to Eq. (2.4) we need a specific transformation that leaves unaltered the conformal time coordinate and the tensor amplitude while the scale factor is simply rescaled through the refractive index itself: This transformation exist and it is nothing but a conformal rescaling. Indeed Eq. (2.5) transforms separately the background and the tensor inhomogeneities but it is not difficult to see that it comes directly from a conformal rescaling that leaves unaltered the tensor fluctuations of the geometry. To make this point more apparent, let us verify explicitly that the transformation (2.5) is just a particular case of the following conformal rescaling of the four-dimensional metric: that implies the transformation (2.5) on the background and on the related inhomogeneities. From Eq. (2.6) the transformation for the background is immediate and it is given by g (E) µν = Ω 2 G µν . In the case of a conformally flat metric of Friedmann-Robertson-Walker type we have g (E) µν = a 2 E (τ ) η µν and this means G (E) µν = a 2 (τ ) η µν . Thus, since a E = Ω a this transformation coincides exactly with the first relation of Eq. (2.5) provided, as anticipated, Ω ≡ n.
The second relation reported in Eq. (2.5) is also a general consequence of the conformal rescaling (2.6). While the proof of this statement is immediate in the case of the tensor modes, it is useful to present the complete argument since there are also symmetric implications in the case of the scalar modes. Neglecting, for simplicity, the vector modes of the geometry the fluctuations of the metric in the Einstein frame are µν ( x, τ ), (2.7) 5 In the case of a conformally flat metric g where δ t g (E) µν and δ s g (E) µν denote respectively the tensor and the scalar fluctuations of the geometry in the Einstein frame. In a conformally flat bacgkground geometry of Friedmann-Robertson-Walker type Eq. (2.7) becomes ij is the (divergenceless and traceless) tensor amplitude appearing in Eq. (2.1). By definition the tensor amplitude h Let us now consider exactly the same decomposition in the conformally related frame defined by Eq. (2.6); as in the case of Eq. (2.7), the G µν can be decomposed into a homogeneous part supplemented by its own tensor inhomogeneities: (2.10) where this time the explicit form of the tensor and scalar fluctuations of the four-dimensional metric will be given by: The tensor amplitude h ij defined in Eq. (2.11) is gauge-invariant while the scalar fluctuations of Eq. (2.12) are not immediately gauge-invariant. To work out the relation between the fluctuations in the two frames we can therefore start with the tensor modes; from Eq. (2.6), recalling the explicit forms of the tensor fluctuations in the two frames (i.e. Eqs. (2.8) and (2.11)) we can write where the projectors and of the four-velocities have been conformally rescaled as If we now posit that the conformal factor with the refractive index itself coincide (i.e. Ω(τ ) = n(τ )), the action of Eq. (2.16) becomes exactly, Let us finally mention, as we close the section, the analog results for the scalar modes of the geometry which are however less central to the discussion of the present investigation. Indeed from Eq. (2.6) we will have that Recalling then the explicit results of Eqs. (2.9) and (2.12), Eq. (2.18) implies a specific relation between the perturbed components in the two frames, i.e.
Equation (2.19) implies that, unlike their tensor counterparts, the scalar inhomogeneities defined in Eqs. (2.9) and (2.12) are neither gauge-invariant nor frame-invariant. Note, however, that the curvature perturbations on comoving orthogonal hypersurfaces are both gaugeinvariant and frame-invariant. Indeed in the two frames they are simply given by Equation (2.20) does not imply that R E = R, as it could be superficially concluded. On the contrary, the mismatch between ψ E and ψ is exactly compensated by the mismatch between H E and H. In fact, from the relation between the background scale factors (i.e. a E = a n) we have 2(H E − H) = n /n so that Eq. (2.20) implies R = R E . We therefore have, as anticipated, that the tensor modes of the geometry discussed in the bulk of the paper and the curvature perturbations on comoving orthogonal hypersurfaces are both frame-invariant and gauge-invariant 6 . We finally mention that after the appearance of Ref. [32], two similar papers [41] pursued the same idea. The two approaches ultimately coincide since they are related by a conformal rescaling involving the refractive index. More specifically, to get from the description of Ref. [32] to the one of Ref. [41] it is sufficient to make a conformal rescaling and to parametrize the propagating speed or the refractive index as a power of the scale factor. Following the suggestion of Ref. [32], the authors of Ref. [41] considered the evolution of the refractive index in an inflating background. This choice is however potentially confusing: since the two descriptions are related by a conformal rescaling the two backgrounds should also be conformally related [42]. This would mean, in practice, that if the background inflates in the Einstein frame, it might not inflate in the conformally related frame. However, since the choice of the pivotal frame where the background inflates is not constrained, the choice of Ref. [41] is, in a sense, mathematically legitimate but physically superficial especially in the light of the previous literature. We are therefore in the situation where the two conformally related actions are simply two complementary parametrizations of the same effect. To cope with this unwanted ambiguity the easiest solution is to define a generalized action for the tensor modes encompassing the various possibilities suggested so far. As we shall see in sections 4 and 5 when γ = 0 (and, in particular, when γ = 1) the spectral index determined in the γ = 0 case is just rescaled by a γ-dependent prefactor that can be reabsorbed in a redefinition of the spectral index.

Effective horizons 3.1 Generalities
According to the results obtained so far the action describing the evolution of the tensor modes of the geometry in the presence of a dynamical refractive index can be parametrized in the following manner: With the redefinition (3.2) of the time coordinate, Eq. (3.1) can be rewritten as The function b(η) plays the role of an effective scale factor: note, in fact, that in the limit n → 1 we have that η coincides with τ and that, consequently, b(η) → a(τ ). When n = 1 the evolution of b(η) defines an effective horizon, namely: where the prime denotes a derivation with respect to the conformal time coordinate τ while the overdot denotes a derivation with respect to the η-time (and not a derivation with respect to the cosmic time coordinate, as in the conventional notations). To clarify this point and to avoid potential confusions the following relations are explicitly given: which can be verified by using Eq. (3.2) and the relation of τ the cosmic time coordinate t, i.e.

The canonical Hamiltonian and the mode functions
In terms of the canonical normal modes Up to a total time derivative Eq. (3.8) can also be written as: Since µ ij ( x, η) is given as the sum over the two polarizations ⊕ and ⊗ the action (3.9) becomes immediately: From Eqs. (3.11) and (3.12) the canonical momenta are π λ =μ λ ; consequently the canonical Hamiltonian associated with Eqs. (3.11) and (3.12) is given by: The commutation relations at equal η-times together with the explicit form of the Hamiltonian (3.13) lead directly to the evolution equations of the operatorsμ λ andπ λ : The Fourier representations ofμ λ andπ λ is: The evolution of the mode functions f k,λ and g k,λ follows from Eq. (3.15) while the normalization of their Wronskian is a consequence of the commutation relations of Eq. (3.14): The equations for the mode functions reported in Eq. (3.18) can be decoupled as: where the polarization index has been omitted since the result of Eq. (3.20) holds both for ⊕ and for ⊗. By recalling thatĥ ij b =μ ij the mode expansion of the tensor amplitudeĥ ij ( x, η) in the η-time is given by: where the explicit form of the two polarizations can be written as: andk i = k i /| k|,m i = m i /| m| andq = q i /| q| are three mutually orthogonal directions andk.
If we now represent the field operatorĥ ij ( x, η) in Fourier space: we also have from Eqs. (3.21) and (3.23): It follows from Eq. (3.24) that the two-point function in real and in Fourier space is given by The two-point functions computed from Eq. (3.24) are simply 7 where j 0 (kr) is the spherical Bessel function of zeroth order [43,44]. The tensor power spectrum of Eqs. (3.25) and (3.26) is then given by

Evolution of the effective horizon
The variation of the refractive index can be measured in units of the Hubble rate in full analogy with what it is customarily done in the case of the slow-roll parameter, namely: Equation (3.29) also implies that the evolution of n(a) could be considered as piecewise continuous across a certain critical value of the scale factor a * ; more specifically the situation we are interested in is the one where while n(a) = 1 for a > a * . It is relatively simple to imagine a number of continuous interpolation between the two regimes but what matters for the present considerations is overall the continuity of n(a), not the specific form of the profile across the normalcy transition. One of the simplest possibilities is given by 8 n(a, ξ) = n i (a/a i ) α e −ξa/a * + 1, going as a α for a < a * and approaching 1 quite rapidly when a > a * and ξ > 1. The typical scale a * (roughly corresponding to the maximum of n(a)) may coincide with the end of the inflationary phase but this possibility is neither generic nor compulsory. The value of a * corresponds to a critical number of efolds N * which is of the order of N t (i.e. the total number of efolds) if a * marks the end of the inflationary phase. This identification is however not mandatory and it will also be relevant, from the physical viewpoint, to consider the case N * < N t or even N * N t . It is relevant to mention, for future convenience, thatḃ ≥ 0; this means that b(η) is always an increasing function of the η coordinate defined in Eq. (3.2). This observation is important for the forthcoming estimates of the cosmic graviton spectrum (see section 4 and discussions therein). More specifically, if we consider separately the cases γ = 0 and γ = 1 Eq. (3.2) implies: where n * = n i (a * /a i ) α . The functions of Eq. (3.31) and (3.32) are always increasing 9 as a function of x = a/a * . Since we shall bound the attention to the case of expanding scale factors, Eqs. (3.31)-(3.32) imply that the explicit evolution of b (either in τ or in η) is always monotonically increasing. The relations between η, τ and the Hubble radius during the refractive phase are affected by the value of the slow-roll parameter. This is a generic consequence of Eq. (3.7) that fixes the relation between η and the conformal time coordinate: If we now integrate Eq. (3.33) by parts we will have: implying, together with Eq. (3.33), that The pump fieldb/b of Eq. (3.20) during the refractive phase can be written as: which can also be written as: The same result can be obtained by assuming a slow-roll phase b = a n γ−1 where a(τ ) = (−τ /τ * ) −β and β = 1/(1 − ). Thus thanks to Eq. (3.7) we have The result of Eq. (3.40) implies: If we now computeb/b from Eq.

Cosmic graviton spectra and thermal histories
The cosmic graviton spectra can be estimated analytically [11,26,14,52,53] by adapting some of the standard methods 10 and by noting that Eq. (3.20) is equivalent to an integral equation whose initial conditions are assigned at the reference time η ex : where η ex is defined as the turning point at which the solution to Eq. (3.20) changes its analytic form: Equation (4.2) can be dubbed by saying that at η ex the given mode k exits the effective horizon defined by the evolution of b; the second requirement of Eq. (4.2) is for the moment pleonastic since the exit always occurs in a regime whereb ex = 0. Even though b(η) never evolves linearly in the vicinity of the exit, this occurrence may arise close to the reentry that defines the second relevant turning point of the problem.

The large-scale power spectra
Neglecting the terms O(k 2 η 2 ) the lowest order solution of Eq. (4.1) is: where, according to Eq. (3.20),ḟ k (η ex ) = g k (η ex ) andḟ k (η) = g k (η). Equations (4.3) and (4.4) determine the approximate form of the power spectrum for wavelengths larger than the 10 These methods must be revisited in a slightly different perspective since the evolution of η and of the conformal time coordinate only coincide, in the present framework, after the end of inflation.
Hubble radius. Since the second term appearing inside the squared bracket at the right hand side of Eq. (4.3) is subleading for typical wavelengths larger than the effective horizon, after inserting Eq. (4.3) into Eq. (3.27) the tensor power spectrum becomes: where Eq. (3.41) has been used to get an explicit expression of b(η) in the regime η < −η * . The amplitude |A| appearing in Eq. (4.5) parametrizes, up to an irrelevant phase, the mismatch between the exact and the approximate solutions at η ex : for k 2 |b/b| the correctly normalized solutions of Eq. (3.20) are f k (η) = e ± ikη / √ 2k. However as soon as η ex is approached the amplitude gets slightly modified and by recalling Eq. (3.38) the exact solution of Eq. (3.20) can be written in terms of Hankel functions [43,44] where H (1) µ (kη) is the Hankel function of the first kind 11 . For wavelengths larger than the Hubble radius the Hankel function of Eq. (4.6) can be expanded in the limit |kη| 1 so that thanks to Eq. (3.27) the tensor power spectrum becomes Since 3 − 2µ = 2(1 − ν) (as implied by Eq. (4.6)), the ratio between Eqs. (4.5) and (4.7) implies that: where ν has been defined in Eq. (3.41). The value of |A(µ)| estimates the theoretical error of the treatment based on Eq. (4.5) and on the approximate form of the mode functions.
While it is often plausible to neglect the complication 12 of A(µ) and simply set A(µ) → 1, at low frequencies the absolute normalization of the cosmic graviton spectrum is however very sensitive to the value of the mode functions for η = O(η ex ). It is then mandatory to use Eq. (4.7) which can also be expressed as: where H * denotes the Hubble rate at η * . Equation (4.9) is the large-scale power spectrum valid for k < a * H * . The scales that exited the Hubble radius for η > −η * have a different spectral slope and, in this respect, we have a twofold possibility. If η * coincides with the end of inflation, then the power spectrum will still be given by Eq. (4.12) where, however, N t = N * . Conversely if the refractive phase terminates before the end of inflation the power spectrum will have a further branch for a * H * < k ≤ a 1 H 1 : where n T = −2 /(1 − ). It is relevant to remark that in the limit α → 0 we have µ → where µ is the Bessel index appearing in Eq. (4.6). Equation (4.10) can be further modified by appreciating that since between −η * and −τ 1 the background inflates we have H * a * = (H 1 a 1 )e N * −Nt , n * = n i (a * /a i ) α ≡ n i e αN * (4.11) Taking into account Eqs. (4.10) and (4.11) the power spectrum (4.9) finally becomes: where M P = √ 8π M P (see also the definitions after Eq. (2.1)). Equation (4.12) determines the tensor to scalar ratio whose explicit form is (4.14) where we defined, for the sake of conciseness, k max = a 1 H 1 .

The power spectra after reentry
Terrestrial interferometers and space-borne detectors operate at the present time and will necessarily measure the cosmic graviton spectrum for typical wavelengths shorter than the Hubble radius. While the largest wavelengths of the problem (i.e. smallest k-modes) reentered after matter-radiation equality the shortest wavelengths (i.e. largest k-modes) crossed the effective horizon at different epochs after the end of the inflationary stage of expansion. The reentry depends on the post-inflationary thermal history and on the expansion rate that can be very different from the one of a radiation-dominated plasma. When the refractive index is not dynamical the previous observation leads to a characteristic class of violet spectral energy distribution [26,28,53] and it will be interesting to see what happens in the present situation.
Provided the reentry occurs when: then kη re = O(1). However, as already remarked above (see Eq. (4.2)), ifb re → 0 in the vicinity of the turning point, then kη re 1. For η ≥ η re the solution of Eq. (3.20) can be expressed as: where f re (η) and g re (η) are the mode functions inside the effective horizon (i.e. quantum mechanically normalized plane waves in the crudest approximation). From the continuity of f k (η) and g k (η), Eqs. (4.3) and (4.4) imply: By continuity Eq. (4.16) evaluated at η re must coincide with Eqs. (4.17)-(4.18); thus c ± (k) can be determined after some simple algebraic manipulation 13 : Inside the Hubble radius the mode functions are plane waves or, more precisely, the limit of Hankel functions for large values of their arguments [43,44]. We can then express directly Eqs. (4.20) and (4.21) bu using the plane wave limit of the correspoding mode functions: Since the coefficients c ± (k) satisfy |c + (k)| 2 − |c − | 2 = 1 it is sufficient to determine just one of the two square moduli. If the exit occurs for η < −η * and the reentry takes place when the refractive index is not dynamical, Eqs. (4.22) and (4.23) can be written more explicitly where this time .
Equation (4.27) allows for a swift determination of the power spectrum and of the spectral energy distribution in the limit kτ 1, i.e. when the relevant wavelengths are all inside the Hubble radius: By taking the ratio between Eqs. (4.28) and (4.29) we recover the standard relation between the power spectrum and the spectral energy density valid when the relevant wavelengths are shorter than the Hubble radius at a given epoch: (4.30) Consequently inside the Hubble radius we can evaluate indifferently either the power spectrum or the spectral energy distribution.

Different thermal histories
The different thermal histories and their effects on the spectral energy distribution can be understood by drawing the salient features of the effective horizon in various physical situations.
In Fig. 1 on the vertical axis we plot the common logarithm of F =ḃ/b and we simultaneously compare it with the wavenumbers of the problem 14 . In practice the conditions kη ex = O(1) and kτ re = O(1) will always be verified except that in the case of a reentry during radiation when kτ re 1 and a re = 0. According to Fig. 1 we have three different classes of modes: i) the modes exiting the effective horizon during the refractive phase and reentering after equality (i.e. k < a eq H eq ); ii) the modes exiting the effective horizon during the refractive phase and reentering during radiation (i.e. a eq H eq < k < a * H * ); iii) the modes exiting the effective horizon after the end of the refractive phase and reentering during radiation (i.e. a * H * < k < a 1 H 1 ). The different regions are separated in Fig. 1 by three horizontal arrows and the typical wavenumbers (i. e. k 1 , k * and k eq ) define the three branches of the spectral energy density (or of the power spectrum). Either the inflationary phase continues after b * or the radiation-dominated epoch suddenly kicks in. Between these two possibilities the former is more generic than the latter which would correspond, in the notation of Fig. 1, to the limit b * → a 1 ; this is why, in Fig. 1, we preferred to distinguish clearly the two scales by assuming b * a 1 . The three different branches of the spectral energy distribution illustrated in Fig. 1 can be deduced from Eqs. (4.27), (4.28) and (4.29); the result of this computation is: a eq H eq < k ≤ a * H * , (4.32) k ≤ a eq H eq , (4.33) where B(b * , n * , , n T ) is given by 15 (4.34) The spectral index n T appearing in Eqs. (4.31)-(4.34) is instead: where the second equality follows in the limit 1. As it must, the exact expression of Eq. (4.35) coincides with Eq. (4.12). The quasi-flat branch of Eq. (4.31) is caused by the modes that exited the effective horizon for a > b * and reentered during the radiation-dominated epoch (i.e. for a > a 1 ). The second branch of the spectrum, reported in Eq. (4.32) involves the modes that exited the effective horizon during the refractive phase and reentered all along the radiation stage. Finally the standard infrared branch corresponds to modes that exiting the effective horizon during the refractive epoch and reentering during the matter-dominated phase (i.e. k < k eq in the terminology of Fig. 1). There are no compelling reasons why the physical situation illustrated by Fig. 1 should be considered preferable to some other potentially viable evolution of the effective horizon. Different thermal histories can be envisaged and cannot be ruled out by present version of the concordance scenario. Prior to nucleosynthesis, there are no direct tests of the thermodynamical state of the Universe and, therefore, the effective equation of state of the primeval plasma can be arbitrarily different than the one of radiation. In Fig.  2 the effective horizon is illustrated in the case where the post-inflationary expansion rate is slower than the one of a radiation-dominated plasma; in this case between a 1 and a s we havė where w denotes the barotropic index of the stiff post-inflationary phase. Note, for comparison, that H ∝ a during inflation while, in the radiation stage, H ∝ 1/a (see aslo Fig. 1). In the case of a stiff post-inflationary phase we have instead that (at most) H ∝ a −2 , as implied by Eq. (4.36) for w → 1. The evolution sketched in Fig. 2 leads to a spectral energy distribution characterized by four different branches: the wavenumbers k < a eq H eq correspond to scales hitting the effective horizon the first time during the refractive phase and reentering after matter radiation equality: their spectral energy distribution will then have the same slope of Eq. (4.33). Following the same way of reasoning Ω gw ∝ |kτ * | n T whenever a eq H eq < k < a * H * : in Fig. 2 this part of the spectrum corresponds to those modes exiting during the refractive phase and reentering during the radiation epoch. The supplementary branch of the spectrum implied by Fig. 2 is caused by those modes exiting in the course of the inflationary phase and reentering during the stiff phase: in this branch the spectral energy density scales as Ω gw ∝ |kτ s | m T where the spectral index m T is now given by: If the expansion rate is slower than radiation the slope in this branch can be very steep (i.e. even violet) with m T = O(1) in the limit w → 1. Incidentally if the expansion rate is faster than radiation 16 it can happen that m T < 0. While the cases illustrated by Figs. 1 and 2 are the most promising from the viewpoint of the potential signals (as we shall see in the following section), there are other possible evolutions of the effective horizon where the resulting spectral energy distribution does not have a flat (or decreasing) plateau and it always increases. In this connection Fig. 3 illustrates a possibility complementary to the one of Fig. 2 but leading to a different spectrum. Both in Figs. 2 and 3 a stiff phase precedes the ordinary radiation epoch. However the refractive and the stiff phases of Fig. 3 are longer then in Fig. 2. This occurrence implies the possibility of modes exiting the effective horizon during the refractive phase and reentering during the stiff phase. This different dynamical situation implies that the intermediate branch of the spectrum (i.e. a * H * < k < a s H s ) is not quasi-flat anymore. Furthermore using Eqs. (4.27) and (4.29), Ω gw scales as |kτ s | s T where now s T is given by increasing (thought at a different rate) all the energy of this spectrum will be concentrated in the highest frequency regime and this is the reason why the detectability prospects are, in this situation, less promising than in the case of a sufficiently long plateau at high frequency.
Various other examples can be analyzed by using the approximate methods described in this section but they are not central to the present discussion.

Basic considerations
The phenomenological signatures of the relic gravitons are customarily assessed by using the comoving frequency that is defined as ν = k/(2π) where k denotes the comoving wavenumber. Four complementary quantities can be used to describe the cosmic graviton background: i) the tensor power spectrum (denoted by P T (ν, τ 0 )), ii) the spectral energy distribution (i.e. Ω gw (ν, τ 0 )), iii) the chirp amplitude h c (ν, τ 0 ) and iv) the spectral amplitude S h (ν, τ 0 ) (measured in units of Hz −1 = sec). Except for S h (ν, τ 0 ) the three remaining variables are dimensionless. The chirp amplitude is, by definition, h 2 c (ν, τ 0 ) = P T (ν, τ 0 )/2. The tensor power spectrum at the present time can be related to the spectral energy distribution as  τ 0 ) and since the detectors of gravitational radiation are operating in the audio band (i.e. between few Hz and 10 kHz) it is useful to stress the explicit relations between the various quantities mentioned above for typical frequencies ν = O(100) Hz where the sensitivities of wide-band detectors to cosmic graviton background are (approximately) maximal 17 :

Pivotal frequencies
The spectral energy distribution is characterized by various typical frequencies which are determined from the wavenumbers appearing in Figs. 1, 2 and 3. The smallest frequency range of the spectrum follows from the pivot wavenumber k p at which the scalar and tensor power spectra are assigned [16,17,18,19]: The frequency associated with the dominance of dark energy is of the same order of Eq. (5.6) and it is fixed by Ω M 0 and Ω Λ ; in the case of the concordance paradigm we have Since the equality wavenumber is k eq = 0.0732 [h 2 0 Ω R0 /(4.15 × 10 −5 )] −1/2 h 2 0 Ω M 0 Mpc −1 the related frequency ν eq is: Hz. (5.8) The frequency ν bbn = O(10 −2 ) nHz enters directly the big-bang nucleosynthesis constraint (see below Eq. (5.13)) and sets the scale for the suppression of the cosmic graviton background due to neutrino free-streaming [12,13]. The explicit expression of the big-bang nucleosynthesis frequency is 18 : ν bbn = 2.252 × 10 −11 g ρ 10.75 Hz. (5.9) The presence of a refractive phase illustrated in Figs. 1 and 2 introduces two further frequencies: Hz, (5.11) where A R denotes the amplitude of the power spectrum of curvature inhomogeneities at the wavenumber k p . Even if Eq. (5.11) suggests that ν max = O(200) MHz, the value of the endpoint frequency of the spectrum may exceed ν max since it depends on the post-inflationary thermal history [32]. For the thermal histories of Figs. 2 and 3 the spectral energy distribution may extend up to ν spike = ν max /σ > ν max (with σ < 1). While a similar spike (with different physical features) may also appear when the refractive index is not dynamical [32], in this particular case two new frequencies appear and they are defined as where H r denotes the Hubble rate at the onset of the radiation dominance, i.e. right after the stiff phase. The difference between ν max and ν spike comes essentially from the redshift during the stiff stage of expansion.

Phenomenological constraints
In the low-frequency range the tensor to scalar ratio of Eq. (4.14) is bounded from above not to conflict with the observed temperature and polarization anisotropies of the CMB; in the present analysis we specifically required r T (ν p ) < 0.06, as it follows from a joint analysis of Planck and BICEP2/Keck array data [16]. As already mentioned in the introduction slightly less restrictive bounds are often used in the current literature and they amount to demanding r T (ν p ) < O(0.1) [18,19,23]. The pulsar timing measurements impose instead the limit Ω gw (ν pulsar , τ 0 ) < 1.9 × 10 −8 at the frequency ν pulsar = O(10) nHz corresponding to the inverse of the observation time along which the pulsars timing has been monitored [54,55,56,57,58,59]. The big-bang nucleosynthesis sets an indirect constraint on the extra-relativistic species (and, among others, on the relic gravitons) at the time when light nuclei have been formed [60,61,62]. This limit is often expressed in terms of ∆N ν representing the contribution of supplementary (massless) neutrino species (see e.g. [63]) but the extra-relativistic species do not need to be fermionic. If, as in our case, the additional species are relic gravitons we will have to demand that: The bounds on ∆N ν range from ∆N ν ≤ 0.2 to ∆N ν ≤ 1 so that the right hand side of Eq. (5.13) turns out to be between 10 −6 and 10 −5 . The basic considerations discussed here can be complemented by other bounds which are, however, less constraining than the ones mentioned above. The same logic employed for the derivation of Eq. (5.13) can be applied at the decoupling of matter and radiation. While the typical frequency of BBN is O(10 −10 ) Hz the typical frequencies of matter radiation equality is O(10 −16 ) Hz (see Eqs. (5.8) and (5.9)). Since the decoupling between matter and radiation occurs after equality we have that (5.14) While the bound itself is numerically similar to the one of Eq. (5.13) the lower extremum of integration is smaller since ν dec ν bbn (see Eqs. (5.8) and (5.9)). The bound (5.14) (discussed in Ref. [64] with slightly different notations) has been also taken into account in the present analysis. However, since we are dealing here with growing spectral energy distributions, Eq. (5.14) is less constraining: for the same (increasing) slope the lower extremum of integration of Eq. (5.14) gives a smaller contribution than the one of Eq. (5.13).

Spectral energy distribution
The analytic estimates of the spectral energy density of section 4 lead to approximate expressions of the spectral energy distribution; however for a more quantitative assessment the cosmic graviton spectrum should be expressed in terms of T eq (ν, ν eq ), T * (ν, ν * ) and T s (ν, ν s ) denoting, respectively, the transfer functions of the energy density at low, intermediate and high frequencies: T eq (ν, ν eq ) = 1 + c eq ν eq ν + b eq ν eq ν 2 , c eq = 0.5238, b eq = 0.3537, (5.15) where the subscripts refer to the typical frequencies involved in each transition, i.e. ν eq , ν * and ν s . To transfer the spectral energy density inside the Hubble radius the procedure is to integrate numerically the equations of the tensor modes; the derivation of T eq (ν, ν eq ) and T s (ν, ν s ), in a different physical situation, has been discussed in detail in [53,66,67]. In the literature it is also customary to introduce the transfer function of the power spectrum [68,69] and the two transfer functions have slightly different numerical features that have been discussed in the past (see e.g. Ref. [53] for a comparison). With these specifications, we have: T (ν, ν eq , ν * , ν s ) = T eq (ν, ν eq ) T * (ν, ν * ) T s (ν, ν s ), (5.19) where n T = [α(3 − 2γ) − 2 ]/(1 + α − ) (see also Eq. (4.35)) and r T (ν p ) is the tensor to scalar ratio of Eq. (4.14) evaluated at the pivot frequency ν p . In the conventional case r T (ν p ) is related to the slow-roll parameter and to the tensor spectral index n T via the so-called consistency relations (see also, in this respect, Ref. [65] where a similar model for the violation of the consistency relations has been discussed). In the present situation r T (ν) and n T do not obey the consistency relations and depend on the rate of variation of the refractive index α and on the critical number of efolds N * . If the refractive index is not dynamical (i.e. α → 0 and γ → 0) we have, as expected, that n T → −2 . In Eq. (5.18) β is a parameters O(1) which depends upon the width of the transition between the inflationary phase and the subsequent radiation dominated phase; for different widths of the post-inflationary transition we can estimate 0.5 ≤ β ≤ 6.3 [66,67]. The numerical coefficients appearing in Eqs.

The constrained parameter space
We shall be predominantly interested in the possibility of a relatively strong signal in the audio and in the mHz bands. We remind that the audio band ranges between few Hz and 10 kHz where the terrestrial wide-band interferometers operate. The mHz band ranges instead between a fraction of the mHz and the Hz; in this range space-borne detectors might one day hopefully within the following score year. The MHz band, extending between 100 kHz and few GHz; this band is immaterial for potential signals coming from conventional inflationary models but could play a relevant role in the present context since, in this case, most of the signal is concentrated exactly in this region. To account for the possibility of a detection in the audio band we shall then impose on the parameter space a further constraint on the chirp amplitude: h c (ν audio , τ 0 ) > 10 −25 , ν audio = 0.1 kHz, (5.21) 19 According to Eq. (5.15), T eq (ν) → 1 for ν ν eq but the realistic situations further suppressions are expected. The neutrino free-streaming produces an effective anisotropic stress leading ultimately to an integrodifferential equation (see, for instance, [12,13]). This aspect will be discussed later when assessing the other minor sources of damping.  Fig. 4 while Fig. 5 concerns the case γ = 1. Different values of γ rescale, in practice, the values of α, as expected from the general relation connecting n T to α and γ. Indeed, to leading order in , the value of n T is the case γ → 0 is roughly thrice its value in the γ → 1 case. We therefore have that the range of α in the two situations must be rescaled by a factor of 3 and this is what we clearly see by comparing the horizontal axes in Figs. 4 and 5.
In the case of a different post-inflationary history the constrained parameter space gets modified and the relevant exclusion plots are illustrated in Fig. 6 for a fiducial choice of the parameters. In the two plots at the right the barotropic index corresponds to 2/3 while in the two plots at the left the barotropic index is maximal (i.e. w = 1). By looking at the inner and at the outer exclusion regions we conclude that a reduction in the sensitivities of the hypothetical detectors drastically reduces the areas of the parameter space. Indeed, as in Figs. 4 and 5 the inner and the outer plots correspond, respectively, to the requirements of Eqs. (5.25)-(5.26) and to the requirements of Eqs. (5.23)-(5.24). The reason for this reduction is a direct consequence of the violet spectral slope in the highest frequency domain. Still, for a given value of σ, the constrained parameter space suggests a potentially interesting signal.
The explicit profiles of different models will now be illustrated. While the same analysis can be easily rephrased either in term of the power spectrum P T (ν, τ 0 ) or in terms of the spectral amplitude S h (ν, τ 0 ) we shall be mainly interested in the chirp amplitude and in the spectral energy distribution. This kind of approach is also useful for the explicit derivation of a template family.
In Fig. 7 we illustrate the chirp amplitude and the spectral energy distribution in the case of a standard post-inflationary history. The scales of the two plots on the horizontal axis are different 20 . The explicit differences among the various models are less pronounced if we look at the chirp amplitude which is is proportional to the square root of the spectral energy distribution and it is further suppressed by one power of the frequency. It is finally relevant to stress that Figs. 7, 8 and 9 have been all derived in the case γ = 0. This is not a limitation, as we saw from the discussion of the constrained parameter space: different values of γ simply log 8/Hz shift the allowed range of α. As N * increases the plateau of Fig. 7 becomes less evident. The results of Fig. 7 can be usefully compared with the ones illustrated in Fig. 8 where the post-inflationary evolution is characterized by a stiff epoch. While in Fig. 8 we considered the case w = 1, in Fig. 9 we took instead w = 2/3. Recalling Eqs. (4.37)-(4.38) and (5.12) the value of the barotropic index controls not only the slope of the cosmic graviton spectrum in the vicinity high-frequency spike but also the frequency range. Larger values of w correspond to more violet slopes in the MHz band while the values of σ are inversely proportional to the frequency of the spike so that as σ gets smaller than 1 the position of the spike exceeds the GHz. Both σ and w determine the length of the stiff phase. In the right plots of Figs. 7, 8 and 9 we can appreciate a minor suppression for frequencies of the order of 0.01 nHz. This suppression is less visible in the chirp amplitude but it is evident from the spectral energy distribution. For ν < ν bbn the slight break in the spectrum is due to the neutrino free streaming. The neutrinos free stream, after their decoupling, and the effective energy-momentum tensor acquires, to first-order in the amplitude of the plasma fluctuations, an anisotropic stress. 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 log 8/Hz parametrized by the function F(R ν ) = 1 − 0.539R ν + 0.134R 2 ν , R ν = r r + 1 , r = 0.681 N ν 3 , (5.27) where, as usual, R ν is the fraction of neutrinos in the radiation plasma; clearly in the concordance model R γ + R ν = 1. In the case R ν = 0 (i.e. in the absence of collisionless particles) 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 due to neutrino free streaming is thus effective for frequencies larger than ν eq and smaller than ν bbn . Besides neutrino free streaming there are other two minor effects taken into account in Figs. 7, 8 and 9: the damping effect associated with the (present) dominance of the dark energy and the suppression due to the variation of the effective number of relativistic species. In the concordance scenario the redshift of Λ-dominance (i.e. (Ω Λ /Ω M 0 ) 1/3 ) determines the numerical value of ν Λ defined in Eq. (5.7). The adiabatic damping of the mode function due to the dominance of the dark energy implies a damping of the order of (Ω M 0 /Ω Λ ) 2 in the spectral energy distribution. This suppression competes with a potential increase of the spectral energy distribution for ν < ν λ and going as (ν/ν Λ ) −2 [14,70]. Finally, for temperatures much larger log 8/Hz 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 and 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 written as [66,67] (g ρ (T k )/g ρ0 )(g s (T k )/g s0 ) −4/3 where, at the present time, g ρ0 = 3.36 and g s0 = 3.90. So, in the case of the minimal standard model 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.

Concluding remarks
The evolution of the refractive index during the early stages of a conventional inflationary phase leads to a spectral energy distribution naturally tilted towards frequencies and different postinflationary thermal histories typically add a further branch to the cosmic graviton spectrum. Overall the spectrum may then exhibit up to four different branches extending between the aHz region and the GHz band. Assuming, in a minimalistic perspective, that the evolution of the refractive index terminates before the end of inflation, the spectral energy distribution involves a quasi-flat plateau at high frequencies that is supplemented by a spike between the MHz and the GHz. Depending on the thermal history the slopes of the spectral energy distribution are red, blue and even violet. After imposing the usual phenomenological constraints, there are still wide portions of the parameter space where the resulting signal could be detectable, at least in principle, either by terrestrial interferometers (in their advanced and enhanced configuration) or by space-borne detectors. In spite of less mundane possibilities leading to growing spectral energy distributions of relic gravitons, the present findings demonstrate that blue and violet spectra are compatible with conventional inflationary scenarios in the presence of a dynamical refractive index.