The propagating speed of relic gravitational waves and their refractive index during inflation

If the refractive index of the tensor modes increases during a conventional inflationary stage of expansion the relic graviton spectrum is tilted towards high frequencies. Two apparently diverse parametrizations of this effect are shown to be related by a rescaling of the four-dimensional metric through a conformal factor that involves the refractive index itself. Non-monotonic spectra with a maximum in the MHz region correspond to a limited variation of the refractive index terminating well before the end of inflation. After exploring a general approach encompassing the ones proposed so far, we estimate the required sensitivity for the direct detection of the predicted gravitational radiation and demonstrate that the allowed regions of the parameter space are within reach for some of the planned detectors operating either in the audio band (like Ligo/Virgo and Kagra) or in the mHz band (like Lisa, Bbo and Decigo).


Introduction
The tensor modes of the geometry are known to be parametrically amplified in the early Universe thanks to the pumping action of the gravitational field, as suggested by Grishchuk well before the formulation of the conventional inflationary paradigm [1,2]. According to the quantum theory of parametric amplification (originally developed in the case of optical photons [3]), relic gravitons are produced from the vacuum or from a specific initial state with opposite momenta [4]. Squeezed states of relic gravitons are characteristic of many scenarios and, in particular, of conventional inflationary models [5,6]. The degree of quantum coherence of the large-scale correlations could be used to disambiguate their origin, at least in principle [7].
Gravitational waves might acquire an effective index of refraction when they travel in curved space-times [8,9] 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 [10]: 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. After the analysis of [10] the same idea has been pursued in two similar papers [11]. The first goal of the present analysis is a specific discussion of the relation between the approaches of Refs. [10] and [11] which ultimately coincide since they are related by a conformal rescaling involving the refractive index. We shall then explore a more general approach encompassing the ones discussed previously and analyze the actual phenomenological constraints together with the associated detectability prospects.
A refractive index varying in the early Universe is implicitly contained in the work of Szekeres [8,9] aiming at the general formulation of a macroscopic theory of gravitation. An effective index of refraction may also arise in modified theories of gravity violating the equivalence principle [10,11]. In the model-independent perspective pursued in Ref. [10] and subsequently followed in [11] (see, in particular, the second paper) the pivotal parameter of the discussion will therefore be the rate of variation of the refractive index in units of the Hubble rate, conventionally denoted by α and defined by: . (1.1) In Eq. (1.1) the overdot denotes a derivation 2 with respect to the cosmic time coordinate t; and c gw (t) are, respectively, the standard slow-roll parameter and the phase velocity of the relic gravitons in units of the speed of light. The phase velocity coincides, in the present case, with the group velocity and if the refractive index increases during inflation neither the phase velocity nor the group velocity will be superluminal. Note finally that the conventional slow-roll dynamics does not necessarily imply the validity of the so-called consistency relations which are instead broken by the presence of the refractive index so that the tensor spectral index and the tensor-to-scalar ratio are not solely determined by , as it happens when the consistency relations are not violated.
During the last year the Ligo/Virgo collaboration published the three-detector observations of gravitational waves from black-hole coalescence [12], the evidence of gravitational waves from neutron star inspiral [13], and the observation of a 50-solar-mass binary black hole coalescence for a redshift z = 0.2 [14]. As gravitational wave astronomy is opening a new window for the observation of the present Universe, the theory is beginning to consider more concretely the potential implications of the current developments for the early Universe. In the course of inflation, the evolution of the space-time curvature induces a stochastic background of relic gravitons [8,9] with a spectral energy density extending today from frequencies O(aHz) (i.e. 1 aHz = 10 −18 Hz) up to frequencies O(GHz) (i.e. 1 GHz = 10 9 Hz). In the conventional scenario the absolute normalization of the tensor power spectrum solely depends on the tensor to scalar ratio Here we shall be assuming the limit r T < 0.07 for the tensor spectral index from a joint analysis of Planck and BICEP2/Keck array data [15] (see also [16]). The typical energy density of inflationary gravitational waves between the mHz (i.e. 1 mHz = 10 −3 Hz) and the kHz is, optimistically, O(10 −16.5 ) in terms of Ω gw (ν, τ 0 ) which will denote throughout the spectral energy density in critical units at the present (conformal) time τ 0 . The corresponding chirp amplitude h c (ν, τ 0 ) will therefore be O(10 −29 ) for a typical comoving frequency of 0.1 kHz as it follows from the general relation: h c (ν, τ 0 ) = 1.263 × 10 −20 100 Hz ν h 2 0 Ω gw (ν, τ 0 ), (1.2) connecting the chirp amplitude to the spectral energy density in critical units. While the h c corresponding to the astronomical signals detected by the Ligo/Virgo collaboration is O(10 −21 ) [12,13,14], various improvements are expected at least in the so-called audio band between few Hz and 10 kHz. In this perspective it is useful to discuss potential deviations from the conventional inflationary signal hoping for larger spectral amplitudes over intermediate and high frequencies.
The plan of this paper is the following. In section 2 we shall demonstrate that the apparently different actions proposed so far to account for the same effects are indeed equivalent up to conformal rescaling of the four-dimensional metric. A general approach to the analysis of the power spectrum is proposed in section 3 where the tensor to scalar ratio and the relevant power spectra are explicitly derived. In section 4 the phenomenological constraints relevant to the present scenario will be studied by charting the corresponding parameter space; the potential signals will be confronted with the hoped sensitivities of the wide-band interferometers on the ground and in space. The concluding remarks are collected in section 5 and some useful details are relegated to the appendix A for the interested readers.

Actions and their parametrization
Let us start with an obvious remark and recall that, in the Einstein frame, the metric can be separated into a background value supplemented by the inhomogeneous contribution: where δ t g (E) µν denotes the tensor fluctuations of the geometry; in the present context the vector modes are anyway suppressed and shall be completely disregarded while the scalar fluctuations are described in terms of the gauge-invariant curvature inhomogeneities (see e.g. Eqs. (A.7)-(A.8) and discussion therein) but will not play a direct role in the forthcoming considerations since the refractive index is fully homogeneous.
Even if the analysis of the present section does not demand the absence of the intrinsic curvature, we shall now restrict the attention to the case where g (E) µν is conformally flat i.e.
where η µν = diag(+1, −1, −1, −1) is the Minkowski metric with signature mostly minus and the three-dimensional rank-two tensor h ij is both divergenceless and traceless. In Eq. (2.2) and throughout the rest of the paper τ denotes the conformal time coordinate. In the frame defined by Eq. (2.1) the action for the tensor modes of the geometry has been derived long ago by Ford and Parker in Ref. [17,18]: 3) can be supplemented by a second term containing the dependence on the refractive index n(τ ) which is the inverse of the phase velocity 3 in natural unitsh = c = 1, i.e. c gw (τ ) = 1/n(τ ): where P µν (E) is the conventional spatial projector tensor orthogonal to u (E) µ : 00 the action of Eq. (2.4) after straightforward modifications assumes the following form [10]: (2.6) 3 In the present case the phase velocity coincides with the group velocity since the relation between frequencies and wavenumbers is linear (as in the case of Refs. [10,11]) but the refractive index is time-dependent.
When the refractive index grows in time [10] c gw (τ ) = 1/n(τ ) decreases but what matters is that thanks either to a conformal rescaling or to a field redefinition the action (2.6) can be brought into the following generalized form: where A(τ ), B(τ ) and C(τ ) depend on the conformal time coordinate τ and are themselves functionals of either the scale factor, or of the refractive index or even of both [10]. A term like C(τ ) may appear by a banal field redefinition of the type a E h = h ij so we shall not dwell about this issue and focus hereunder on A(τ ) and B(τ ). After the results of Ref. [10], the authors of Ref. [11] suggested the following form of the action: 8) and argued that it leads to a blue spectrum of tensor modes 4 . It is relevant to point out that Eq. (2.6) has been used in Ref. [10] to argue that the spectra of relic gravitons may indeed increase in the presence of an inflating background. The authors of Ref. [11] did consider instead Eq. (2.8), reached a similar conclusion, but did not realize (or acknowledge) that Eq.
(2.8) is conformally related to Eq. (2.6). As a consequence Ref. [11] analyzed the case where the background is inflating but in the frame defined by Eq. (2.8). The choice of Ref. [11] is not particularly consistent and it is even confusing in the light of the previous analyses of Ref. [10]. We shall propose a generalized parametrization which may hopefully address and clarify the potentially confusing statements present in the current literature (see, in particular, Eq. (3.1) and discussion therein).
To show the conformal equivalence of Eqs. (2.6) and (2.8) let us now define a frame conformally related to Eq. (2.1) by positing Exactly as in the case of Eq. (2.1), the G µν can be decomposed into a homogeneous part supplemented by its own tensor inhomogeneities: In the frame defined by Eqs. (2.9) and (2.10) the tensor modes of the geometry obey: so that, thanks to Eqs. (2.12), (2.13) and (2.14) the action (2.1) in the conformally related frame becomes: If we now identify the conformal factor with the refractive index itself Ω(τ ) = 1 c gw (τ ) = n(τ ), (2.16) the action of Eq. (2.15) becomes which is exactly the action anticipated in Eq. (2.8) and discussed in Ref. [11]. All in all the action of Eqs. (2.4) and (2.6) (originally discussed in [10]) and the action (2.17) (subsequently used in Ref. [11]) are one and the same action since they are simply related by a conformal rescaling.
To make even more transparent the equivalence of the two descriptions it is useful to write down explicitly the evolution of the tensor modes in the two frames. From Eq. (2.8) it follows that the evolution of h (E) ij is given by: We must finally remark that n(τ ) and c gw (τ ) are completely homogeneous. While this is the simplest situation, if the spectral index is not fully homogeneous there will be additional scalar fluctuations mixing with the scalar adiabatic mode. At the end of appendix A it is demonstrated that the curvature perturbations share the same property of the tensor amplitudes since they are gauge-invariant and also invariant under conformal rescaling.
3 Generalized form of the tensor to scalar ratio

Diagonal form of the action
If we conventionally define the pivotal frame as the one where the background experiences a phase a quasi-de Sitter expansion, we should conclude that the choice of the pivotal frame is to some extent ambiguous, as already suggested. This ambiguity is not of a fundamental nature but simply stems from the current literature. Indeed, as already mentioned, in Ref. [10] the action (2.6) has been used to argue that the dynamics of the refractive index during an inflationary stage induces a growing spectrum of the tensor modes. After Ref. [10] the authors of Ref. [11] suggested instead the action (2.8) and considered (in their own framework) an inflating background. The authors of Ref. [11] did not acknowledge the results of Ref. [10] and did not notice or appreciate that the actions (2.6) and (2.8) are one and the same action, up to a conformal rescaling: this is the main reason of the frame ambiguity mentioned above. It is not a surprise that if we impose that the background in inflating in both frames defined by Eqs. (2.6) and (2.8), growing spectra of the tensor modes will be obtained in both frames.
In what follows the aim will be first to quantize the problem in general terms and then to deduce the spectral energy density in various cases that have not been addressed so far. For this purpose and to account of the perspectives of Refs. [10,11] it is convenient to analyze simultaneously the two situations by considering the following diagonal form of the action for the tensor modes 6 : where γ is a further parameter varying between 0 and 1. If γ = 0 the action (3.1) reproduces exactly Eq. (2.6) with c gw (τ ) = 1/n(τ ). Conversely when γ = 1 Eq. (3.1) reproduces Eq. (2.17) (or Eq. (2.8)). According to a purely heuristic logic the values of γ may also lie outside the interval 0 ≤ γ ≤ 1: in this perspective γ can be considered as a further parameter that simply rescales the tensor spectral index from its value for γ = 0. However, if γ = 0 the sign of the tensor spectral index does not change, as we shall see. The presence of the refractive index in the second term inside the squared bracket suggests that the action can be diagonalized in terms of a new time coordinate, namely dτ = n(η)dη. In the η-parametrization the action of Eq. (3.1) becomes: While the relation of b(η) to a(η) and n(η) may change, what matters is the overall dependence of the power spectrum on the background fields. The parametrization of Eq. (3.2) will be used to demonstrate, in a unified approach to the problem, that both Eqs. (2.8) and (2.17) lead to a growing power spectrum for typical wavelengths larger than the Hubble radius.

Quantization and power spectrum
For a correct normalization of the power spectrum the quantization of the system is essential and for this purpose the η-parametization is more convenient, as already emphasized in Ref.
[10] without specific details. This gap will now be bridged by only using the action of Eq. (3.2) for the explicit derivation of the power spectrum, of the tensor to scalar ratio and of all the other relevant quantities. We start by noticing that from the explicit form of the tensor polarizations given in Eq.
where the overdot denotes a derivation with respect to η, i.e.ḣ λ = ∂ η h λ . From the definition of the canonical momenta, the Hamiltonian in the η-parametrization becomes 7 : The canonical form of the commutation relations at equal η-times and the explicit form of the Hamiltonian (3.5) lead directly to the evolution equations forĥ λ andπ λ : Using Eq. (3.7) the Fourier representations of the operatorsĥ λ andπ λ is: The evolution equations of the mode functions can be decoupled as: where the polarization index has been omitted since the result of Eq. (3.12) holds both for ⊕ and for ⊗. The quantum mechanical initial conditions correspond to (3.13) In the limit n(η) → 1 the η-time turns into the conformal time while b(η) → a(τ ): therefore the initial data of Eq. (3.13) generalize the standard quantum mechanical initial condition to the case where the refractive index is dynamical. Note that since b(η) has a power-law dependence for η < −η * (see below Eqs. (3.23) and (3.24)), Eq. (3.12) can be exactly solved in the η-time (see Eq. (3.25)) with the initial conditions dictated by Eq. (3.13).
In conclusion we have that the full mode expansion of the tensor amplitudeĥ ij ( x, η) in the η-time parametrization easily follows from the combination of Eq. (3.3) and (3.8) and it is given by:ĥ 14) The two-point function computed from Eqs. (3.14) is simply 8 where j 0 (kr) is the spherical Bessel function of zeroth order [22,23]. The tensor power spectrum of Eq. (3.15) is then given by

Explicit evolution of the mode functions and related matters
As already mentioned in Eq. (1.1) the variation of the refractive index is measured in units of the Hubble rate and for the sake of concreteness we shall consider the general situation where while n(a) = 1 for a > a * . It is not difficult to imagine various continuous interpolation between the two regimes but what matters is the continuity of n(a), not the specific form of the profile across the normalcy transition. One of the simplest possibilities is given by 9 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. These profiles are illustrated in Fig coincide with the end of the inflationary phase. The value of a * corresponds with 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 interesting to consider the case N * < N t or even N * N t . It is finally useful to bear in mind that the scale factor during the different stages of the evolution are continuous and differentiable. An explicitly continuous parametrization of the scale factors is given by:

18)
9 Note that n i ≥ 1 but we shall always consider the case n i = 1 as representative of the general situation.
where τ 1 coincides with the end of the inflationary phase and τ 2 coincides with the time of matter-radiation equality; note that β → 1 in the case of a pure de Sitter phase and β = 1−O( ) in the quasi-de Sitter case. The transition to the domination of the dark energy does not affect the slope and it has a mild effect on the amplitude of the spectrum [19,20,21]; this effect will be included in the actual estimates of section 4 by following the approach already pursued in [21]. A similar comment holds for other effects which are customarily associated with potential suppression of the amplitude of the power in the intermediate frequency branch such as the neutrino free streaming [24,25] or the effective evolution of number of relativistic species. Equations (3.18), (3.19) and (3.20) are all continuous with their first derivatives at the transition points 10 . If the scale factor is continuous at the transition points and if n(a) is continuous and differentiable as a function of the scale factor, then the relation between η and τ is also continuous and differentiable. All in all we can say that whenever n(a) is a well behaved function interpolating between the two regimes illustrated in Fig. 1 it will always be possible to consider the evolution as piecewise continuous without worrying about the evolution in the transition regime. By then recalling Eq. (3.2) and the explicit relation between η and τ , we have that for a < a * the η-time and the conformal time are related as: so that, in the conformal time parametrization, b(τ ) turns out to be: . (3.23) In the regime of Eq. (3.24) we can easily solve Eq. (3.12) by defining the auxiliary mode function F k (η)b(η) = f k (η). The resulting equation can be solved in terms of Hankel functions and F k (η) is then given by 11 : where H (1) µ (kη) is the Hankel function 12 of the first kind [22,23]. Note that with the characteristic choice of the phase appearing in N the limit of Eq. (3.25) for η → −∞ is exactly the one already mentioned in Eq. (3.13).

Power spectra and tensor to scalar ratio
Inserting Eqs. (3.24) and (3.25) into Eq. (3.16) the explicit form of the tensor power spectrum becomes: In the simplest situation η * = η 1 (i.e. N * = N t ) and when the relevant modes are larger than the Hubble radius Eq. (3.26) can be modified by recalling Eq. (3.22); in this limit we obtain [22,23]: Since the limits k aH and kη 1 describe the same physical regime, after some algebra Eq. (3.27) becomes 13 : From Eq. (3.28) the tensor spectral index is defined as n T = 3 − 2µ; using then Eq. (3.25) the explicit expression of n T depends on α, β and γ in the following manner: The result of Eq. (3.28) can be further modified by including the slow-roll corrections 14 and by recalling, from Eq. (3.22), that η 1 and τ 1 are related as The solution (3.25) can be expressed in the τ -time by using Eq. (3.22); depending on the convenience, both parametrizations can be consistently discussed. Note, finally, that when β = 1, α = 0 and γ = 0 the index µ equals 3/2, as expected. 12 This solution reminds of the standard solution for the mode function; note, however, that in the present case the argument of the Hankel function is not kτ but rather kη. In terms of the τ -parametrization the evolution is more complicated but equally solvable. 13 The η-dependence in Eq. (3.28) becomes: (3.30) Denoting with Ω R0 the present value of the critical fraction of radiative species (in the concordance paradigm photons and neutrinos) and with A R the amplitude of the scalar power spectrum at the pivot scale (see Eq. (A.8)) the value of k max in Mpc −1 units is: In Eq. (3.31) the δ accounts for the possibility of a delayed reheating terminating at an Hubble scale H r smaller than the Hubble rate during inflation. In what follows, as already mentioned in section 1, we shall rather stick to the conventional case where the reheating is sudden and δ = 1/2 (or H 1 = H r since the end of the inflationary phase coincides with the beginning of the radiation epoch 15 ).
As suggested in Ref. [10], whenever α > 0 the spectral index gets blue since it is always greater than the slow-roll correction. Using Eq. (3.29) this statement can be made more precise by expanding n T (α, β, γ) in the limit < 1: When γ = 0 we are in the situation of Ref. [10] and Eqs. (3.32)-(3.33) imply that the tensor spectral index is: When γ = 1 we are in the situation of Ref. [11] and Eqs. (3.32)-(3.33) implies instead: it is immediate to appreciate that the spectrum is blue and n T > 0 provided: these two conditions are always verified in practice by definition of slow-roll parameter (implying 1). All in all we can then say that within the physical range of the parameters the presence of γ just shifts the value of the tensor spectral index (see Eqs. (3.34) and (3.35)) but does not change qualitatively the whole picture 16 .
If the transition to normalcy is sufficiently sudden, the details of the transition will be immaterial for the large-scale spectrum. We then have that the tensor to scalar ratio is given by 17 : (3.37) Equations (3.32) and (3.37) demonstrate that the consistency relations do not hold in the case where the gravitons have a refractive index. A similar phenomenon has been suggested in [26] as a consequence of the opposite physical situation where the initial thermal bath of primordial phonons (i.e. large-scale quanta of curvature perturbations) is characterized by a sound speed different from the speed of light. The large-scale power spectra are not modified in the large-scale limit as long as the evolution of b(η) is continuous. For this purpose, without indulging in the explicit analytic expression in the transition region, we can consider the situation where before −η * the large-scale power spectrum is given by the solution (3.25). The resulting power spectrum at a later time is therefore directly computable since, in the large-scale limit, Eq. (3.12) becomes 18 ∂ η [b 2 (η)∂ η F k (η)] = 0. This equation can then be explicitly integrated from −η 2 onwards, namely: where the time −η 2 is a generic reference time before −η * and very close to it. With this observation we have that the power spectrum after −η * can be written as (3.39) Equation (3.39) shows that the only thing that counts is the continuity of b(η) across the transition. Well after −η * we also have that b(η) → a(η) = a(τ ) and the background still evolves like quasi-de Sitter with slow-roll parameter given by . It is then clear that the second term inside the square bracket is subleading and this shows, as expected, that the large-scale power spectrum at large-scale is not affected by this kind of continuous transitions. 16 From Eq. (3.33) we have that n T > 0 iff γ < 3/2: it is interesting to appreciate that large values of γ lead to a red spectrum so that, even in the heuristic approach mentioned after Eq. (3.1) (which is not the one followed here) we should always demand 0 ≤ γ < 3/2. 17 Denoting with η i the initial time of the evolution (for instance at the onset of inflation) we shall have that n i = n 1 (a i /a 1 ) α = n 1 e −αN * with n i ≥ 1. When n i = 1, n 1 → exp (αN * ) (see also Fig. 1 and discussion  therein). 18 Note that in the large-scale limit the term k 2 F k (η) is negligible in Eq. (3.12).

The spectral energy density
The independent variation of the tensor spectral index n T and of the tensor to scalar ratio r T (k p ) can be investigated from the absolute normalization of the spectral energy density [21]. When the relevant modes are today inside the Hubble radius 19 the relation between the spectral energy density and the power spectrum at the corresponding time is given by: (3.40) There are two possibilities for an accurate assessment of the absolute normalization appearing in Eq. (3.40). In the first approach we can compute the transfer function for the amplitude as [27,21,28,29,30] where Ω M0 denotes the present critical fraction of dusty matter in the concordance paradigm. By using the transfer function for the tensor amplitude, the spectral energy density for frequencies ν ν eq can be simply given by: where d A is the (comoving) angular diameter distance to decoupling and Ω R0 is the fraction of critical energy density attributed to massless particles (photons and neutrinos in the concordance paradigm). The second strategy is to use the direct integration of the mode functions and discuss the transfer function directly in terms of the spectral energy density [21,29]: where ν eq denotes the equality frequency: where T ρ (ν/ν eq ) is the transfer function of the spectral energy density: T ρ (ν/ν eq ) = 1 + c 2 ν eq ν + b 2 ν eq ν 2 , c 2 = 0.5238, b 2 = 0.3537. , the amplitude for ν ν eq differs, roughly, by a factor 2. This coincidence is not surprising since Eqs. (3.43)-(3.44) have been obtained by averaging over the oscillations (i.e. by replacing cosine squared with 1/2). These manipulations are less accurate than the procedure used to derive the transfer function for the spectral energy density, as stressed in the past [21].
Recalling Fig. 1 and the related discussion, when the critical number of efolds equals the total number of efolds (i.e. N * = N t ) the form of the spectral energy density is the same as in the standard case with the difference that the tensor to scalar ratio and the spectral index dapend on the parameters of the model and are given, respectively, by Eqs. (3.29) and (3.37). If N * < N t the spectral energy density is characterized by two branches. More specifically the full expression of the spectral energy density can be written as: where ν * and ν max are given, respectively, by: In general terms, as we shall see, ν * will always be larger than O(10 −11 ) Hz but, depending on the values of α and N * it can be either larger or smaller than the Hz.

Phenomenological implications
The rate of variation of the refractive index α and the critical number of efolds N * are the two pivotal parameters for the analysis reported in the present section. If N * is of the order of N t the transition to normalcy occurs at the end of inflation; conversely when N * < N t (as in the case of Fig. 1) the transition to normalcy takes place well before the onset of the radiation-dominated epoch (i.e. when the background is still inflating deep inside the quasi-de Sitter stage of expansion). These two different situations will now be separately examined by assuming, for illustrative purposes, a fiducial number of inflationary efolds N t = O(70). Larger values of N t do not affect the main conclusions inferred below.

Basic constraints
The scale-dependent constraint on the tensor to scalar ratio is conventionally applied at the pivot scale k p = 0.002 Mpc −1 corresponding to the comoving frequency ν p = k p /(2π). While the Planck data alone would suggest an upper limit r T (k p ) < 0.11, the joint analysis of Planck and BICEP2/Keck array data [15] already demands r T (k p ) < 0.07. In a conservative perspective we shall therefore require the limit:  nucleosynthesis [31,32,33] imply instead a constraint on the integral of the spectral energy density of the gravitons: (4. 2) The lower limit of integration in Eq. (4.2) is given by the frequency corresponding to the Hubble rate at the nucleosynthesis epoch: where g ρ denotes the effective number of relativistic degrees of freedom entering the total energy density of the plasma and T bbn is the temperature of big-bang nucleosynthesis. The limit of Eq. (4.2) sets an indirect constraint on the extra-relativistic species (and, among others, on the relic gravitons). Since Eq. (4.2) is relevant in the context of neutrino physics (when applied to massless fermionic species), the limit is often expressed for practical reasons in terms of ∆N ν representing the contribution of supplementary neutrino species. The actual bounds on ∆N ν range from ∆N ν ≤ 0.2 to ∆N ν ≤ 1; the integrated spectral density in Eq. (4.2) is thus between 10 −6 and 10 −5 . Finally the pulsar timing measurements determine a local bound on the spectral energy density [34,35] Ω gw (ν pul , τ 0 ) < 1.9 × 10 −8 , ν pul 10 −8 Hz, (4.4) which is an upper limit at a typical frequency, roughly corresponding to the inverse of the observation time along which the timing has been monitored. Various applications and refinements of the pulsar timing cosntraint have been discussed through the years and also recently (see e.g. Refs. [36,37,38,39]).

Target sensitivities and speculations
Generally speaking the phenomenological constraints of Eqs.  .2) and (4.4) are satisfied in the shaded area which has been obtained by assuming (incorrectly, as we shall see in a moment) that the spectral index and the tensor to scalar ratio change independently. The labels appearing on the curves of the plots of Figs. 2 and 3 denote the common logarithm of h 2 0 Ω gw at the four frequencies indicated above each of the plots. The two frequencies discussed in Fig. 2 correspond, respectively, to 10 and 100 Hz: they illustrate the spectral energy density in the so-called audio band ranging between few Hz and 10 kHz. Terrestrial wide-band interferometers, depending on their respective physical dimensions, operate within this band. The two frequencies illustrated in Fig. 3 fall instead within the mHz band ranging between a fraction of the mHz and few Hz. Space-borne detectors operate within this band. We can finally define a MHz band, extending between 100 kHz and few GHz; this band is irrelevant for the signal 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.
While the current upper limits on stochastic backgrounds of relic gravitons are still far from the final targets [40,41], the advanced Ligo/Virgo projects are described in [42,43]. We shall then suppose, according to Refs. [42,43] 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 two further interferometers (with different designs) operational in the audio band namely the Japanese Kagra 20 (Kamioka Gravitational Wave Detector) [51,52] and the Einstein telescope [53] both with improved sensitivities in comparison with the advanced Ligo/Virgo interferometers.
The target sensitivity to detect the stochastic background of inflationary origin should be h c = O(10 −29 ) (or smaller) in the chirp amplitude and h 2 0 Ω gw = O(10 −16 ) (or smaller) in the spectral energy density. These figures 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 hopefully operational in the future. The space-borne interferometers, such as (e)Lisa (Laser Interferometer Space Antenna) [44], Bbo (Big Bang Observer) [45], and Decigo (Deci-hertz Interferometer Gravitational Wave Observatory) [46,47], might operate between few mHz and the Hz hopefully within the following score year. While the sensitivities of these instruments are still at the level of targets, we can say that they should probably range between h 2 0 Ω gw = O(10 −12 ) and h 2 0 Ω gw = O(10 −15 ). In spite of the minute signals related to conventional inflationary scenarios, if the spectrum increases in frequency (i.e. n T > 0) Figs. 2 and 3 seem to suggest potentially large signals both for terrestrial interferometers (in their advanced version) and for space-borne detectors. However, this way of reasoning is misleading, at least partially. Indeed, while it is true that in the present situation the consistency relations are violated, it is equally true that r T and n T depend on , N * , α and γ in a specific way (see e.g. Eqs.  Consider first the regions r T < 0.06 (from Eq. (3.37)) and n T > 0.4 (see Eq. (3.32)): Fig.  4 shows that these two corners of the parameter space are mutually exclusive. If we relax a bit the requirement on the spectral index and impose n T > 0.3 the two regions get a bit closer (see right plot in Fig. 4). Whenever N t = O(N * ) the potential signals due to relic gravitons in the audio band are excluded and this result holds on spite of the value of γ. By looking at Fig. 5 we see that the only appreciable overlap between the areas n T > 0.07 and r T < 0.06 is a tiny slice representing the allowed region of the parameter space. By going back to Fig.  3 we see that the spectral energy density corresponding to the allowed region of Fig. 5 is . These values of the spectral energy density will be hopefully probed by Bbo/Decigo interferometer. We therefore conclude that in the case N t = N * the variation of the refractive index does lead to a blue spectrum which is however only relevant in the mHz band (but not in the audio band). In particular in the mHz band this signal can be larger (by a couple of orders of magnitude) than the spectral energy density coming from the conventional inflationary models.

Charting the parameter space in the N * N t case
When N * < N t (or even N * N t ) the spectral energy density consist of three branches: the usual infrared branch is supplemented by an intermediate regime (where the spectral slope is blue) terminating with a quasi-flat plateau at high frequencies. In the high-frequency region the spectrum is quasi-flat with a red slope determined by the value of . The frequency ν * divides  the intermediate branch of the spectrum from the high-frequency tail and it only depends on N * , N t and α. There are two complementary physical possibilities. If ν bbn < ν * < Hz the transition regime takes place between the frequency of the big-bang nucleosythesis and the mHz band where space-borne interferometers will be operating; this region correspond to the shaded area appearing in the left plot of Fig. 6. In the complementary regime we have instead that ν * > Hz and the break of the spectrum is in the audio range. This region corresponds to the shaded area appearing in the right plot of Fig. 6. As long as 25 < N * < 50 and 0 < α < 2, the results of Fig. 6 confirm that ν bbn < ν * < Hz, as expected from qualitative arguments. For the same range of α the condition ν * > Hz will be verified for a comparatively larger N * .
We are now in conditions of charting the parameter space in the (α, N * ) plane by requiring that the constraints of Eqs. (4.1) Eqs. (4.2) and (4.4) are simultaneously satisfied. In Fig. 7 we enforced all the relevant bounds and also required 22  of γ can also be considered continuous in the purely heuristic perspective invoked after Eq. (3.1), we shall just illustrate the cases γ = 0 and γ = 1 corresponding to the two equivalent parametrizations of the action discussed in sections 2 and 3. In this connection we can remark that in the case γ = 0 (left plot in Fig. 7) the allowed region is larger than in the case γ = 1 (right plot in Fig. 7). Note that in Fig. 7 the frequency ν a = 0.1kHz has been chosen as representative of the typical signal in the audio-band where the sensitivity of wide-band detectors is approximatively larger. A different choice in the range Hz < ν a < kHz will lead to similar conclusions 23 . In Fig. 8 the same analysis of Fig. 7 is extended to mHz band. The dashed areas in Fig. 8 have the same meaning of those depicted in Fig. 7 with the difference that we are now imposing a more stringent requirement on the spectral energy density namely:  If we would require a larger sensitivity in the mHz band (such as the the ones often quoted for the Bbo or Decigo projects) we would be led to demand, for instance, In Fig. 9 the Bbo/Decigo target sensitivities are illustrated in the (α, N * ) plane. As it is clearly visible from the direct comparison of Figs. 8 and 9 the dashed areas get larger. Still the allowed region in the γ = 1 case is always larger than for γ = 0. By comparing Figs. 7, 8 and 9 in a unified perspective, the sweet spots now get larger if we demand that space-borne interferometers will be able to probe values of the spectral energy density as low as the ones of Eq. (4.7). In particular the dashed area of Fig. 7 is now all contained within the dashed area of Fig. 9: this means that by remaining in the dashed slice of Fig. 7 the resulting signal will be, a fortiori, visible both in the audio band and in the mHz band. Even if the target sensitivities of Eq. (4.7) look by definition quite optimistic they might not be totally insane if we bluntly assume that the minimal detectable chirp (or strain) amplitudes are the same in the mHz and in the audio band. The argument goes in short as follows. A given chirp (or strain amplitude) of the signal in the audio band (say for ν = ν LV = O(0.1)kHz) probes a spectral energy density which is comparatively larger than the one achievable, with the same sensitivity, in the mHz band. To appreciate this we can consider, for instance, Eq. (1.1) and suppose that the spectral energy density of the signal will be, in critical units, h 2 0 Ω gw (ν a , τ 0 ) = O(10 −10 ). The required sensitivity in terms of h c will therefore be of the order of 10 −25 . Let us now pretend that the achievable sensitivity (in terms of the chirp amplitude and in the mHz band) is the same as in the audio band, i.e. h c = O(10 −25 ). Since the frequency in the mHz band is O(10 −3 ) Hz, Eq. (1.1) can be inverted in order to obtain the minimal spectral amplitude detectable in the mHz band; the result will be h 2 0 Ω gw (ν a , τ 0 ) = O(10 −15 ). We shall now select the values of α and N * within the dashed areas of Figs. 7, 8 and 9 and present some interesting example. In Fig. 10 we illustrate the spectral energy density by fixing α = 0.14 (for γ = 0); similarly in Fig. 11 we demand α = 0.42 (for γ = 1). Both in Figs. 10 and 11 we choose different values of N * (i.e. N * = 40, 50, 60). As N * → N t the results discussed at the beginning of this section are recovered: the high-frequency plateau gets narrower and asymptotically disappears. Conversely, when the value of N * decreases, ν * also diminishes and the quasi-flat branch of the spectrum gets wider, as it can be seen from Figs. 10 and 11. The other values have been selected within the allowed regions of the parameter space and they are α = 0.14 and γ = 0.
Let us now consider more closely the values of α and γ. In the left plot of Fig. 9 (i.e. γ = 0) the value of α = 0.14 is within the shaded region and the same is true for α = 0.42 in the right plot (i.e. γ = 1) of the same figure. The value of α in Fig. 10 is three times smaller than the one of Fig. 11 and this occurrence has a simple explanation 24 which is directly related to Eqs.  Figure 11: The spectral energy density is illustrated for N * < N t ; to compare this plot with the one of Fig. 10 we stress that α = 0.42 and γ = 1.
n T in the cases γ = 0 and γ = 1, the spectral energy density in the case γ = 1 is smaller than in the case γ = 0 for the same choice of N * . All in all, when N * < N t the spectral energy density develops two branches for ν > ν eq . The intermediate branch corresponds to those modes leaving the Hubble radius during inflation when the refractive index is still dynamical. These modes will then lead to an increasing spectral slope for ν eq < ν < ν * . The second branch of the spectrum corresponds to those modes that left the Hubble radius when the refractive index already ceased to be dynamical: in this branch the spectrum decreases as ν −2 for ν * < ν < ν max .In Fig. 10 the audio branch of the spectrum and the mHz branch have been illustrated, respectively, with trapeziodal and triangular shapes encompassing the respective frequency ranges. The heights of the two triangles roughly matches the hoped sensitivities of (e)Lisa and Bbbo/Decigo. Similarly the height of the trapezoid in the audio branch approximately correspond to the sensitivities of terrestrial wide-band detectors (in their advanced version). Not all the curves of Figs. 10 and 11 match the required sensitivities of future detectors.
In the (α, N * ) plane there exist regions that are potentially visible both in the audio band and in the mHz bands. These regions lie within the shaded areas of Figs. 7 and 9 and are illustrated in Fig. 12 in the case γ = 0. The trapezoidal shape of Fig. 12   encompasses the best sensitivities in the spectral energy density for terrestrial and space-borne detectors in the audio and mHz bands. The parameters illustrated in Fig. 12 belong to the same fine-tuning line: if the parameters lie along that line in the (α, N * ) plane we can be reasonably sure that the corresponding spectral energy density in critical units will be potentially visible both by terrestrial and by space-borne detectors. The parameters leading to a comparatively large spectral energy density in the case γ = 0 correspond to a much smaller value of h 2 0 Ω gw in the case γ = 1. This conclusion is clear by comparing the two plots of Fig. 12: the two plots share exactly the same parameters but correspond to different values of γ. As explained above in connection with Figs. 10 and 11 this feature has a simple explanation since different values of γ = 0 simply shift the tensor spectral index of the intermediate branch in comparison with the case γ = 0.

Concluding remarks
We showed that a rather generic scenario for the evolution of the refractive index during the early stages of a conventional inflationary phase leads to power spectrum of relic gravitons naturally tilted towards frequencies. The spectral slope is determined, in this context, by the competition between the rate of variation of the refractive index and the expansion rate of the quasi-de Sitter stage. We argued that it is always possible to redefine the time coordinate and to compute the evolution of the field operators in terms of an effective horizon incorporating the dynamical effects due to the variation of the refractive index.
The slope of the tensor spectrum differs for the modes crossing the effective horizon in the subsequent phase where the gravitons propagate again at the speed of light and the refractive index goes back to 1. In this branch the spectral slope is quasi-flat but it decreases at a rate that depends on the slow-roll parameter. If the evolution of the refractive index stops at the end of inflation, the spectral energy density consists of the usual low-frequency branch (between few aHz and 0.1 fHz) supplemented by a growing branch extending up to the MHz or even GHz regions. Conversely if the evolution of the refractive index terminates before the end of inflation the previous branches are supplemented by a quasi-flat plateau at high frequencies. The frequency range of the plateau depends on the rate of variation of the refractive index (in units of the expansion rate) and on the critical number of efolds during which the refractive index evolves. After imposing the usual phenomenological constraints, it turns out that the plateau may extend from the mHz band to the GHz band. In a sweet spot of the parameter space all the known bounds are satisfied and the resulting signal could be detectable, in principle, either by terrestrial interferometers (in their advanced and enhanced configuration) or by space-borne detectors.

A Tensor, scalars and their frame-invariance
The two polarizations of the gravitational wave are defined as: The tensor power spectra P T (k, η) and Q T (k, η) determine the two-point function at equal times: ĥ ij ( k, η)ĥ mn ( p, η) = 2π 2 k 3 P T (k, η) S ijmn (k)δ (3) ( k + p), (A.5) where P T (k, η) is the power spectrum in the η parametrization while: ij (k) e (λ) mn (k)/4 = 1 4 p mi (k)p nj (k) + p mj (k)p ni (k) − p ij (k)p mn (k) . (A.6) Equation (A.5) follows the same conventions used when deriving the spectrum of curvature perturbations on comoving orthogonal hypersurfaces (customarily denoted by R( x, τ )) which is the quantity employed to set the initial conditions for the evolution of the temperature and polarization anisotropies of the Cosmic Microwave Background. According to the standard convention, the scalar power spectrum is assigned as: where k p is called the pivot scale, n s is the scalar spectral index and A R is the amplitude of the scalar power spectrum at the pivot scale. The conformal rescaling of Eq. (2.1) assumes that the refractive index is homogeneous. If this is not the case the results for the tensor modes will be exactly the same while the scalar modes will be modified since the fluctuations of the conformal factor contribute to the scalar and not to the tensor problem. More specifically, the scalar modes of the geometry in the Einstein frame: will be related to the scalar modes in the conformally related frame via the scalar fluctuation of Eq. (2.9), i.e. δ s g (E) µν = δ s n G  .16) imply H E = H. However, the mismatch between ψ E and ψ is exactly compensated by the mismatch between H E and H. In fact, from Eq. (2.13) we have 2(H E − H) = n /n so that Eq. (A.13) 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.