Flat spectra of cosmic gravitons in the nHz and audio bands

The spectra of the relic gravitons are customarily normalized in the low-frequency domain where the signal of the concordance paradigm is expected to peak and this is why their contribution to the temperature and polarization anisotropies of the microwave background is only described by the tensor to scalar ratio. If the consistency relations are broken, the same strategy is accomplished by introducing the tensor spectral index as a further independent parameter. When the dominant component of the spectral energy density is distributed for frequencies much larger than the aHz, the logic behind this conventional approach is much less compelling. The improved bounds in the audio band and the current data from the pulsar timing arrays in the nHz region motivate a new strategy for the absolute normalization of the cosmic background of relic gravitons. After introducing a general four-dimensional action for the analysis of the relic gravitons the new approach is illustrated in the case of conventional and unconventional inflationary models.


Introduction
Since the evolution of the tensor modes of the geometry is not Weyl-invariant [1,2,3] the production of relic gravitons is expected, with different phenomenological signatures, in a variety of scenarios and, in particular, during an isotropic phase of quasi-de Sitter expansion [4]. Besides the vanilla ΛCDM paradigm 2 , the simplest version of the concordance scenario includes only one further free parameter, namely the ratio r T (k p ) describing the tensor component of the large-scale inhomogeneity at a conventional pivot scale (see e.g. [5,6,7]) that coincides, in what follows, with k p = 0.002 Mpc −1 . This scale corresponds 3 to a comoving frequency ν p = k p /(2π) = 3.09 aHz and this is why, in the conventional lore, the limits on r T (k p ) translate into constraints on the spectral energy density of the relic gravitons in the aHz range. The spectral energy density in critical units (and at the present time) is denoted hereunder by Ω gw (ν) ≡ Ω gw (ν, τ 0 ) where ν is the comoving frequency. In the concordance paradigm Ω gw (ν) approximately scales as ν −2 between the aHz and 100 aHz [8,9,10] while it is flat (or slightly decreasing) for larger frequencies.
Two tacit assumptions are implicitly associated with r T : the first one is that the early completion of the concordance paradigm is provided by the conventional inflationary scenario, the second is that the consistency relations are not violated 4 . Between these two assumptions (which are rarely stressed) the former seems stronger than the latter but they are instead equally essential if the only final objective is a stringent set of bounds on r T from the temperature and polarization anisotropies of the microwave background. In practice the consistency relations stipulate that the tensor spectral index n T and the slow-roll parameter are both determined by the value of r T according to the following (approximate) chain of equalities n T −r T /8 −2 . The consistency relations are valid in the case of single-field inflationary scenarios (see e.g. [11]) but they can be otherwise broken. Since the vanilla ΛCDM only represents a useful compromise between the available data and the number of ascertainable parameters, the addition of a tensor component (only described by r T ) allows for a rather accurate set of limits implying, in a conservative perspective, that r T ≤ 0.06 [5,6,7]. If the consistency conditions are broken the accuracy on r T becomes comparatively smaller and, for this reason, some time ago it has been suggested that the bounds on r T must be viewed in a broader perspective where the analysis of the three standard cosmological data sets (i.e. cosmic microwave background anisotropies, large-scale structure and supernovae), is combined with the bounds of wideband interferometers in the so-called audio band (i.e. between few Hz and few kHz) [12].
Since direct measurements are now available it seems appropriate to reconsider the high-frequency normalization of the cosmic backgrounds of relic gravitons as a possible alternative to the conventional approach based on the analysis of the aHz region. There are two complementary motivations corroborating this suggestion and they are associated with the recent claims of the pulsar timing arrays (PTA) in the nHz band and with the improved constraints provided by the wide-band interferometers in the audio band. Indeed, in the last thirty years the pulsars provided a series of relevant upper limits on the 2 Λ stands for the dark-energy component while the CDM denotes the Cold Dark Matter contribution. The simplest scenario where the neutrinos are massless, the dark-energy does not fluctuate and the tensor modes are absent is customarily referred to as the vanilla ΛCDM model. 3 The unitsh = c = 1 will be consistently employed throughout. Furthermore in this investigation the scale factor is normalized as a(τ0) = a0 = 1. Finally, the standard prefixes shall be consistently employed, e.g. 1 aHz = 10 −18 , Hz, 1 mHz = 10 −3 Hz, 1 MHz = 10 6 Hz and so on and so forth. 4 For the sake of conciseness the argument of rT (k) will be dropped when not strictly necessary and we shall therefore employ the following shorthand notation rT = rT (kp).
spectral energy density of the relic gravitons at intermediate frequencies [13,14,15,16,17]. Broadly speaking the previous results 5 suggested h 2 0 Ω gw (ν, τ 0 ) < 10 −10 for ν = ν P = O(nHz). More recently the PTA reported a series of effects that could be attributed to the relic gravitons [18,19,20,21]. Barring for the slight differences between results reported by the four collaborations, the current evidence might suggest that h 2 0 Ω gw (ν, τ 0 ) = O(10 −9 ) for a narrow slice of frequencies approximately ranging between a fraction of the nHz and 30 nHz. The observations of Refs. [18,19,20,21] are, at the moment, very preliminary and the key property of a PTA is that the signal from relic gravitons, unlike a potential noise, is correlated across the baselines. These correlations have not been observed but it is nonetheless interesting to consider more carefully the relic gravitons in the nHz domain; this is the perspective already conveyed some years ago [22,23] before the evidences of the PTA. The potential signals in the nHz band must anyway be complemented by a series of more robust limits in the audio band coming from the direct measurements of wide-band detectors. Depending on the slope of the spectral energy density, the bounds from wide-band interferometers slowly improved through the years [24,25,26,27,28,29,30] and finally led to the joint analysis of the LIGO, Virgo and KAGRA collaborations [31] suggesting that we can even have Ω gw (ν, τ 0 ) ≤ O(10 −9 ) for typical frequencies between 20 Hz and 80 Hz. The current measurements in the nHz and in the audio band seem then to point towards a quasi-flat spectral energy density of the relic gravitons when the comoving frequency encompasses the nHz and the audio bands.
In this paper we argue that even in the context of conventional inflationary models, the spectral energy density of the relic gravitons can be quasi-flat at high-frequency with typical amplitudes much larger than the ones of the concordance scenario. The potential signal must then be compatible, in this context, with the PTA observations and consistent with the limits of wide-band interferometers in the audio band. Instead of worrying about r T (as customarily done in the standard approach [5,6,7]) we can directly impose the high-frequency normalization 6 . To pursue this possibility we first propose a general four-dimensional action for the consistent analysis of the relic gravitons evolving in conformally flat background geometries and generalizing the Ford-Parker action of Ref. [3]. After introducing the different parametrizations of the action and its quantization, the conditions leading to a flat slope at high-frequency are analyzed by computing the spectral energy density within the Wentzel-Kramers-Brillouin (WKB). The strategy for the high-frequency normalization is then explained by considering the PTA evidences, the limits from the wide-band detectors and the other constraints customarily associated with the high-frequency gravitons. In short the layout of the paper is the following. In section 2 the basic action of the problem is analyzed in its different forms. As suggested by the current data, the normalization of the spectral energy density in the nHz and audio bands is discussed in section 3. Section 4 is focused on the conditions for a flat high-frequency spectrum with amplitudes potentially much larger than the signal of the vanilla ΛCDM scenario. In section 5 we present some phenomenological considerations that complement the results of section 4. Finally section 6 contains the concluding remarks. 5 As usual h0 is the Hubble rate expressed in units of 100 Hz km/Mpc and since Ωgw(ν) denotes the spectral energy density in critical units, h 2 0 appears in its denominator. For this reason it is common practice to phrase the discussions directly in terms of h 2 0 Ωgw(ν) that is independent of the specific value of h0. 6 The approach discussed here does not exclude that rT is drastically smaller than the current bounds stemming from the temperature and the polarization anisotropies of the microwave background.
2 General parametrization of the action

The general action and its parametrizations
In the single-field case the effective action of generic inflationary models involves all the terms that include four derivatives and are suppressed by the negative powers of a large mass scale [32]. In nongeneric models of inflation the higher-order corrections may assume a specific form either because the inflaton has a particular symmetry or because the rate of inflaton roll remains constant (and possibly larger than 1). Examples along this direction are certain fast-roll scenarios [33,34,35] or the higher-order curvature corrections given in terms of the Gauss-Bonnet combination and weighted (in four space-time dimensions) by inflaton-dependent couplings [36,37,38]. Similar modifications of the evolution of the tensor modes occurs in the case of Einstein-aether models [39,40,41] or in the case of compact extra-dimensions [45,46]. From the physical viewpoint a common aspect of different parametrizations is that the gravitational waves may acquire an effective index of refraction [22,23], as suggested long ago [42,43] without any reference to the inflationary dynamics (see also [44]). Even if the geometry undergoes a stage of conventional accelerated expansion the intermediate slope of the spectral energy density increases depending on the evolution of the refractive index [22,23]. The different contributions to the evolution of the tensor modes of the geometry in the case of conformally flat background geometries are summarized as follows: In Eq. (2.1) the terms that would break parity and that are associated with quadratic combinations involving either dual Riemann tensor or the dual Weyl tensor have been neglected; both terms would appear in the effective action [32] (see also [47]) and might in principle polarize the relic gravitons. Equation (2.1) contains three undetermined functions A(τ ), B(τ ) and B c (τ ). The coefficient B(τ ) refers to the expanding dimensions while the presence of B c (τ ) is related to a possible mass term that arises from the internal (compact) dimensions 7 . In the case of a toroidal compactification (which is the one discussed here) γ A B = δ A B [45,46]. There are two equivalent ways in which Eq. (2.1) can be phrased in terms of an appropriate refractive index. The first possibility is to factor A(τ ): where n(τ ) and n c (τ ) denote the refractive indices associated, respectively, with the expanding and with the compact dimensions: The explicit form of the action (2.3) simplifies by changing the time parametrization and by rescaling the background dependence: If we again introduce the η-time defined as n(η) dη = dτ , Eq. (2.6) becomes where c(η) = B(η) n(η

Quantization and spectra in the η-time parametrization
Since the two approaches mentioned above are equivalent we can introduce, as usual, the normal modes of the system µ i j ( x, η) = c(η) h i j ( x, η) so that the action takes the form: where F denotes the rate of variation of c(η): In Eq. (2.9) the overdot and the prime denote, respectively, a derivation with respect to η and with respect to τ . Three complementary time parametrizations are relevant for the present analysis and their features can be summarized, in short, as follows 8 . The η-time parametrization and the conformal time are related as n(η) dη = dτ (see also Eqs. (2.2)-(2.3)); the dictionary between the two is: (2.10) The variation of the background geometry is typically expressed in terms of the cosmic time coordinate t that is related to τ as a(τ )dτ = dt and, as usual, the connection between the rates of variation of the background is: (2.11) 8 Since the usual Hubble rate H = ∂ta/a is defined in the cosmic time parametrization we remind that the overdot often denotes the derivative with respect to the cosmic time coordinate t but, in the present context, the overdot will be reserved for the derivation with respect to η-time, as suggested in Eq. (2.9).
The rate of variation of the refractive index in units of the Hubble rate and the rate of variation of H itself are then defined as: where is the usual slow-roll parameter. Since the phase velocity coincides with the group velocity, the refractive index must increase during inflation (i.e. α ≥ 0) to prevent a superluminal propagation of the signal; there are, in practice, two relevant physical situations depending on the value of α, i.e. α < 1 and α = O(1): in the first case α and are of the same order while in the second case α . In what follows α is kept generic however, as we shall see, the tensor spectral index is determined by the competition of α and and the physical range corresponds to α < 1. We finally remark, as already mentioned, 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 enforced. After these necessary specifications w can define the canonical momenta from Eq. (2.8) so that the canonical Hamiltonian becomes: (2.14) From Eq. (2.14) the Hamilton's equations are: where the Laplacian associated with the internal dimensions has been denoted by We can now quantize the system in the standard manner but, for the sake of accuracy, we repeat here the main steps. It is first useful to express the quantum field operators µ i j ( x, η) and π m n ( x, η) directly in Fourier space: The field operators in Fourier space can then be expanded in the basis of the two (linear) tensor polarizations 9 : i j (q) µ λ (q, η), π m n ( p, η) = λ=⊕, ⊗ e (λ) m n (p) π λ (p, η).
(2.17) 9 The explicit form of the two linear polarizations is given by e The field operators of Eq. (2.17) can be directly written in terms of the creation and annihilation operators obeying [ a q, λ , a † p, λ ] = δ (3) ( q − p) so that, ultimately, µ i j ( q, η) and π m n ( p, η) are: It can be directly checked that the commutation relations between µ i j ( q, η) and π m n ( p, η) are given by: where S i j m n (q) and p i j = (δ i j −q iqj ). Note that Eq. (2.20) holds provided the Wronskian normalization condition is verified: It is finally practical to deduce the two-point functions ofĥ ij ( k, η) and of ∂ ηĥij ( k, η) in Fourier space: where the power spectra P T (k, τ ) and Q T (k, τ ) are defined, respectively, by: All in all, putting everything together, we have that the field operators are expressed as: where, according to Eq. (2.15), the evolution of the mode functions f k,λ and g k,λ obeys: Thanks to the coupling among the scalar and the tensor modes the gravity wave evolution equation get what looks like a massive contribution [45,46]. In terms of the eigenstates of the Laplace operators appearing in Eq. (2.15) we have that ∇ 2 µ i j = −k 2 µ i j and ∇ 2 µ i j = −q 2 µ i j . While k denotes the external momentum, q is the momentum associated with the extra-dimensions which can be viewed as a massive contribution as it can be appreciated by decoupling Eqs. (2.28)-(2.29): the polarization index has been omitted since the result of Eq. (2.33) holds both for ⊕ and for ⊗.

Different physical limits
Depending on the values of c(η), n(η) and r(η), the action of Eq. (2.8) describes a number of relevant situations that are however physically different. In the standard limit the refractive index is absent and the internal dimensions disappear: In the limit (2.31) the rescaled normal mode becomes µ i j = a(τ )h i j and a(τ ) is the scale factor appearing in the four-dimensional line element where η µν is the Minkowski metric. According to Eq. Eq. (2.8) coincides with the original Ford-Parker action [3]; furthermore, as stressed in Eq. (2.31), the η-time and the conformal time coordinates coincide. In the case of Eq. (2.31) the only possibility of getting a flat spectrum at high-frequency is to modify the standard inflationary dynamics by considering, for instance, the possibility of bouncing backgrounds. In this case Eq. (2.30) becomes: Let us now suppose that n → 1 in the presence of d internal dimensions characterized by the scale factor b 2 (τ ) so that the line element is given by: where A, B = 1, . . . , d runs over the d internal dimensions and the dimensionality of the space-time is D = 4 + d. As already suggested, we mainly consider the case γ A B = δ A B . In the case of Eq. (2.34) the internal volume is b d/2 and we then have that the various parameters A scenario of dimensional decoupling based on Eq. (2.36) has been discussed in Refs. [45,46] and the considerations reported here can be easily extended to that situation. We finally consider the framework that is more realistic, at least for the present ends: (2.37) In the limit (2.37) the evolution of the mode functions (see Eq. (2.30)) becomes: In the present discussion we also assume that the refractive index is exactly 1 at the present time. It is useful to remark that the evolution of the refractive index is specified unambiguously by assigning n(a). Even though the phase velocity of the relic gravitons is not required to be sub-luminal we consider here the situation where n(a) ≥ 1. When n(a) changes appreciably during inflation and it goes to 1 in the standard decelerated stage of expansion 10 : Equation (2.39) defines, in practice, three successive physical regimes. For a a 1 the refractive index goes to 1 and the standard situation is recovered depending on the value of γ ≥ 1 which controls the sharpness of the transition. When a * < a < a 1 the refractive index is practically constant but still larger than 1, i.e. n(a) n * > 1. Finally for a < a * we have the truly refractive stage where n(a) n * (a/a * ) α .

Chirp amplitude, spectral amplitude and spectral energy density
When discussing the relic gravitons at the present time the observational collaborations typically assume that the space-time is flat and that the frequency of the gravitons is always larger than the rate of variation of the geometry which means, in terms of the previous notations, that kη Furthermore, since the scale factor is normalized to 1, the comoving and the physical frequencies are (today) concident. To introduce the spectral amplitude the simplest approach is to expand the tensor amplitude as where ν = k/(2π) is the comoving frequency and d k = d cos ϑ dϕ. As before h i j (ν, k) can be expanded in the basis of the linear polarisations ⊕ and ⊗: The spectral amplitude S h (|ν|) is defined as the expectation value of the tensor amplitudes expressed as a function of ν and k: where A c is an overall constant that parametrizes the different choices currently adopted by different authors 11 . Consistently with the previous notations, the angular delta function appearing in Eq.
If we now compute the expectation value of two tensor amplitudes with different indices we obtain where S i j m ( k) has been already introduced in Eq. (2.21). The expectation value of the tensor amplitudes at equal time becomes: From the direct comparison of Eq. (2.46) with the analog expression computed in terms of P T (ν) it follows that the relations between the chirp amplitude, the spectral amplitude and the power spectrum is: The specific values of A c can be used to rationalize the obtained expressions. The LIGO/Virgo collaboration is normally setting A c = 1/(16 π) so that Eq. (2.47) becomes The PTA collaborations express their results in terms of the chirp amplitude h 2 c (ν). In this respect we just note that, up to a numerical factor, the square of the chirp amplitude coincides with the power spectrum so that its relation with the spectral energy density may be easily determined and it is: 3 Normalization in the audio band and in the nHz domain

Wide-band interferometers and the audio band
The observations of wide-band detectors led through the years to a number direct upper limits on the backgrounds of relic gravitons for a frequency interval ranging between few Hz and 10 kHz [24,25,26,27,28,29,30,31]. The most relevant upper limits are summarized in Tab. 1 and they depend on the spectral slope of the signal. For the purposes of Tab. 1 the spectral energy density has been parametrized with a power-law slope of the type: Ref. [31] where ν ref is a (conventional) frequency while Ω(σ) is the constant amplitude that differs 12 depending on the slope σ. The results of Tab. 1 show that the constant amplitudes associated with the various σ are constrained at different levels. The scale-invariant limit represents the simplest signal grossly compatible with the concordance paradigm; conversely, if σ = 3 in Eq. (3.1), the factor Ω gw (ν) 2 /ν 6 is practically constant and, in this case, the estimate of the integral appearing in the signal-to-noise ratio gets simpler (see e.g. [49,50,51,52]). In Refs. [24,25,26,27,28,29] the LIGO/Virgo collaboration presented the upper limits for Ω(0) and Ω(3); a third case has been subsequently analyzed and it corresponds to σ = 2/3 [30,31]. In the case of an exactly scale-invariant spectrum the constraints obtained in Refs. [29,30] imply that Ω(0) < 6 × 10 −8 . If this value is compared with the analog result obtained in Ref. [27], the upper limit is O(100) times more constraining for the same slope and for the same frequency band. For σ = 2/3 the constraints of Ref. [29] imply Ω 2/3 < 4.8 × 10 −8 with 95% confidence within the 20-95 Hz frequency band with ν ref = 25 Hz. The slope σ = 2/3 may actually parametrize a potentially interesting foreground for the relic gravitons 13 . The most constraining limit to date in the case σ = 2/3 has been obtained in Ref. [31] by the Kagra-Ligo-Virgo collaboration and it requires that Ω(2/3) < 3.4 × 10 −9 . The most constraining set of bounds appearing in Tab. 1 corresponds to the one obtained by the LIGO, Virgo and KAGRA collaborations [31]. In the case of a flat spectral energy density the bound reads 14 : Ω gw (ν L ) < 5.8 × 10 −9 , 20 Hz < ν L < 76.6 Hz, (3.2) where ν L denotes the LIGO-Virgo-KAGRA frequency; we shall commonly refer to this limit as the LVK bound. Even if in Eq. (3.2) we just quoted the most constraining limit, the LIGO-Virgo-KAGRA collaboration actually reports a threefold bound for three different values of σ. When the value of σ increases the bound becomes more restrictive once the reference frequency is kept fixed. The three results are unified in the following interpolating formula that will be used to impose the LVK bound at high-frequency.  Ref. [21] tions [18,19,20,21]. Before going through the details we can already see that the determinations of the spectral energy density seem to be, naively, at the same level of the LVK bound but in a different frequency domain. The two results of Ref. [20] differ because of the different values of β and the same comment also holds in the case of Ref. [21]. The various PTA collaborations [18,19,20,21] express the chirp amplitude at a pivot frequency ν ref = 31.68 nHz corresponding to yr −1 :

Pulsar timing arrays and the nHz band
The pivotal model analyzed so far assumes β = −2/3 and this is the case preferentially reported in Tab. 2. The EPTA [20] and IPTA [21] collaborations also consider a more general class of scenarios where Q and β vary simultaneously. For instance the EPTA finds that the most favoured model to be the common uncorrelated red noise described by Q = 5.13 +4.20 −2.73 × 10 −15 with γ = 3.78 +0.69 −0.59 where we recall that, within the present notations, β = (3 − γ)/2. If the spectral index is instead fixed as γ = 13/3 (i.e. β = −2/3) Ref. [20] suggests Q = 2.95 +0.89 −0.72 ×10 −15 . We not that γ will is not employed hereunder and are it is only mentioned for the sake of accuracy since some of the PTA collaborations introduce this notation which is a bit contrived from the viewpoint of the present discussion.
For the different estimates of Refs. [18,19,20,21] the value of Q is always O(10 −15 ) it is therefore useful to express Q = q 0 × 10 −15 . Using this notation and Eq. (3.4) the spectral energy density in the nHz band can be finally expressed as 15 : (3.5) The potential signal of the PTA recently reported evidence of a potential signal in the nHz band.
Using the spectral energy density in critical units as a pivotal variable the features of this purported signal would imply, in the present notations, that:

Audio band and nHz band: which is the most restrictive?
To avoid potential confusions it is relevant to compare the limits in the audio band and in the nHz to understand which are the most restrictive. Before discussing this comparison it is equally important to stress that the two classes of measurements are qualitatively different: while the limits from the audio band are specific to the case of relic gravitons [30,31], the property of a PTA is that the signal from relic gravitons should be correlated across the baselines while that from the other noise will not. Since these correlation have not been observed so far, the interpretation suggested in by the PTA is still preliminary and the measurements of Refs. [18,19,20,21] might not have anything to do with relic gravitons: the correlation signature of an isotropic gravitational wave background follows the so-called Hellings and Downs curve which depends on the angle between a pair of Earth-pulsars baselines. As already mentioned this correlation has not been observed yet by admission of the various experimental collaborations [18,19,20,21]. If we assume that the PTA are indeed related with relic gravitons the constraint of Eqs.
where the limit has been referred to a fiducial value of h 0 that follows from the CMB data [5,6,7]. By comparing Eqs. (3.6) and (3.7) it seems that the former is superficially more constraining than the latter: by choosing q 0 = 1 we would have that h 2 0 Ω gw (ν) < 10 −8.86 . However q 0 is not 1; on the contrary if we take the average of the four measurements presented so far (see Tab. 2 in the case β = −2/3) we obtain and averaged value given by q 0 = 2.467 which implies (3.8) Note that Eq. (3.7) is always more constraining than Eq. (3.8) even if we choose the smallest value of q 0 which is the one associated with the NANOgrav estimate [18]: if q 0 = 1.92 we get from Eq.

Big-bang nucleosynthesis limits
While the PTA measurements constrain the spectral energy density at intermediate frequencies, the bounds coming from big-bang nucleosynthesis [55,56,57] imply a constraint on the integral where Ω γ0 is the (present) critical fraction of CMB photons. The limit of Eq. (3.9) sets an indirect constraint on the extra-relativistic species possibly present at the time of nucleosynthesis. Since Eq. (3.9) is relevant in the context of neutrino physics, 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. (3.9) is thus between 10 −6 and 10 −5 . It is relevant to point out, as we shall see, that the upper limit of integration (labeled by ν max ) depends on the specific post-inflationary evolutions 16 . Conversely, the lower limit of integration in Eq. (3.9) is given by the frequency corresponding to the Hubble rate at the nucleosynthesis epoch: Hz 0.01 nHz, (3.10) where N ef f 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. We finally remark that the bound of Eq. (3.9) could be relaxed if the nucleosynthesis takes place in the presence of matter-antimatter domains [56]. This possibility will not be specifically considered hereunder and we shall instead enforce the bound of Eqs. (3.9)-(3.10) in its conservative version. As we shall see when the quasi-flat spectrum is normalized in the audio band the limit (3.9) is always satisfied.

Conditions for flat spectra at high-frequency
In the conventional situation the spectral energy density for typical frequencies larger than the nHz is always smaller than 10 −15 . If all the sources of late-time suppression are taken into account we approximately have h 2 0 Ω gw (ν) = O(10 −16.5 ). This conclusion can be however evaded in, at least, two complementary situations that are simultaneously illustrated in Fig. 1 where the common logarithm of | F =ḃ/b | is reported. The blobs appearing in the cartoon represent the transition regimes between Figure 1: We schematically illustrate the evolution of F when the conventional inflationary phase is preceded by a further evolutionary stage.
the different stages and they are immaterial in the approach based on WKB approximation that is adopted in the present section. Figure 1 illustrates two of the three limits discussed in Eqs. (2.31)-(2.35) and (2.37). In the case n → 1 the η-time coincides with the conformal time coordinate and F = H = a H. Furthermore for a > a * the refractive index is not dynamical, and, in this limit, c = a. The most conservative case is the one associated with Eq. (2.37) where the refractive index evolves and we shall consider, in particular, the situation where n(a) increases and then gets back to 1 for c > c * (see Eq. (2.39) and discussion therein). For k < k * the amplitude and the slope of Ω gw (k, τ ) depend on the profile of F for c < c * ; this means that for k eq < k < k * , Ω gw (k, τ ) increases provided the rate of variation of n is sufficiently large in Hubble units. In the language of Eq. (2.12) this implies α > 0 even if, as we shall see, α cannot exceed 1 if the results of the wide band detectors are used to set the high-frequency normalization of the spectral energy density. According to Fig. 1 the spectral range k * < k < k max is determined by the evolution of F for c * < c < c 1 where F ∝ a −1 . This scaling actually corresponds to the conventional situation since, in this interval, F ∝ H aH so that the Hubble rate is roughly constant. In the minimal situation where the refractive index flattens out for c > c * the spectral energy density is quasi-flat (or slightly decreasing) simply because, in this range, we get back to the conventional case where Ω gw (k, τ ) is determined by the wavelengths crossing the Hubble radius during the inflationary stage and reentering when the background is dominated by radiation. Indeed for c > c 1 the evolution is completely standard and F = H ∝ a −1 . Finally for c > a eq we have that F = H ∝ a −1/2 . In Fig. 1 we did not include the stage dominated by the dark energy simply because its contribution to the spectral energy density is immaterial for the considerations of this section. A more detailed account can be however found in section 5 where the late-time suppression of the spectral energy density will be more specifically investigated. For the normalization in the high-frequency regime it is useful to have a fairly general expression of the spectral energy density that can be deduced rather simply within the WKB approximation; a key role is played, in this context, by the structure of the turning points η ex and η re . As illustrated in Fig. 1, a generic wavenumber k < k max crosses |F(η)| twice in c ex and c re and these two values correspond to the moments where a given wavelength exists and reenters the effective Hubble radius |F(η)| −1 . From a technical viewpoint η ex and η re are defined as the turning point at which the solution to Eq. (2.38) changes its analytic form. In particular the exit corresponds to: Equation (4.1) actually defines a regular turning point and it tells that, at η ex , kη ex 1. It can happen, however, that the turning point is singular and this happens, in particular, whenc → 0 in the vicinity of the turning point. In the problem at hand we have thatc ex = 0 however, if the reentry takes place during a radiation-dominated stage of expansion, we typically havec re → 0 since, in a radiation stage, η = τ and a = 0. In this case, as we shall see in a moment, we must recall that the condition kη re 1 is not verified and it must be replaced by kη re 1. In what follows we shall first discuss the expression of the spectral energy density for the different ranges of wavelengths and then analyze the way the high-frequency normalization must be implemented.

Wavelengths larger than the effective horizon
To deduce the general evolution of the mode function for large wavelengths it is practical to transform the relevant differential equations in a set of integral equations. The evolution of the mode function of Eq. (2.38) is equivalent to 17 an integral equation whose initial conditions are assigned at the reference time η ex : where, by definition, c ex = c(η ex ). As argued in Eq. (4.1) the first turning point is always regular so that kη ex 1. Neglecting then the terms O(k 2 η 2 ) the lowest order solution of Eq. (4.2) is: Equation (4.3) determine the approximate form of the power spectrum for wavelengths larger than the 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, the explicit expression of the tensor power spectrum follows from Eq. (2.25) by recalling that, for η < −η * , The expression of c(η) follows, in its turn, from the definition of the η-time and from Eq. (2.39) evaluated in a conventional inflationary background with slow-roll rate given by . In the limit | k η | 1 the tensor power spectrum is in fact constant and it is approximately given by: where A(α, ) and n T (α, ) are defined, respectively, by: Within the approximation scheme leading to Eq. (4.3) we have that in Eq. (4.5) |A| = √ 2k |f k (η ex )| = 1. In a more general perspective 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 |c/c| the correctly normalized solutions of Eq. (2.33) are f k (η) = e ± ikη / √ 2k. However as soon as η ex is approached the amplitude gets slightly modified and the exact solution of Eq. (2.33) is given in terms of Hankel functions [53,54]. The power spectrum of Eq. (4.5) can be further simplified by observing that: 1 (4.7) If we now use Eqs. (4.6)-(4.7) into Eqs. (4.5) we obtain the following form of the power spectrum for k < a * H * (4.8) For the modes k > a * H * the spectrum (4.8) is modified since now the evolution of the mode functions is given by: Recalling that a /a = a 2 H 2 (2 − ) and that a H = −1/[(1 − )τ ] Eq. (4.9) becomes If the refractive phase terminates before the end of inflation the power spectrum inherits a further branch for a * H * < k ≤ a 1 H 1 :

Wavelengths shorter than the effective horizon
All the cosmic gravitons measured at the present time are inside the Hubble radius and provided the reentry occurs whenc re = 0 we have that, approximately, kη re = O(1). Conversely ifb re → 0 in the vicinity of the turning point, then kη re 1. For η ≥ η re the solution of Eq. (2.33) is where f re (η) are the mode functions inside the effective horizon (i.e. e −ikη / √ 2 k in the crudest approximation). The coefficients C ± (k, η ex , η re ) are . (4.13) If the reentry takes place, as we are considering here, when the refractive index is not dynamical, Eq. (4.13) can be simplified even further since b(η) → a(τ ) and F → H: where −τ 1 marks, as before, the end of the inflationary stage. Because c(η) always increases, in Eqs. (4.14)-(4.15)the terms proportional to |c ex /c re | can be neglected in comparison with |b re /b ex |. Since C ± (k) are both complex but subjected to the condition |C + (k, η ex , η re )| 2 − |C − (k, η ex , η re )| 2 = 1 it is sufficient to estimate the approximate form of |C − (k, η ex , η re )| 2 : Equation (4.16) 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: From the ratio between Eqs. (4.17)-(4.18) the standard relation between the power spectrum and the spectral energy density is recovered and it is generally valid when the relevant wavelengths are shorter than the Hubble radius at a given epoch. Inside the Hubble radius we can evaluate indifferently either the power spectrum or the spectral energy distribution. It is finally useful to estimate more explicitly k * and k max . Since k * = 1/η * we also have that As usual in Eq. (4.20) N * = ln (a * /a i ) denotes the number of e-folds during the refractive stage N t = ln (a 1 /a i ) is the total number of e-folds. Recalling that ν max = k max /(2π) we then have, in explicit terms: MHz, (4.21) where A R is the amplitude of the scalar power spectrum at the pivot scale k p = 0.002 Mpc −1 and Ω R0 is the total fraction of relativistic species at the present time in the concordance paradigm. Thanks to Eq. (4.20) we also have, by definition, that We note that in Eq. (4.21) we introduced the slow-roll parameter and not r T since we did not assume the consistency relations so that r T might be smaller (or even much smaller) than 0.06.

Flat spectra at high frequencies
The spectral energy density for typical wavenumbers a * H * < k < a 1 H 1 is quasi-flat and to show this point in general terms it is useful to go back to Eqs. (4.16) and (4.18): (4.23) As already mentioned twice, in Eq. (4.23) we have two complementary possibilities depending on the nature of the turning point 18 . When the wavelength renters the Hubble radius during radiation the correct limit is kτ re 1 and the spectral energy density at high-frequency can then be estimated as: For the slice of wavenumbers a * H * < k < a 1 H 1 the exit takes place after the refractive phase where a ex = n −1/(1− ) * |k τ 1 | 1/(1− ) . From Eq. (4.24) we then obtain: where, as usual, ∆N = (N t − N * ). The same reasoning for lower frequencies leads to: (4.28) 18 If the reentry takes place during the radiation phase we have that, in the vicinity of τre a → 0 so that kτre 1 in Eq. (4.23). Conversely the modes exit during the inflationary phase when ω 2 τex 1 (i.e. k 2 τ 2 ex n 2 * ).

Equations (4.25)-(4.26) and (4.27)-(4.28)
give Ω gw (k, τ ) in the three spectral regions defined in Fig.  1. In summary the lowest frequency region involves the modes exiting the Hubble radius during the refractive phase and reentering after equality (i.e. k < a eq H eq ). The intermediate region concerns the modes exiting the effective horizon during the refractive phase and reentering during radiation (i.e. a eq H eq < k < a * H * ). The highest frequency domain encompasses the modes that exit the Hubble radius after the end of the refractive phase and reenter during radiation (i.e. a * H * < k < a 1 H 1 ). As we shall see the relevant physical regime is the one where the rate of variation of n is larger than ; in this limit we can expand the spectral index for 1: where the second equality follows in the limit 1. The spectra obtained so far in the case of a dynamical refractive index are fully compatible with a conventional inflationary stage and this is ultimately the reason for the flatness of Ω gw (k, τ ) at high-frequencies. It is however interesting to mention that the parametrization of Eqs. (4.25)-(4.26) and (4.27)-(4.28) holds also, with the appropriate differences, when the spectrum at intermediate frequencies is not dictated by the evolution of the refractive index. Two cases are particularly important in the light of the current data: the bounces of the scale factor and the curvature bounces. For a bounce of the scale factor the scale factor first contacts in an accelerated manner 19 and then undergoes a stage of decelerated expansion (i.e.ȧ > 0,ä < 0). In the case of an accelerated contraction the scale factor can be parametrized as a power-law a(t) (−t/t 1 ) δ with 0 < δ < 1. When δ < 0 we have instead an accelerated expansion with growing curvature (i.e.Ḣ > 0). The spectral energy density of the relic gravitons produced in this kind of scenarios can be estimated as in the previous case with the relevant similarity that the intermediate spectral index is also blue. In the case of accelerated contraction we have In the case of accelerated expansion with δ < 0 we have instead For ν > ν * the slope s controlled by m T even if, in this case, the quasi-flat slope is not related to a slow-roll dynamics as in the conventional inflationary case. We can therefore conclude that Ω * is generally given by the product of two separate contributions where Ω late denotes the late-time contribution that only depends on the parameters of the concordance scenario; in Eq. (4.32) Ω early is instead the early contribution that is generally model-dependent. In particular we have that, in the present case, Ω * = Ω * (α, , N * , N t ).

High-frequency normalization of the spectral energy density
The high-frequency normalization can be studied accurately by considering all the late-time sources of suppression but this analysis is postponed to the following section since, in what follows, the attention is focused on the general logic that can be more directly appreciated from the analytic results deduced above. The first observation is that in the conventional situation (H 1 /M P ) appearing in Eq. (4.25) is fixed from the amplitude of the scalar power spectrum A R and from the slow-roll parameter: (4.33) When the normalization is set at low-frequencies is related to the r T whose upper limits fix the spectral energy density in the aHz range. In the present case, however, the situation is different and there are, in purely abstract terms, two complementary possibilities: • the first logical possibility is to fix the normalization by requiring that the spectral energy density matches the value measured by the PTA at the typical frequency ν = ν P = O(30) nHz; • the second possibility is instead to normalize the potential signal in the audio band by enforcing the LVK bound at the frequency ν = ν L = O(60) Hz.
Between these two possibilities the latter is more plausible than the former for the simple reason that, generally speaking, ν * < ν L so that the quasi-flat branch of the spectrum falls in the audio band.
On the other hand it is plausible to have that ν * is either larger or smaller than ν P . This point is illustrated in Fig. 2 where the darker area defines the region where 100 aHz < ν * < ν P . The various labels appearing on the curves define the common logarithm of ν * expressed in Hz. We have that, approximately, ν P = O(10 −7.5 ) Hz. Figure 2 explains why the high-frequency normalization has ben fixed by requiring According to Eqs. (3.2)-(3.3) the absolute upper limit for the quasi-flat spectral energy density should be O(5.8) × 10 −9 and, for the sake of simplicity, we are going to require, in a conservative perspective 20 , that Ω gw (ν L ) = 5.8 × 10 −9 . To illustrate the procedure (and to avoid the possible complications that are addressed in the next section) we assume, in the notations of Eq. (4.32), that , 100 aHz < ν * < ν L . The two plots differ because of the total number of e-folds but are otherwise qualitatively similar. The darker region in the central part of both plots corresponds to 100 aHzν * < ν P . The various curves are the contours where the values of ν * remain the same and in the labels we report the common logarithm of ν * expressed in Hz.
As we mentioned in Eq. (3.2) the value of ν L ranges approximately between 20 and 76 Hz. This range is related to the frequency region which is more sensitive to the backgrounds of relic gravitons; for the sake of concreteness we therefore posit ν L = 60 Hz. To pass from Eq. (4.37) we first employed Eq. (4.34) and also noted that ν L /ν * = (ν L /ν max )(ν max /ν * ) since the ratio (ν max /ν * ) can be directly estimated from Eq. (4.22). In practice all the terms at the left of Eq. (4.37) contain the parameters of the model while the quantities ate the right-hand side are either directly measured or can be determined as late-time parameters of the concordance paradigm. For a swift estimate of the parameters the right-hand side of Eq. (4.37) can be evaluated in the limit m T ( ) 1. If we now take the logarithm of both sides of Eq. In the limit m T ( ) 1 Eq. (4.38) roughly implies 3 α N * 15.62. This mans that α and N * are inversely proportional: a longer refractive phase imposes a smaller α and vice-versa. Having determined the normalization according to Eq. (4.38) we can directly write spectral energy density as , ν eq < ν < ν * , where we stress that B(ν L ) = 10 −8.61 has been fixed by normalising the spectral energy density to the largest value compatible with the LVK bound. Since we already saw from Fig. 2 that the values of ν * can be either larger or smaller than ν P it is natural to parametrize ν * in terms of ν P by setting ν * = f 0 ν P where f 0 can be either smaller or larger than 1. The discussion developed so far assumes PTA LVK a conventional inflationary stage where H 1 /M P is determined from the scalar power spectrum and, in this case, Eq. (4.33) holds. The same strategy adopted above can however be applied also in the complementary situation where the universe bounces across a typical curvature scale of the order of H 1 . In this case the spectral energy density for ν > ν * is given by: The main difference is that in the case of Eq. (4.40) the scale H 1 cannot be estimated as in Eq. (4.33). On the contrary it is exactly the LVK bound that sets the scale of H 1 : (4.41) As already mentioned the results of Eqs. (4.38) and (4.40)-(4.41) do not include the late-time damping that has the effect of shifting a bit the above estimates as we shall see in the following section. The logic of this section has been to deduce analytically the high-frequency plateau and to normalise its value to the LVK bound.

Phenomenological considerations
The previous section illustrates how the PTA data and the constraints provided by the LVK collaboration set the normalization of the spectral energy density without any reference to the low-frequency data of the aHz region. The general expression of the high-frequency amplitude consists of two contributions that are formally distinguished in Eq. (4.32). In the previous section, for the sake of simplicity, only the leading contribution to the late-time suppression has been included while some of the phenomena that may further suppress the high-frequency plateau have been neglected. Even though the considerations presented hereunder are relevant, from a quantitative viewpoint, for an improved determination of the theoretical template, the general logic pursued so far remains unaltered.
If the late-time contribution (i.e. Ω late in Eq. (4.32)) is reduced, the early-time contribution may become comparatively larger whenever the high-frequency normalization is imposed (see also Eqs. (4.35)-(4.36)). In this respect it is useful to answer three separate questions. The first one concerns the overall magnitude of the suppression in the high-frequency domain where the normalization is set. The second issue calls for a distinction between the damping effects that depend on the frequency and the ones that are instead frequency-independent. It is finally useful to gauge the relative error on the spectral energy density when these effects are simply neglected, as tentatively assumed in the previous section.

Impact of the neutrino free streaming
The effect of neutrino free streaming on the spectral energy density of the relic gravitons has been first pointed out by in Ref. [58] (see also [59,60,61,62]) where the various authors first argued and then confirmed that the correction to the spectral energy density is mildly dependent on the frequency and it is overall of the order of the 10 %. Since the effect of neutrino free-streaming is fully operational for ν < ν bbn it is comparatively less relevant for the high-frequency normalization. The logic is that, after neutrino decoupling, the neutrinos free stream and the effective energy-momentum tensor acquires, to first-order in the amplitude of the plasma fluctuations, an anisotropic stress. In the present case the spectral energy density at intermediate frequencies is not quasi-flat (as in the standard case) but it increases as a function of the frequency but still the magnitude of the effect depends on R ν , i.e. the neutrino fraction in the radiation plasma: In Eq. (5.1) 3 counts the degrees of freedom associated with the massless neutrino families, (7/8) arises because neutrinos follow the Fermi-Dirac statistics; the factor (4/11) 4/3 stems from the relative reduction of the neutrino (kinetic) temperature (in comparison with the photon temperature) after weak interactions fall out of thermal equilibrium. Assuming that the only collisionless species in the thermal history of the Universe are the neutrinos and recalling Eq. (5.1), the amount of suppression can be parametrized by the function This suppression is effective for relatively small frequencies which are larger than ν eq and smaller than ν bbn . In what follows we shall stick to the case of the ΛCDM paradigm [5,6,7] but more complicated examples do not change the nature of the effect. In the context of the concordance paradigm, the first bracket at the right-hand side of Eq. (5.3) is evaluated during the radiation stage and it does not introduce any suppression as long as a 4 H 2 is strictly constant. This conclusion is, however, slightly inaccurate because of the evolution of the relativistic species. In local thermal equilibrium the total entropy and energy densities of a relativistic plasma can be written, respectively, as s t = 2π 2 g s (T )T 3 /(45) and and as ρ t = π 2 g ρ (T )T 4 /(30) where T is the common temperature of all the species of the plasma. At the plasma cools down the effective number of spin degrees of freedom appearing in the total entropy density and in the total energy density (i.e. g ρ (T ) and g s (T )) eventually decrease in a computable manner depending on the microscopic description of the plasma. If all the species of the plasma are in local thermal equilibrium at the same temperature we have that g s and g ρ coincide and this is what happens in the standard model for temperatures larger than the top quark mass where g s = g ρ = 106.75. The evolution of the relativistic species suggests that, during the radiation stage,

Impact of the decoupling of relativistic species
In principle if a given mode k reenters the Hubble radius at a temperature T k the spectral energy density of the relic gravitons is (kinematically) suppressed by a factor which can be estimated as [g ρ (T k )/g ρ0 ][g s (T k )/g s0 ] −4/3 [63,64,65]; at the present time g ρ0 = 3.36 and g s0 = 3.90. In general terms the effect parametrized by Eq. (5.4) will cause a frequency-dependent suppression, i.e. a further modulation of the spectral energy density Ω gw (k, τ 0 ). The maximal suppression one can expect follows by inserting into Eq. (5.4) the largest g s and g ρ . So, in the case of the minimal standard model this would imply that the suppression (on Ω gw (k, τ 0 )) will be of the order of O(0.38).
In popular supersymmetric extension of the minimal standard model g ρ and g s can be as large as O(230) reducing the previous estimate to O(0.29). These considerations demonstrate, a posteriori, the overall correctness of the assumptions behind Eq. (5.4) and its descendants.

Impact of the dominance of dark energy
The largest wavelengths of the spectrum crossed Hubble radius about 65 e-folds before the end of inflation and reentered when dark energy was already dominant. In the vanilla ΛCDM scenario the role of the dark energy is played by the cosmological constant and since he evolution of |a H| is monotonically decreasing for a eq < a < a Λ , all the wavelengths that left the Hubble radius during inflation will remain within the Hubble radius after they reenter either during radiation or during matter dominance. If we also consider the evolution for a > a Λ , the behaviour of |a H| is, overall, non-monotonic: a bunch wavelengths that reentered after equality will again exit after dark energy dominates at the redshift We can now consider the mode k Λ = a Λ H Λ , i.e. the mode approximately reentering the Hubble radius at τ Λ . Since H is approximately constant for a > a Λ we have that H Λ ≡ H 0 where H 0 is the present value of the Hubble rate. A typical wavenumber that is today of Hubble size (i.e. k 0 = a 0 H 0 ) is larger than k Λ even if, according to a superficial intuition, it should be smaller 21 and this is because of the non-monotonic evolution of aH. This range of wavelengths is currently inside the Hubble radius and their power spectrum is given by: The spectral energy density when the relevant wavelengths are inside the Hubble radius can be obtained from Eq. (5.6) and the result is: where we used that, by definition, k Λ /k 0 = (Ω M 0 /Ω Λ ) 1/3 . If the dominance of dark energy is completely neglected, Eq. (5.7) can be written in the same form where, however, the term (Ω M 0 /Ω Λ ) 2 is absent. We conclude that, in this branch of the spectrum, the dominance of dark energy suppresses the spectrum by a factor (Ω M 0 /Ω Λ ) 2 = (0.44) 2 = 0.193. The wavenumbers falling instead in the interval k Λ < k < k 0 correspond to wavelengths reentering the Hubble radius during the matter-dominated epoch and exiting again the Hubble radius when dark energy is already dominant. These wavelengths are currently outside the Hubble radius. The same logic leading to Eq. (5.6) implies that the power spectrum in the interval k Λ < k < k 0 is given by P T (k, τ re ) (a re /a ex ) 2 or, more precisely, Equation (5.8) gives the power spectrum for typical wavelengths that are larger than the Hubble radius at the present time. The wavelengths belonging to the interval k Λ < k < k 0 reenter the Hubble radius but then exit again and while they are larger than the Hubble radius their power spectrum remains unchanged. This is why, in Eq. (5.8), the power spectrum is not suppressed by a further factor (a re /a 0 ) 2 : these wavelengths are larger than the Hubble radius at the present time and therefore the power spectrum remains constant exactly as in the case of the scales exiting the Hubble radius at the onset of inflation. Using the identity k Λ /k 0 = (Ω M 0 /Ω Λ ) 1/3 the spectrum of Eq. (5.8) can be written as: We can finally compute the spectral energy density in this interval and we obtain Ω gw (k, τ 0 ) = P T (k, τ re ) 12 If the presence of dark energy would be neglected, Ω gw (k, τ 0 ) would scale as (k/k 0 ) −2 and the suppression going as (Ω M 0 /Ω Λ ) 2 would be totally absent. Note finally that when k = k Λ Eq. (5.9) implies that Ω gw (k Λ , τ 0 ) is only suppressed as (Ω M 0 /Ω Λ ) 4/3 .

Numerical discussion of the spectral energy density
After having assessed the main sources of damping at high-frequency we now go back to Eqs. (4.35)-(4.36) and (4.37) with the purpose of improving the quantitive evaluation of the late-time suppression.
Since the neutrino free-streaming is only effective at comparatively low-frequencies, at the right-hand side of Eq. (4.37) we should add a further contribution, namely In Fig. 4 we illustrate the normalization curves deduced in Eq. (4.37). In particular, while the full lines refer to the case when the late-time damping is neglected, the dashed curves include instead the different sources of damping previously analyzed in this section. The contribution associated with the neutrino free-streaming affects frequencies that are approximately smaller than ν bbn and, for this reason, they are less relevant for the high-frequency normalization. At intermediate frequencies the spectral energy density is instead affected and this observation is relevant for the potential signal of the PTA. This aspect is illustrated in Fig. 5 where the spectral energy density has been reported for an illustrative choice of the parameters that is now specifically discussed. The pivotal parameters that determine the spectrum are primarily α, N * and N t . When N * and N t are of the same order the transition to normalcy occurs at the end of inflation but in this case it is impossible to get a large signal in the nHz range without jeopardizing the big-bang nucleosynthesis constraint of Eqs. (3.9)-(3.10). If we ought to address the PTA measurements we must require N * < N t since, in this case, the transition to normalcy takes place 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). In Fig. 5 we choose N t = 65 and set = 0.003. Even if this value has been deduced from the upper limits on r T [5,6,7] and from the consistency we could easily consider values 0.003 and r T 0.06 since they are immaterial for the overall normalization but only for the slope m T ( ) of the high-frequency plateau. The values of α appearing in Fig. 5 are around α 2/7. In this case, the spectral index at intermediate frequencies is given by n T 2/3, up to slow-roll corrections O( ) which are negligible since < 10 −3 . This value of α actually corresponds to β = −2/3 (see, in this respect, Tab. 2 and discussions thereafter). As discussed the PTA results are often reported in terms of a chirp amplitude scaling as ν −2/3 for a typical reference frequency O(yr −1 ). The value α = 2/7 0.28 corresponds to β = −2/3 and n T 2/3. More precisely we have that α = (2 + 4 )/7 which can be approximated as α = 2/7 + O( ) for < 10 −3 . In Fig. 5 the box illustrates the PTA measurements. In Fig. 5 the spectral energy density in critical units has been normalized by using the limits from the wide-band detectors given in Tab. 1 and by also imposing, as a particular choice, that α = O(0.28) as suggested by the PTA measurements of Tab. 2. It is interesting that these two independent choices lead to a large signal in the nHz range when the variation of the refractive index occurs sufficiently early during the inflationary stage and anyway not beyond the first 20 e-folds. It is finally important to remark that, in the present context, then constraints on the integral of the spectral energy density of the relic gravitons (see Eqs. (3.9)-(3.10)) are automatically satisfied after imposing the normalization deduced from the direct limits set by the operating interferometers (see Tab. 1 and discussion thereafter).

Concluding remarks
Within the concordance paradigm the absolute normalization of the spectral energy density of the relic gravitons in critical units is assigned in the aHz region and it depends on r T (i.e. the tensor to scalar ratio) that should not exceed O(0.06), at least according to the limits set by the temperature and polarization anisotropies of the microwave background. Since the ΛCDM is a compromise between the available data and the number of ascertainable parameters, the most stringent limits on r T hold when the consistency relations between the scalar and tensor power spectra are enforced; in practice this only happens in the case of single-field inflationary models. If the spectral energy density of the cosmic gravitons is predominantly distributed for frequencies much larger than the aHz the logic leading to the low-frequency normalization is less compelling. In the nHz region the PTA recently reported a potential excess even if a bona fide signal coming from the relic gravitons should be correlated across the baselines and so far no indications along this direction have been obtained. Motivated by the improved limits in the audio band and by the current data from the pulsar timing arrays in the nHz domain, we analyzed the conditions for a quasi-flat spectrum of relic gravitons at intermediate and high-frequencies by introducing an improved physical strategy for the absolute normalization of the cosmic background of relic gravitons.
After proposing a general four-dimensional action for the discussion of relic gravitons in spatially flat backgrounds, we concentrated on the classes of scenarios where a large signal can be expected between the nHz and kHz ranges. While the most promising possibilities involve a dynamical refractive index and the bouncing dynamics, between these two cases the former is slightly more conservative than the latter insofar as it is compatible with the presence adiabatic and Gaussian initial data for the temperature and polarisation anisotropies of the microwave background. In both situations the spectral energy density increases over intermediate frequencies and then flattens out. Since the data from the wide-band detectors set the normalization of the spectral energy density at high-frequencies, a scheme based on the WKB method is preferable for a general estimate. Within this approach the early contributions can be easily distinguished from the late-time effects that are evaluated, depending on the convenience, in different approximations. The results obtained here also suggest an effective mechanism for the origin of a flat spectrum of relic gravitons with typical amplitudes that are even six or seven orders of magnitude larger than in the case of conventional inflationary models.
The results obtained here also suggest that the signal of the PTA can be explained by a background of relic gravitons of inflationary origin without conflicting with the bounds coming from big-bang nucleosynthesis which are automatically satisfied as long as the data from wide-band interferometers set the normalization of the spectral energy density in critical units between few Hz and 0.1 kHz. In the present framework the low-frequency constraints can also be imposed a posteriori as a limit on the intermediate spectral slope but they are overall less crucial and they should be applied, strictly speaking, only when the consistency conditions between scalar and tensor modes are enforced.