Relic gravitons and pulsar timing arrays: a theoretical viewpoint

During the last three years the pulsar timing arrays reported a series of repeated evidences of gravitational radiation (with stochastically distributed Fourier amplitudes) at a benchmark frequency of the order of $30$ nHz and characterized by spectral energy densities (in critical units) ranging between $10^{-8}$ and $10^{-9}$. While it is still unclear whether or not these effects are just a consequence of the pristine variation of the space-time curvature, the nature of the underlying physical processes would suggest that the spectral energy density of the relic gravitons in the nHz domain may only depend on the evolution of the comoving horizon at late, intermediate and early times. Along this systematic perspective we first consider the most conventional option, namely a post-inflationary modification of the expansion rate. Given the present constraints on the relic graviton backgrounds, we then show that such a late-time effect is unable to produce the desired hump in the nHz region. We then analyze a modified exit of the relevant wavelengths as it may happen when the gravitons inherit an effective refractive index from the interactions with the geometry. A relatively short inflationary phase leads, in this case, to an excess in the nHz region even if the observational data coming from competing experiments do not pin down exactly the same regions in the parameter space. We finally examine an early stage of increasing curvature and argue that it is not compatible with the observed spectral energy density unless the wavelengths crossing the comoving horizon at early times reenter in a decelerated stage not dominated by radiation.


Introduction
The low-frequency gravitons notoriously affect the the propagation of electromagnetic signals so that the temperature and the polarization anisotropies of the cosmic microwave background (CMB) are in fact detectors of relic gravitational waves in the aHz range 2 (see, for instance, [1]).Even though the physical nature of the effects remains practically the same, the nHz range exceeds the Cosmic Microwave Background (CMB) frequencies by roughly 10 orders of magnitude and, as suggested long ago [2,3,4], the millisecond pulsars can be employed as effective detectors of random gravitational waves for a typical domain that corresponds to the inverse of the observation time during which the pulsar timing has been monitored.More specifically, Sazhin [2] suggested that the arrival times of pulsar's pulses could be used for direct searches of gravitational radiation and shortly after Detweiler [3] made a similar point by deriving one of the first upper limits on the relic gravitational waves in the nHz range.This aspect has been more accurately reinstated in Ref. [4] where the authors also derived an important property of the isotropic backgrounds of gravitational radiation stipulating that the signal coming, in our context, from relic gravitons should be correlated across the baselines while that from others noises will not.This effect depends on the angle between a pair of Earth-pulsars baselines and it is often dubbed by saying that the correlation signature of an isotropic and random gravitational wave background follows the so-called Hellings-Downs curve [4].If the gravitational waves are not characterized by stochastically distributed Fourier amplitudes the corresponding signal does not necessarily follow the Hellings-Downs correlation.
By using this logic, since the mid 1990s a series upper limits on the spectral energy density of the relic gravitons in the nHz range has been obtained [5,6,7,8] and during the last three years the pulsar timing arrays (PTA) reported an evidence that could be attributed to isotropic backgrounds of gravitational radiation [9,10,11,12].The PTA collaborations 3 customarily assign the chirp amplitude at a reference frequency ν ref = 31.68nHz that corresponds to yr −1 : The pivotal models analyzed in Refs.[9,10,11,12] assumed β = −2/3 or, which is the same 4 , γ = 13/3.In Eq. (1.1) the value of Q inferred from the observational data should be O(10 −15 ) so that we find it natural to parametrize the amplitude at ν ref as Q = q 0 × 10 −15 .In the previous data releases the q 0 ranged between 1.92 and 5.13 depending on the values of β [9,10,11,12].Even more recently the new data of the PTA collaborations have been released [13,14,15] together with the first determinations of the Chinese Pulsar Timing array (CPTA) [16].The Parkes PTA collaboration considered 30 millisecond pulsars spanning 18 years of observations; they estimated q 0 = 3.1 1.3 −0.9 with 2 The scale factor is normalized at the present conformal time as a(τ0) = a0 = 1 so that comoving and physical frequencies coincide for τ → τ0.Moreover the standard prefixes of the International System of Units are consistently employed, e.g. 1 aHz = 10 −18 , Hz, 1 nHz = 10 −9 Hz and similarly for all the other frequency domains of the spectrum. 3A generic pulsar timing array is in fact a series of millisecond pulsars that are monitored with a rhythm that depends on the specific experiment.We refer here, in particular, to the NANOgrav collaboration [9], to the Parkes Pulsar Timing array (PPTA) [10], to the European Pulsar Timing array (EPTA) [11]; the PTA data have been also combined in the consortium named International Pulsar Timing array (IPTA) [12]. 4Even if the various collaborations use indifferently γ and β we shall refrain from this practice since γ could be confused with a different physical quantity (denoted simply by γ) that shall be introduced later on in the discussion.
For the record we recall that the relation between γ and β is simply given by since β = (3 − γ)/2.
Even if the interpretation of various collaborations preferentially focuses on a background produced by diffuse astrophysical sources at the present time 5 , in this paper we are going to examine the qualitatively different hypothesis that the excesses summarized by Eqs.(1.1)-(1.2) are caused by the parametric amplification of the quantum fluctuations of the tensor modes of the geometry.Since the relic gravitons are only coupled to the evolution of the space-time curvature Grishchuk [18,19], Ford and Parker [20] and some others first outlined that the formation of a diffuse backgrounds of gravitational radiation are a generic consequence of quantum mechanics in cosmological spacetimes.Depending on the early variation the space-time curvature, relic gravitons are then expected in different dynamical situations and, in particular, during a stage of de Sitter and quasi-de Sitter expansion [21].
Considering that, at early times, the conventional expanding stages of the concordance paradigm are complemented by an inflationary epoch, h 2  0 Ω gw (ν, τ 0 ) is a monotonically decreasing function of the comoving frequency from the aHz range to the MHz domain.More precisely between few aHz and 100 aHz the spectral energy density in critical units scales approximately as ν −2 [25,26,27]; for higher frequencies h 2 0 Ω gw (ν, τ 0 ) exhibits a quasi-flat plateau whose precise slope depends in fact upon the parameter r T measuring the ratio between the tensor power spectrum and its scalar counterpart.In the case of single-field inflationary scenarios (often taken as the benchmark class of models by some observers [22,23,24]) the consistency relations stipulate that the tensor spectral index m T and the slow-roll parameter ϵ k are both determined by the value of r T according to the following (approximate) chain of equalities m T ≃ −r T /8 ≃ −2ϵ k that are customarily referred to as the consistency conditions.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 that coincides, in what follows, with k p = 0.002 Mpc −1 .The addition of a single tensor component (only described by r T ) allows for an accurate set of limits implying that r T ≤ 0.06 or even r T ≤ 0.03 [22,23,24].Moreover, the pivot wavenumber k p corresponds to ν p = k p /(2π) = 3.09 aHz and this is why the limits on r T (k p ) may be translated into constraints on the spectral energy density of the relic gravitons in the aHz range.
For typical comoving frequencies of the order of the nHz the current bounds on r T suggest that in the concordance paradigm h 2 0 Ω gw (ν, τ 0 ) must be O(10 −17 ) (or smaller); in this case all the relevant wavelengths exit in the inflationary stage and reenter during the radiation phase.Furthermore, the accurate estimate of h 2 0 Ω gw (ν, τ 0 ) in the nHz region involves a number of late-time effects that may reduce even further the overall amplitude of the spectral energy density: the corrections due to neutrino free-streaming [28] (see also, for instance, [29,30,31]) suppress h 2 0 Ω gw (ν, τ 0 ) for ν < ν bbn where ν bbn = O(10 −2 ) nHz is the typical frequency associated with big-bang nucleosynthesis.Thus in the nHz range the relic gravitons produced within the conventional lore are nine or ten orders of magnitude smaller than the figures reported in Eq. (1.2).This means that if the nHz excess is caused by relic gravitons amplified by the evolution of the space-time curvature the relevant time-scale of the problem is primarily given by τ k defining the moment at which the wavelength associated with ν ref crossed the comoving Hubble radius after the end of inflation.A simple estimate suggests that τ k is just a fraction of the time-scale of big-bang nucleosynthesis6 In what follows ν P T A denote the frequency of the PTA and it can be either slightly larger or smaller than ν ref ; when ν P T A > ν ref the corresponding wavelength crossed the comoving horizon even earlier.
The second relevant scale of the problem is given by the ratio between ν P T A and the expansion rate at the end of inflation: where A R denotes the amplitude of the curvature inhomogeneities at the pivot scale k p .From Eqs.
(1.3)-(1.4) it follows that any modification of the post-inflationary evolution is unlikely to produce a hump for frequencies O(ν ref ); interesting enhancements may arise for ν > mHz ≫ ν ref implying that the post-inflationary evolution is indeed able to increase the spectral energy density between the mHz and the MHz [32,33,34] but not in the nHz range.As we shall see Eqs. (1.3)-(1.4)ultimately imply that an excess comparable with Eq. (1.2) cannot be related to a post-inflationary modification of the comoving horizon.Before plunging into the analysis it is probably useful to stress that already after the first data releases of the PTA a large number of different (and sometimes opposite) explanations for the nHz excess have been proposed.A common characteristic distinguishing the hypothesis pursued in this paper and the ones propounded in the current literature is that the nHz excess is typically due the presence of late-time sources of anisotropic stress.An example along this direction is represented by cosmic strings whose oscillating loops may emit gravitational waves at different epochs ultimately producing a stochastic background [35] with quasi-flat spectral energy density which is typically larger than the inflationary signal.For the largest values of the string tension in Planck units there is the possibility of an excess in the nHz region (see e.g.[36,37]).Various subsequent analyses based on these observations have been proposed more recently with various degrees of success.Another late-time source of gravitational radiation is represented by strongly first-order phase transitions that are known to produce spikes at low and intermediate frequencies because of the partial breaking of homogeneity due to the nucleation of bubbles of the new phase.The amount of gravitational radiation produced by the phase transition depends chiefly on the difference between the energy density in the broken and in the symmetric phase.This energy density may be comparable with the energy density of the ambient plasma (and in this case the phase transition experiences a strong supercooling) or smaller than the energy density of the surrounding radiation (and in this case the phase transition is mildly supercooled).When the gravitational radiation is produced from the collisions of the bubbles of the new phase [38,39,40] the spectral energy density scales like ν 3 , reaches a maximum and then decreases with a power that may be faster than ν −1 .The spectral energy density inherits also contribution from the sound waves of the plasma [41] and this second component may be even larger than the one due to bubble collisions.There are two points that make this explanation difficult.The first one is that the powers of the hump are typically steeper than the ones suggested by the PTA observations.The second observation concerns the physics of the phase transition.The bursts of gravitational radiation coming from the TeV are absent at frequencies much lower than the µHz also because the electroweak phase transition is not strongly first-order given the values of the Higgs mass.This implies the necessity of more contrived explanations (see e.g.[42] and references therein).There exist also more exotic possibilities where the nHz excess is caused by modifications of gravity in alternative cosmological scenarios; see, in this respect, [43] and references therein.
The layout of this investigation is the following.In section 2 the post-inflationary modifications of the expansion rate (related to the reentry of λ P T A ≃ 1/ν P T A ) are analyzed and confronted with the nHz excesses.In section 3 we instead examine the exit of the PTA scale during inflation and focus on the possibility of a refractive index coming from the interactions of the relic gravitons with the background geometry.In section 4 the general features of the bouncing scenarios (modifying both the exit and the reentry of λ P T A ) are specifically examined in the light of a potential nHz signal.Section 5 contains our concluding remarks.Appendix A covers the main notations employed in the paper and appendix B is instead focussed on the general form of the effective action of the relic gravitons and on its different parametrizations.

The comoving horizon after inflation 2.1 General considerations
When all the wavelengths of the relic gravitons are shorter than the comoving Hubble radius at the present time, the spectral energy density in critical units can be expressed as: where H = a H and H = a ′ /a; the notations employed here are summarized in appendix A and, in particular, in Eqs.(A.2) and (A.10)-(A.11).Equation (2.1) has been deduced in the limit a re ≫ a ex which is valid for all conventional and unconventional inflationary scenarios even if, in some cases, Ω gw (k, τ ) may also depend on the intermediate evolution between τ ex and τ re since Q k (τ ex , τ re ) (whose square modulus enters Eq. (2.1)) contains an integral over the conformal time coordinate τ : From the evolution of the mode functions of Eq. (A.15), τ ex and τ re are the roots of the equation that ultimately defines the different regimes of the evolution of the mode functions.In a stricter mathematical perspective, τ re and τ ex define the turning points where the solutions of Eq. (A. ) are both negligible in the limit kτ ≫ 1 when all the relevant wavelengths are inside the comoving horizon.Therefore, if all the wavelengths exit during inflation and reenter in a stage dominated by radiation (as it happens in the concordance paradigm and in Fig. 1) Eq. (2.2) implies that Q k (τ ex , τ re ) → 1.There can be however physical contingencies where the contribution of the integral in Eq. (2.2) gets larger than 1 and, as long as the exit is a regular turning point that the integrals appearing in Eq. (2.2) can be estimated as: where x = τ /τ ex and H ex = a ex H ex ≃ k.In Eq. (2.4) we also included two subscripts in the integrand with the purpose of stressing that the corresponding contributions may arise during the inflationary and in the post-inflationary stage; we also conventionally assumed, for the sake of concreteness, that the accelerated stage lasts up to −τ 1 and it is replaced for τ > −τ 1 by a decelerated evolution.

The post-inflationary evolution of the comoving horizon
The ratio between the comoving horizon and λ P T A is illustrated in Fig. 1 when the relevant bunch of wavelengths exits the comoving horizon during inflation (i.e.ϵ ex ≪ 1) and reenters in a radiationdominated stage of expansion (i.e.ϵ re → 2); in this case the spectral slope of Ω gw (ν, τ 0 ) follows from Eqs. (2.1)-(2.4)and by just focussing, for simplicity, on the k-dependence we can deduce the spectral slope (2.5) In Eq. (2.5) we used that a ex H ex = −1/[(1 − ϵ ex )τ ex ] and kτ re ≪ 1; moreover, by enforcing the validity of the consistency relations (i.e.r T = 16ϵ k ), the spectral index m T is given by7  The ratio between the comoving horizon and the wavelength λ P T A = O(0.30)pc is illustrated; as already mentioned λ P T A corresponds to the comoving frequency ν P T A = 31.68nHz.Common logarithms are employed on both axes and the two blobs indicate the exit and to the reentry of the corresponding wavelength.The full line is associated with ν P T A while the dot-dashed and the dashed lines illustrate, respectively, the situations where the frequency is slightly larger (i.e. 10 ν P T A ) or slightly smaller (i.e.0.1 ν P T A ) than the benchmark value provided by ν ref .
The timeline reported in this cartoon is characteristic of the concordance paradigm where all the wavelengths (including λ P T A ) do their first crossing during inflation (i.e. for a < a 1 ) and reenter in the radiation-dominated stage.
Equation (2.6) is a consequence of Fig. 1 where the first crossing of λ P T A occurs during inflation and the reentry takes place for ϵ re → 2; this is however not the most general situation even assuming that the left portion of Fig. 1 (i.e. a < a 1 ) is not modified.Indeed for a > a 1 , the post-inflationary expansion rate can be modified as (a H) −1 ∝ a 1/δ with δ ̸ = 1 and this case is illustrated in Fig. 2 when the usual radiation epoch (with δ → 1) is either preceded by a stage expanding faster than radiation (i.e.δ > 1) or by one expanding at a rate slower than radiation (i.e.δ < 1).While more complicated possibilities are swiftly examined at the end of the present section, already these two opposite situations lead to a modified Ω gw (k, τ ) above ν bbn [32,33,34] and may potentially lead to an excess in the nHz range.The dashed curve of Fig. 2 corresponds to the comoving horizon in the limit δ → 1 while for the two remaining profiles the radiation phase is preceded by a stage where δ ̸ = 1.In the three cases the minima depend both on the length of the post-inflationary stage and on its expansion rate: where H r and H 1 denote, respectively, the Hubble rates at the onset of the radiation stage and at the end of inflation (see e.g.Eq. (A.2)).Since, by definition, ξ < 1 the comoving horizon at its minimum is comparatively larger for δ > 1 than for δ → 1.For the same reason the opposite is true when δ < 1 and (a 1 H 1 ) −1 /λ P T A gets smaller than in the limit δ → 1.In Fig. 2 the λ P T A crosses the comoving Hubble radius when the expansion rate is different from radiation and this is why, according to Eq. (2.1) the spectral energy density gets modified.1 the ratio between the comoving horizon and λ P T A is illustrated with the difference that, in the present cartoon, the post-inflationary evolution is not immediately dominated by radiation; the curve δ → 1 (dashed line) is also reported for comparison.The postinflationary stage with δ ̸ = 1 lasts until the crossing time of λ P T A .On the one hand this choice maximizes the deviations from the standard form of the spectral energy density, on the other hand the modified post-inflationary evolution cannot last much longer; this is because the crossing time of λ P T A is just a fraction of τ bbn (see Eq. (1.3) and discussion therein).
If the λ P T A reenters the Hubble radius when δ ̸ = 1 (as it may happen in Fig. 2) ϵ re = O(1) while, as before, ϵ ex ≪ 1.For this reason, unlike the standard case illustrated in Fig. 1, |kτ re | −2 = O(1) in Eq. (2.1) and the spectral slope is not given by m T (as in Eq. (2.6)) but rather by the intermediate spectral index n T : The consistency relations remain valid since, for a < a 1 , Figs. 1 and 2 share the same inflationary evolution.We note however that for δ > 1/2 the integral appearing in Eq. (2.4) gives a subleading contribution and Q k (τ ex , τ re ) → 1. Conversely, when δ < 1/2, we have instead Q k (τ ex , τ re ) ≃ 1 + |k τ 1 | 1−2δ so that the second contribution dominates for all the frequencies of the spectrum (i.e. for kτ 1 < 1).Finally, when δ → 1/2 the integral of Eq. (2.4) inherits a logarithmic correction which is relevant in specific models [32,33,34].
Unless the relic gravitons would lead exactly to the same slope of the astrophysical foregrounds associated with black-hole binary systems, the value β = −2/3 is not particularly compelling in a cosmological setting.In the general case (i.e. when the special value β = −2/3 is not preliminarily selected) the Parkes PTA collaboration [13] suggests that β = −0.45± 0.20.This determination is only marginally compatible with the value of Eq. (2.11) in the limit δ ≥ 1/2.The discrepancy between the observational determination of β and the values predicted by Eq. (2.11) becomes even more significant if we look at that NANOgrav data suggesting [14] β = −0.10 ± 0.30.Thus if we stick to the general situation suggested by the observational collaborations the limit δ < 1/2 should also be adequately considered in Eq. (2.8) so that Eq. (2.11) is ultimately replaced by: From Eq. (2.11) we have δ = −β+O(r T ) while in the case δ < 1/2 Eq.(2.12) implies δ = 1+β+O(r T ).
From the profile of Fig. 2 with the help of Eqs.(2.1)-(2.2) and (2.8) the spectral energy density at high-frequencies becomes: where N (r T ) includes the effect of the suppressions associated with the low-frequency transfer function, with the neutrino free-streaming [28,29,30] and with the other late-time sources of damping (like the one associated with the dark-energy dominance [17]).For r T = 0.03 we can numerically estimate that N (r T ) = 10 −16.8 [51]; we then safely consider 8 for the present ends N = O(10 −17 ) with a theoretical amplitude at ν ref that ultimately depends upon ν r (see Eq. (2.13)) ν r = ν max ξ = 271.88g s, eq g s r 1/3 g ρ, r g ρ eq ξ MHz, ( where the ratio ξ = H r /H 1 has been already introduced in Eq. (2.7) while ν max given by: ν max = 271.88g s, eq g s r 1/3 g ρ, r g ρ eq MHz. (2.15)In spite of the specific value given in Eq. (2.14), ν r cannot be smaller than ν bbn and this remark turns out to be quite relevant for the comparison between the observed excesses and the theoretical expectations.Note that A R is the amplitude of curvature inhomogeneities at the pivot scale k p (see also Eq. (1.4) and discussion thereafter).The two filled rectangles represent the regions probed by the Parkes PTA and by NANOgrav in the plane (log q 0 , β).Since the two diagonal lines do not overlap with the shaded areas appearing in the lower portion of the plot, the amplitudes and the slopes of the theoretical signal cannot be simultaneously matched with the corresponding observational determinations.

Theoretical expectations and observed excesses
Common logarithms have been employed on the horizontal axis since the range of the observational determinations of q 0 spans nearly one order of magnitude.
Even if the theoretical and the observed slopes can be compatible for specific values of δ, the corresponding amplitudes involve rather different orders of magnitude and to analyze this essential aspect we may impose that Eqs.(2.10) and (2.13) coincide at ν ref In spite of the equality sign, the left-hand side of Eq. (2.16) turns out to be systematically smaller than the right-hand side and the two sides of the equation are in agreement only if ν r is much smaller than ν ref while, at the same time, n T = 2 + 2β is sufficiently large and positive.A large enough value of n T guarantees a sharp increase of the spectral energy density while the condition ν r ≪ ν ref makes wider the frequency range of the potential growth.Since the minimal value of ν r is provided by ν bbn we can select the most favourable situation and posit ν r = O(ν bbn ).For different values of N (r T ) (see Eq. (2.13) and discussion thereafter) Eq. (2.16) leads therefore to a specific relation between β and log q 0 : To close the circle, the result of Eq. (2.18) must then be compared in the plane (log q 0 , β) with the ranges of β and q 0 determined by the PTA collaborations [13,14].The two filled rectangles in Fig. 3 correspond, in this respect, to the observational ranges of q 0 and β; in the same plot the relation between β and log q 0 has been illustrated as it follows from Eq. (2.18) for two neighbouring values of N (r T ).The two diagonal lines in Fig. 3 imply that the values of β required to obtain h 2 0 Ω gw (ν ref , τ 0 ) of the order of 10 −8 or 10 −9 should be much larger than the ones determined observationally and represented by the two shaded regions.Since the full and dashed lines of Fig. 3 do not overlap with the two rectangles in the lower part of the cartoon, we can conclude that the excess observed by the PTA collaborations cannot be explained by the modified post-inflationary evolution suggested of Fig. 2. For the sake of accuracy we may separately analyze the case β = −2/3 and, in this situation, Eq. (2.16) becomes ) follows from the profile of Fig. 2 when the postinflationary expansion rate is slower than radiation since only in this case the slope of h 2 0 Ω gw (ν, τ 0 ) increases, as required by Eq. (2.10).There is however a complementary possibility stipulating that the comoving horizon prior to radiation dominance consists of two successive stages expanding at different rates and this situation is illustrated in Fig. 4 where the two consecutive post-inflationary stages (characterized by the rates δ 1 and δ 2 ) precede the ordinary radiation-dominated epoch.Even the standard sources considered before (satisfying δ i ≥ 1/2 for i = 1, 2) imply that when δ 1 < 1 and δ 2 > 1 (see the left plot of Fig. 4) the spectral energy density h 2 0 Ω gw (ν, τ 0 ) develops a trough at the intermediate frequency ν 2 ≫ ν r ≥ ν bbn .In other words for ν < ν 2 the spectral energy density decreases while it increases above ν 2 .The presence of a trough in h 2 0 Ω gw (ν, τ 0 ) corresponds to an expansion rate that is first slower (i.e.δ 1 < 1) and then faster (i.e.δ 2 > 1) than radiation.In this case, however, λ P T A crosses the comoving horizon when δ 2 > 1 and the spectral index for the corresponding bunch of wavelengths is negative whereas Eq. (2.10) suggests it should be positive.Instead of a single spectral index n T there will now be two different spectral indices n 1 and n 2 .For the profile appearing in the left plot of Fig. 4 the spectral index n 1 is positive since δ 1 < 1. However for ν < ν 2 the spectral index will be instead given by: This means that the spectral index is positive for ν > ν 2 and negative when ν < ν 2 with a trough in ν = O(ν 2 ).Consequently the spectral energy density is even more suppressed than in the case of the concordance paradigm and the corresponding timeline of the comoving horizon is not relevant for the present discussion where would instead aim at a large signal in the nHz range.The evolution of the comoving horizon is again illustrated in units of λ P T A but the postinflationary evolution is now characterized by two different expansion rates prior to the radiationdominated phase.In the plot at the left we have that the expansion rate is initially slower than radiation (i.e.δ 1 < 1) and then faster than radiation (i.e.δ 2 > 1).In the right plot the timeline is inverted and δ 1 > 1 while δ 2 < 1.In the left plot λ P T A crosses the comoving horizon when the plasma expands faster than radiation and this is why, for ν = O(ν ref ) the spectral energy density decreases in frequency.Conversely, in the right plot, the expansion rate is slower than radiation when λ P T A crosses the comoving horizon and, in this case, h 2 0 Ω gw (ν, τ 0 ) inherits a growing frequency spectrum.
For the present ends the relevant timeline is illustrated in the right plot in Fig. 4 where δ 1 > 1 and δ 2 < 1: the expansion rate is initially faster than radiation and then it gets slower so that h 2 0 Ω gw (ν, τ 0 ) decreases above ν 2 while it increases for ν < ν 2 .This means that the spectral energy density in critical units develops a hump for ν = O(ν 2 ); the signs of the spectral indices appearing in Eqs.(2.20) and (2.21) are exchanged: since δ 1 > 1 we have that n 1 < 0 while n 2 > 0 because δ 2 < 1 (see always the right plot of Fig. 4).In this case the wavelengths smaller than λ P T A cross the comoving horizon when the the expansion rate is controlled by δ 2 < 1 and, for frequencies ν = O(ν ref ), the spectral energy density increases as: where, as already stressed, n 2 > 0. The frequency ν 2 follows from profile of the comoving horizon and it is given by where, by definition, ξ 1 < 1 and ξ 2 < 1 denote the ratios of the curvature scales at the end and at the beginning of each expanding stage that precedes the ordinary radiation phase.The expressions of ν r and ν max coincide, respectively, with Eqs.(2.14) and (2.15) since, by definition, As before the analytic expression of ν 2 is not essential for the present purposes.Indeed, in spite of the expression of ν 2 the largest signal will be obtained when ν 2 = ν ref since for ν > ν 2 the spectral energy density decreases as: where we introduced the absolute value of the spectral index n 1 to stress that h 2 0 Ω gw (ν, τ 0 ) is suppressed for ν > ν 2 with a negative spectral index n 1 = 2(1 − δ 1 ) < 0 (i.e.δ 1 > 1).We can now try to match the amplitude and the slope of Eq. ( 2 Again the widest frequency range corresponds to the case when ν r → ν bbn and after some algebra we get back exactly to the same condition of Eq. (2.18).Due to the smallness of the potential signal, the direct bounds coming from wide-band interferometers [52,53] and the indirect constraints from big-bang nucleosynthesis [55,56,57] do not play a relevant rôle in the argument of this section 9 .The results obtained so far can then be summarized in the following manner: • if the post-inflationary evolution is modified prior to radiation dominance h 2 0 Ω gw (ν, τ 0 ) may increase in comparison with the concordance paradigm for typical frequencies ν = O(ν ref ); this happens only if the wavelengths O(0.3) pc (roughly corresponding to comoving frequencies O(30) nHz) cross the comoving horizon when the expansion rate is slower than radiation; • the PTA signals imply a growing spectral energy density (i.e. 2 + 2β > 0) and this is consistent with an expansion rate that is slower than radiation at least when the wavelengths O(λ P T A ) reenter the comoving horizon; • however the amplitudes and the slopes of the theoretical signal do not match simultaneously the observational determinations of the PTA in the (log q 0 , β) plane.
A hump in h 2 0 Ω gw (ν, τ 0 ) for ν = O(ν ref ) indeed follows when the post-inflationary expansion rate is slower than radiation but the theoretical signal is not consistent with the slope and with the amplitude of the observational spectrum.In specific cases the theoretical and the observed slopes are compatible but the corresponding amplitudes differ by 9 or even 10 orders of magnitude.As usual the agreement between theoretical and observed slopes is just a mere indication that is, per se, irrelevant for the final conclusion.
3 The comoving horizon during inflation

Modified evolution of the comoving horizon
If the refractive index n(a) of the relic gravitons is dynamical the condition of Eq. (2.3) defining the exit of a given wavelength is now replaced by: where the overdot denotes, as usual, a derivation with respect to the cosmic time coordinate; in Eq. (3.1) a(τ ) is actually rescaled as b = a/ √ n and ϵ is, in practice, the generalization of the slow-roll parameter ϵ already introduced in the previous section and in the appendices.The details on the connection between F and H are also discussed in Eqs.(B.5)-(B.6) of appendix B and the analog of (a H) −1 is now represented by (b F ) −1 whose evolution in units of λ P T A is illustrated in Fig. 5.By looking at the profiles of Fig. 5 there are in fact two complementary possibilities: • if λ P T A crosses the comoving horizon when the refractive index is dynamical h 2 0 Ω gw (ν, τ 0 ) may inherit a growing spectrum comparable with the PTA excess and in Fig. 5 the two curves (labeled, respectively, by 1 and 2) illustrate this possibility; while in the case 2 the first crossing occurs during the refractive phase, for the curve labeled by 1 the crossing occurs nearly at the end of it; • if the first crossing takes place after the end of the refractive phase (see the curve labeled by 3 in Fig. 5) the spectral energy density does not show any appreciable excess and the resulting spectral energy density is quasi-flat.
Even if the curve 3 in Fig. 5 does not lead to a growing spectrum for λ P T A it is anyway relevant for shorter wavelengths λ ≪ λ P T A ; for the corresponding frequencies h 2 0 Ω gw (ν, τ 0 ) may have a flatter slope but also further spikes caused by the post-inflationary evolution.Both structures are significantly constrained in the audio and in the MHz bands.

Refractive index and effective action
When the refractive index of the relic gravitons increases during a conventional stage of inflationary expansion the spectral energy density is blue at intermediate frequencies (typically above the fHz) and then flattens out after a knee that is generally smaller than the mHz [58].The general shapes of h 2 0 Ω gw (ν, τ 0 ) suggest that this possibility is particularly interesting in the light of the PTA excesses.The comoving horizon is illustrated in units of λ P T A when its early evolution is modified during the inflationary stage; as usual common logarithms are employed on both axes.The dynamics of the refractive index of the relic gravitons is responsible for the modified evolution illustrated in this cartoon.If the crossing of λ P T A occurs during the refractive phase (or at the end of it) the spectral energy density inherits a blue or violet spectrum that may eventually explain, under very specific conditions, the PTA excess.
When the refractive index is dynamical the action of the relic gravitons (see e.g.Eq. (A.3) and discussion therein) is modified as described in appendix B and it can be written as: 2) The analysis of Eq. (3.2) is greatly simplified if the conformal time coordinate is redefined from τ to η where the relation between the new and the old time parametrizations implicitly follows from n(η) dη = dτ .Equation (3.2) becomes then canonical in terms of a redefined scale factor conventionally denoted by b(η) [58]: The result of Eq. (3.3) generalizes the standard Ford-Parker action discussed in appendix A to the case of a dynamical refractive index and it explains how and why the evolution of the tensor modes is modified even during a conventional stage of inflationary expansion.A number of different physical reasons may lead to an effective index of refraction of gravitational waves in curved spacetimes [58,59,60].For instance the effective action of single-field inflationary models involves all the terms that include four derivatives 10 and are suppressed by the negative powers of a large mass scale [61].Another possible origin of a refractive index are non-generic models of inflation where the higher-order corrections assume a specific form since the inflaton has some particular symmetry (like a shift symmetry φ → φ + c) or because the rate of inflaton roll remains constant (and possibly larger than 1), as it happens in certain fast-roll scenarios [65,66,67].There are also the cases where the higher-order curvature corrections are given in terms of the Gauss-Bonnet combination weighted by some inflaton dependent-coupling [68,69,70].In [58] (see also [71,72,73]) it has been argued that in all these situations the effective action of the relic gravitons can be modified and ultimately assumes the general form discussed in appendix B.
For the present purposes what matters is not so much the origin of the refractive index but rather the possibility that its dynamical evolution could lead to a nHz excess.We then assign n(a) and even if the phase velocity of the relic gravitons is not required to be sub-luminal we impose, for consistency, that n(a) ≥ 1.Moreover, since the contributions to n(a) arise from diverse physical considerations we prefer to reverse the problem by focussing the attention on those profiles that eventually lead to potential excesses in the nHz range.Along this perspective we are led to consider an appreciable change of the refractive index during inflation with the concurrent requirement that n(a) reaches 1 in the standard decelerated stage of expansion 11 : where a i and a 1 denote, respectively, the beginning and the end of the inflationary epoch; a * indicates the boundary of the refractive stage.The three successive physical regimes described by Eq. (3.4) are, in some sense, more relevant than the specific analytic form that is however quite useful for numerical estimates.When a ≫ a 1 we have that n(a) → 1 and the sharpness of the transition is controlled by the parameter d ≥ 1.In the range a * < a < a 1 n(a) is constant but still larger than 1 (i.e.n(a) ≃ n * > 1) and, finally, when a < a * the refractive index is truly dynamical since n(a) ≃ n * (a/a * ) α .

The spectral energy density in critical units
We start by examining the spectral energy density for typical wavenumbers k < a * H * where Ω gw (k, τ ) is actually increasing.Using Eqs.(3.3)-(3.4) the spectral energy density in critical units can be obtained from Eqs. (2.1)-(2.2) with few relevant modifications: (3.5) Equation (3.5) is valid under the assumption that η re = τ re so that the reentry of the relevant wavelengths occurs when the refractive index is not dynamical; furthermore if τ re falls within the radiation phase (i.e. a ′′ → 0) we also have kτ re ≪ 1 in Eq. (3.5).Since any wavelength exiting for η < −η * does its first crossing during the inflationary phase, the corresponding refractive index is n = n * (a/a * ) α ; moreover the explicit expression of Q k (η ex , η re ) is also slightly more general than in the case of Eq. (2.2) 11 Some other possibilities have been considered in [58] and these cases can be easily added but they will not be examined here for the sake of conciseness.In this sense we shall regard the profile of Eq. (3.4) as the minimal example that potentially leads to a nHz excess.
Finally, from Eq. (3.3) b(η) can be expressed in the following manner: The result of Eq. (3.7) is valid if all the wavelengths O(λ P T A ) exit while the refractive index is still dynamical (as illustrated in the curves 1 and 2 of Fig. 5).Inserting now Eq.(3.7) into Eq.(3.5) a more explicit form of the spectral energy density can be deduced: g s, eq g s, r where g ρ and g s are, respectively, the effective number of relativistic species associated with the energy and with the entropy density.As usual Ω M 0 and Ω Λ denote the present critical fractions of matter and dark energy.It is well known that the dominance of dark energy suppresses the spectrum by a factor (Ω M 0 /Ω Λ ) 2 = O(0.1)(see, for instance, [17]).In Eq. (3.8) ν * denotes the frequency of the spectrum associated with η * and since k * = 1/η * the corresponding comoving frequency is: In Eq. (3.10) N * = ln (a * /a i ) is the number of e-folds during the refractive stage while N t = ln (a 1 /a i ) denotes the total number of e-folds; finally ν max indicates the maximal frequency of the spectrum and it coincides with Eq. (2.15) since, so far, the radiation dominance starts right after the end of inflation.The spectral index n T appearing in Eq. (3.8) depends on α and on ϵ k and it is: The tensor spectral index of Eq. (3.11) applies in the intermediate frequency range when the corresponding wavelengths exit during inflation and reenter in the radiation phase; in Eq. (3.11) α is always much larger than ϵ k ≃ r T /16 ≤ 0.03/16 ≪ 1 so that the exact result can be accurately evaluated in the limit ϵ k ≪ 1.While Eqs. (3.8) and (3.11) hold for ν < ν * , the spectral energy density can also be evaluated in the range ν * < ν < ν max (i.e. a * H * < k < a 1 H 1 ) corresponding to wavelengths that exited the comoving horizon when the refractive index was already constant (see the curve 3 in Fig. 5).Since the corresponding wavelengths do their first crossing when the refractive index is already constant (i.e.n → n * ) we have that η ex = τ * /n * and the spectral energy density becomes: where the spectral index is given by m T = −2ϵ k = −r T /8 and g s, eq g s, r Equation (3.12) evaluated for ν = ν * corresponds exactly to Eq. (3.8) computed at the same reference frequency and the equivalence of the two expressions ultimately follows from Eq. (3.10).Furtheremore, in Eqs.(3.8) and (3.11) (H 1 /M P ) 2 can be traded for π ϵ k A R where A R is the amplitude of curvature inhomogeneities at the pivot scale k p (see also Eqs. (1.4) and (2.15)).It is finally worth recalling that, for a standard thermal history, g s, eq = 3.94 while g ρ, r = g s, r = 106.75 in Eqs.(3.9) and (3.13).
Before proceeding further it is useful to comment on the possible range of variation of N t and N * appearing in Eq. (3.10).The value of N * measures the range of variation of the refractive index during inflation and, for this reason, N * < N t .The total number of e-folds is, in principle, arbitrary but a useful benchmark value is notoriously given by N t = O(60).This figure coincides, approximately, with the number of e-folds elapsed from the moment where the CMB wavelengths crossed the Hubble radius during inflation.If we estimate these wavelengths with k −1 p (where k p = 0.002 Mpc −1 ) then we have that N k (i.e. the number of e-folds elapsed since the crossing of k −1 p ) is, approximately, O(60) for r T = 0.06; more precisely it can be shown that where, as before, ln denotes the natural logarithm and ϵ k is the value of the slow-roll parameter at the crossing of the given set of wavelengths.Equation (3.14) assumes that after inflation the evolution is always dominated by radiation even if in Fig. 12 this assumption will be relaxed.The value of N k is also of the order of N max (i.e. the maximal number of e-folds presently accessible to largescale observations); see in this respect the appendices of Ref. [54].This value follows by requiring that the redshifted inflationary event horizon fits within the present Hubble patch; in practice this means that where H i denotes the expansion rated during the initial stages of inflation.Neglecting for simplicity the evolution of the relativistic species of the plasma we get that N max = 61.55 for the same fiducial set of parameters employed in Eq. (3.14).Typical values of N t of the order of N k and N max define, in practice, the minimal duration of the inflationary stage.Conversely values of N t smaller than N k (or N max ) characterize the durations of inflationary stages that are comparatively shorter than the benchmark value of Eq. (3.14).As we shall see in a moment, the relatively short inflationary stages (where N t ≤ O(61)) seem to be preferred for a potential explanations of the PTA excesses.

Accounting for the PTA excesses?
Equations (3.5), (3.8) and (3.12) are now compared with the parametrizations of the PTA signal given in Eqs.(1.2) and (2.10).Since, by definition, the intermediate spectral index is given as 2 + 2β = n T Eq. (3.11) implies a relation that determines α as a function of ϵ k (or r T ) and β: Moreover, given that q 0 depends on all the other parameters determining the amplitude of Ω gw (ν, τ 0 ) (see Eqs. (3.8) and (3.12)), we can demand that β and q 0 fall within the phenomenologically allowed ranges and check if the results of Eqs.(3.8) and (3.12) are compatible with the empirical determinations of the PTA.According to the Parkes PTA the values of β and q 0 fall, respectively, in the We also remind that N * controls the length of the phase where the refractive index is effectively dynamical.In both plots we traded α for β at a fixed value of ϵ k = 0.0018 (see Eq. (3.15) and discussion therein); the value of ϵ k is related to r T = 0.03 by the consistency relations.While in the left plot we illustrated q 0 (β, N * ) for different values of N t , in the right plot β has been fixed to −2/3; in this case the allowed range of q 0 (N t , N * ) is slightly different from the one of Eq. (3.16) and it is given by 1.82 < q 0 (N t , N * ) < 2.29.Finally, the shaded regions in both plots are consistent with the higher-frequency bounds coming from the audio band.
(3.16) Equation (3.16) constrains the spectral energy density and the corresponding region of the theoretical parameters is illustrated in the left plot of Fig. 6 where we report q 0 (β, N * ) for different values of N t ; the shape of each shaded region directly follows by requiring 2.2 < q 0 (β, N * ) < 4.4 for the various N t mentioned in the plot.On a technical side we note that Eq. (3.15) has been used with the purpose of trading directly α for β at a fixed value of ϵ k .In the right plot of Fig. 6 β is however fixed (i.e.β → −2/3) and, for such a choice, also the range of q 0 must be adapted following the observational determinations (i.e.q 0 = 2.04 +0.25 −0.22 [13]).For this reason in the right plot of Fig. 6 the constraints can be directly examined in the plane (N t , N * ).For the sake of accuracy, in Fig. 7 we illustrated the region pinned down in the plane (α, log r T ) by the different values of β.The same analysis illustrated in the case of the Parkes PTA can be repeated for the NANOgrav determinations with slightly different results; the analog of Eq. (3.16) is now given by [14] −0.40 ≤ β ≤ −0.20, 3.7 < q 0 < 10.6. (3.17) While the range of β given in Eq. (3.17) is narrower than in Eq. (3.16), in the case of q 0 we observe the opposite: the allowed values of q 0 of Eq. (3.17) are comparatively larger than the ones  A second class of constraints determining the shaded regions of Figs. 6 and 8 is related to the direct bounds from the operating wide-band detectors.In particular we remind that the LIGO, Virgo and Kagra collaborations (LVK) reported a constraint [52] implying: Ω gw (ν, τ 0 ) < 5.8 × 10 −9 , 20 Hz < ν L < 76.6 Hz, (3.18) in the case of a flat spectral energy density; in the present notations ν L indicates the LIGO-Virgo-Kagra frequency.The limit of Eq. (3.18) improves on a series of bounds previously deduced by the wide-band interferometers (see Ref. [17] for a review of the older results); in particular in Ref. [53] the analog of Eq. (3.18) implied Ω gw (ν, τ 0 ) < 6 × 10 −8 for a comparable frequency interval and always in the case of a flat spectral energy density.The bound of Eq. (3.18) can be used also in Eq. (3.8) since at high-frequency the spectral energy density is nearly scale-invariant.The results of Ref. [52] report however a threefold bound for a handful of spectral slopes; in particular, if the spectral energy density is parametrized as the limits of Ref. [52] read Ω(0) < 5.8 × 10 −9 (valid in the case σ = 0), Ω(2/3) < 3.4 × 10 −9 (when σ = 2/3) and Ω(3) < 3.9 × 10 −10 (when σ = 3).As the value of σ increases the bound becomes more restrictive for a fixed reference frequency and the three previous results are summarized by the following interpolating formula: (3.17) for different values of N t .The right plot is obtained by fixing β → −2/3 and the shaded region corresponds to the range 1.8 < q 0 (N t , N * ) < 3.1; this is because for β → −2/3 the NANOgrav collaboration suggests q 0 = 2.4 0.7 −0.6 .As in Fig. 6 the high-frequency bounds from the audio band have been also imposed.

Shapes of the spectra and further possibilities
Instead of using the approximations employed above it is instructive to compute numerically the spectral energy density from the exact form of the mode functions 12 .
As an example in the two plots of Fig. 9 we considered two different values of β (i.e.β = −0.65 and β = −0.55).Given a specific value of β within the observational range, the previous results (see, in particular, Figs. 6 and 8) lead directly to the allowed values of N * and N t .If N * and N t are of the same order the refractive index stops evolving when inflation approximately ends and, in this case, it is impossible to get a large signal in the nHz range without jeopardizing the big-bang nuclosynthesis constraint [58,73].Conversely, when N * < N t the refractive index stops evolving well before the onset of the post-inflationary stage, i.e. when the background is still inflating deep inside the quaside Sitter stage of expansion.In both plots of Fig. 9, to ease the comparison, we selected N t = 61 while different values of N * are illustrated.In both plots, for the same choice of the parameters, we also illustrated (with an arrow) the PTA excess and the Ligo-Virgo-Kagra bound [52,53].The PTA The discussion of section 2 does not exclude the possibility of two concurrent modifications of the comoving horizon operating before and after the end of inflation [71,72].This viewpoint is explored in Fig. 10 where we consider the possibility that the refractive index stops its evolution well before the end of inflation (i.e.N * ≪ N t ); however, unlike the case of Fig. 9, the post-inflationary evolution includes a long phase expanding at a rate slower than radiation.The spectral energy density in critical units will therefore have three different slopes for ν > ν eq .In both plots of Fig. 10, at intermediate frequencies h 2 0 Ω gw (ν, τ 0 ) has the same intermediate slopes appearing in Fig. 9 (see also Eqs. (3.11)-and (3.15)).However, after the quasi-flat plateau, the spectral energy density exhibits a further increasing branch before the maximal frequency.The corresponding wavelengths left the comoving Hubble radius during inflation and reentered in the post-inflationary stage before radiation dominance.In Fig. 10 the high-frequency spectral slope is O(1) since during the post-inflationary stage the evolution is described by a stiff fluid with δ ≃ 1/2 implying that (a H) −1 ∝ a 2 .The main difference between the plots of Figs. 9 and 10 ultimately comes from the high-frequency shape.While in the case of Fig. 9 the most relevant constraint comes from the operating interferometers [52,53], in the case of Fig. 10 the bounds coming from big-bang nucleosynthesis [55,56,57] must be taken into account since they imply: where Ω γ0 is the (present) critical fraction of CMB photons.As it is well known, the limit of Eq. (3.21) also sets an indirect constraint on all the extra-relativistic species (and, among others, on the relic gravitons).If applied to massless fermionic species, the limit is often expressed for practical Figure 10: As in Fig. 9 we illustrate the common logarithm of the spectral energy density as a function of the common logarithm of the comoving frequency.In the two plots the value of r T is the same but the values of N t are slightly dissimilar.In the plot at the left N * = 14 while the three spectra correspond slightly different values of β.In the plot at the right β = −0.63 and the three curves illustrate the variation of N * .Since the effect of neutrino free-streaming has been included, in both plots R ν denotes the neutrino fraction.
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.21) is thus between 10 −6 and 10 −5 .It is interesting to point out that the spectra of Fig. 10 are sensitive both to the interferometric bounds [52,53] and to the the limits of Eq. (3.21): the wide-band detectors constrain the height of the intermediate plateau while Eq.(3.21) sets a bound on the integrated h 2 0 Ω gw (ν, τ 0 ) and, ultimately, on the (global) maximum of the spectral energy density.
4 Bounces of the scale factor and curvature bounces

Basic considerations
In bouncing scenarios the spectral index can be positive between the fHz and the Hz [17] as it happens in the presence of a dynamical refractive index.However, while in the previous section both the amplitudes and the slopes could be predicted with reasonable accuracy thanks to the underlying inflationary dynamics, for the bouncing case it is comparatively easier to estimate h 2 0 Ω gw (ν, τ 0 ) in the intermediate frequency region rather than in the high-frequency domain which is often associated with a regime of strong curvatures.Even in the absence of a detailed theoretical derivation of the corresponding slopes, the high-frequency normalization can be disambiguated by employing the constraints of the audio band [52,53] and the big-bang nucleosynthesis limits [55,56,57].
The bouncing dynamics is sometimes associated with a contracting stage (i.e.ȧ < 0 and H < 0), as historically suggested, along slightly different perspectives, by Tolman and Lemaître [74,75] (see also [76,77]).For the present purposes it is useful to distinguish (at a purely kinematical level) the bounces of the scale factor from the ones involving the extrinsic (Hubble) curvature.While H and Ḣ change sign (at least) once in the bounces of the scale factor, for the curvature bounces H is always positive but Ḣ changes sign (at least) once.Broadly speaking the spectral slope at intermediate frequencies (which is the relevant one for the discussion of the PTA excesses) is related to the wavelengths that crossed the comoving Hubble radius before either H or Ḣ changed their sign the first time 13 .In both situations we can express the spectral energy density in terms of a common template where m T is the high-frequency spectral index and Ω * accounts for the corresponding normalization.The notations employed in Eq. (4.1) are purposely similar to the ones of the previous section since, in both cases, Ω * controls the high-frequency normalization.The physical meaning of the two quantities is however slightly different since Ω * now depends on the maximal curvature scale at the bounce (of the order of H 1 ) and on a number of other late-times parameters; in the simplest situation (where all the wavelengths reenter during a radiation dominated stage) the value of Ω * is given by g s, eq g s, r In inflationary scenarios the analog of (H 1 /M P ) is fixed by the amplitude of the (adiabatic and Guassian) curvature inhomogeneities whereas in the case of Eq. ( 4.2) it is more productive to determine Ω * directly from the available phenomenological constraints and to confront the obtained templates with the PTA observations.In this respect the most relevant limits are the ones coming from the audio band (see Eqs.

Growing spectral slopes at intermediate frequencies
While Eq. (4.2) applies in the high-frequency regime where the amplitude of h 2 0 Ω gw (ν, τ 0 ) is constrained by the current phenomenological bounds, in the intermediate frequency range (i.e. for ν < ν * ) the slope of the spectral energy density is instead denoted by n T = n T (γ, δ): In Eq. (4.3) γ and δ control, respectively, the expansion (or contraction) rates at early and late times; the general form of the spectral index can be expressed as: where the value of γ is related to the expansion (or contraction) rate in the cosmic time coordinate as a(t) = a 1 (−t/t 1 ) γ and this parametrization is valid before the first zero of H or Ḣ.As in the previous sections also in Eq. (4.4) δ denotes the expansion rate in the conformal time coordinate during the decelerated stage after the bounce.For 0 < γ < 1 Eq.(4.4) gives the spectral slope when, prior to the bouncing regime (taking place for |t| < t 1 ), the background experiences a stage of accelerated contraction (i.e.ȧ < 0 and ä < 0).Conversely when γ < 0 the background expands and accelerates (i.e.ȧ > 0 and ä > 0) with growing Hubble rate (i.e.Ḣ > 0).If δ > 1/2 the contribution to n T (γ, δ) only comes from the term k 4 |a re /a ex | 2 in Eq. (2.1) while which is much larger than 1 for all the amplified modes of the spectrum (i.e. for kτ 1 ≪ 1).Since the second term in Q k (τ ex , τ re ) dominates for δ < 1/2 and for kτ 1 < 1 we have that the spectral energy density Ω gw (k, τ ) ∝ k 4 |a re /a ex | 2 |Q k (τ ex , τ re )| 2 has overall a spectral slope (valid for a generic value of δ) proportional to |δ − 1/2|.After matter-radiation equality the same analysis implies: Below O(100) aHz h 2 0 Ω gw (ν, τ 0 ) is much smaller than in the case of the concordance paradigm.In

The slopes and amplitudes of the PTA excesses
If all the wavelengths of the order of λ P T A do the second crossing in a radiation-dominated stage (i.e. when δ → 1) the spectral index of Eq. (4.4) becomes: which is pinned down by the Parkes PTA [13].In the left plot 0 ≤ γ < 1 and this values describe a stage of accelerated contraction, as it happens for the bounces of the scale factor.In the right plot we discuss the case γ < 0 (i.e.accelerated expansion with growing expansion rate) that occurs for the bounces of the scale factor.
When γ → 1/3 in Eq. (4.6) the spectral index n T (1/3, 1) → 3; in this case the stage of accelerated contraction may be driven, in the simplest situation, by the kinetic energy of a scalar field as it happens, for instance, in the dilaton-driven phase of pre-big bang scenarios [78,79,80,81].Similarly when γ → 0 we would have, from Eq. (4.6), that n T (γ, 1) → 2; this is the situation preferred in the context of ekpyrotic models [82,83].Whenever the expansion rate at the second crossing is eventually different from radiation (i.e.δ ̸ = 1) also the form of Eq. (4.6) is different.For this reason it is relevant to illustrate the general situation and this has been done in Fig. 11 for the two complementary ranges 0 < γ < 1 and γ < 0. Grossly speaking the results of Fig. 11 imply that the spectral slopes suggested by the PTA collaborations can be reproduced either when δ > 1 (if 0 ≤ γ < 1) or when 0 < δ < 1 (provided γ < 0).To scrutinize this point in further detail it is easier to express β directly as a function of γ and δ: Using Eq. (4.7) we can therefore limit the range of variation of β(γ, δ) and deduce the allowed domain of the parameters.Along this perspective, the shaded areas of Fig. 12 are determined by requiring −0.65 < β(γ, δ) < −0.25 as suggested by the Parkes PTA [13].In Fig. 12 the line δ = 1 intersects the allowed region only if γ > O(0.6).A similar conclusion follows from Fig. 13 where we examine the NANOgrav observations [14]; in this case, as already stressed in the previous sections, the allowed range of β is narrower (even if the corresponding range of q 0 is larger).The shaded area in both plots of Fig. 13 correspond to the interval −0.40 < β < −0.20 and the correct values of β seem to be reproduced for 0 < γ < 1 but only when δ > 1.The opposite is true in the case γ < 0 where the region δ < 1 seems comparatively wider.Figure 13: The contours are determined exactly as in the case of Fig. 12 but in the present situation we consider the NANOgrav intervals [14] for β rather than the ones of the Parkes PTA.Since the interval of β is narrower in the case of Ref. [14], the shaded regions of both plots are narrower than the ones illustrated in Fig. 12.As usual, in both figures the values of β(γ, δ) remains the same on each contour.
So far we discussed the spectral slopes and we must now consider the corresponding amplitudes.For this purpose we then examine Eq. (4.1) and express it in terms of the parametrization of Eq. (1.1): 6.287 × 10 −10 q 2 0 (ν/ν ref Since the slope n T (γ, δ) coincides with 2 + 2β(γ, δ) the amplitude q 0 (which is determined experimentally) not only depends on β but also on (ν * /ν ref ) and Ω * .The high-frequency amplitude of the spectral energy density appearing in Eqs.(4.1)-(4.3)can then be fixed to the largest value compatible with the limits of wide-band interferometers in the audio band [52,53].From the theoretical viewpoint q 0 becomes a function of β and of (ν * /ν ref ): Since we know from Figs. 11 and 12 that there exist regions in the (γ, δ) plane where β falls within the observed range, from Eq. (4.9) we can determine the range of variation of ν * .A careful analysis reported in Fig. 14 shows that ν * must be slightly smaller than ν ref .More specifically in Fig. 14 we illustrate the different contours of q 0 (β, ν * ) and the shaded areas correspond, respectively, to the regions pinned down, respectively, by the Parkes PTA (i.e.3.7 ≤ q 0 ≤ 10.6 in the left plot of Fig. 14) and by the NANOgrav PTA (i.e.2.2 ≤ q 0 ≤ 4.4 in the right plot of Fig. 14).

Excluding a single blue slope
The allowed region of the parameter space pins down a range ν * = O(ν ref ) but there are theoretical models where the transition to the decelerated regime is very short (i.e.ν * ≫ ν ref ) so that ν * can eventually be of the order of ν max .To investigate this situation we can therefore write the spectral energy density in the following form: and demand that n T (γ, δ) falls in the interval of slopes associated with the PTA; if this is the case, then n T (γ, δ) = 2(1 + β) > 0. If we combine Eq. (4.10) with the requirements of Eq. (3.21) we also have that Ω max is constrained as: Thanks to Eq. ( 4.11) we can first trade Ω max for Ω bbn and then we can also express q 0 in terms of β and y = (ν ref /ν max ): Equation (4.12) is illustrated in Fig. 15 where q 0 is viewed as a function of β and ν max ; on the vertical Figure 14: We illustrate the contours of constant q 0 (ν * , β) for a fixed value of Ω * .The labels appearing on the various contours correspond to the common logarithms of q 0 (ν * , β).The intervals of β match the ones of the corresponding observations.The shaded regions correspond, as indicated in each of the plots, to the ranges of q 0 determined by the Parkes [13] and by the NANOgrav [14] PTA.
axis of both plots we report the range of β compatible with each of the corresponding experiment (i.e. the Parkes [13] and the NANOgrav [14] PTA).The shaded regions define the allowed range of q 0 for each of the two experiments.The results of Fig. 15 are incompatible with Eq. (4.10) since, grossly speaking, we can expect ν max between few MHz and the THz [17].In the parametrization of Eq. (4.10) the common logarithm of (ν ref /ν max ) (reported on the horizontal axes of the plots of Fig. 15) should be between −20 and −15.On the contrary the allowed region in Fig. 15 is located for y = O(10 −2 ).)) we illustrate the contours of constant q 0 (β, ν ref /ν max ).The intervals of β match the ones of the corresponding observations and the shaded regions correspond, as indicated in each of the plots, to the ranges of q 0 determined by the Parkes [13] and by the NANOgrav [14] PTA.
If the PTA excesses are the result of a bouncing stage it is necessary that the wavelengths of the order of λ P T A do their second crossing during a decelerated stage not dominated by radiation.Moreover a bouncing model leading to a spectral energy density with a single blue slope (possibly ranging between the equality frequency and ν max ) is unable to account for the PTA excesses if all the phenomenological bounds are concurrently satisfied.Various bouncing scenarios have been constructed with the aim of either complementing or even challenging the conventional inflationary ideas.Bouncing scenarios appear by physical premises that are quite different and various reviews are currently available (see, for instance, [84,85,86,87,88]).It is not uncommon that the solutions obtained in a given framework are recycled in an entirely different situation.Recently it has been argued that conventional inflationary models are generically in tension with the swampland criteria [89] and this is often seen as a further motivation for bouncing dynamics.The indications of the PTA can be very precious in this context and they might even contribute to the long standing problem of bouncing scenarios, i.e. the origin of a Gaussian and adiabatic mode of large-scale curvature inhomogeneities [17].

Concluding remarks
During the last two decades a series of significant limits from the millisecond pulsars at intermediate frequencies (roughly corresponding to the inverse of the time-scale along which the various pulsars have been monitored) severely constrained the isotropic and random backgrounds of gravitational radiation.More recently different ensembles of millisecond pulsars have been scrutinized with regular cadence by the pulsar timing arrays (PTA in the bulk of the paper).Despite the different conclusions on the expected correlations in the pulse arrival times between pairs of pulsars (and despite the slightly dissimilar determinations of the spectral parameters), the competing experiments seem to suggest concurrent evidences of gravitational radiation with stochastically distributed Fourier amplitudes at a benchmark frequency O(30) nHz and with h 2 0 Ω gw (ν, τ 0 ) ranging between 10 −8 and 10 −9 .While the origins of the PTA excesses are still perplexing, in this paper we speculate that the relic gravitons are responsible of the observed signal.Even if a collection of late-time sources may eventually lead to a diffuse background of gravitational radiation, the relic gravitons are instead produced, by definition, at much earlier times and solely because the rapid variation of the space-time curvature.
The theoretical perspective explored in this investigation strongly suggests that the problem is not yet to fit (more or less reliably) the existing data in terms of a series of preferred scenarios but to understand preliminarily whether or not the observed excesses in the nHz range are compatible with a modified evolution of the comoving horizon since this is the only way the spectrum of relic gravitons at intermediate frequencies can be affected.The goal of this study is therefore not to endorse a specific model (or to pin down the likely values of a hypothetical spectral index) but to see, more modestly, if and how the relic gravitons could be associated with the nHz excesses.All in all the systematic approach developed in this paper propounds three complementary physical possibilities that have been carefully perused and that should be further analyzed in the near future.
• The most conventional option stipulates that the timeline of the comoving horizon is not modified during inflation so that the nHz excess is caused by the drastic change of the postinflationary expansion rate prior to big-bang nucleosynthesis.
• A second alternative implies a modified evolution of the tensor modes during a conventional inflationary stage as it happens, for instance, when the gravitons inherit an effective refractive index from the interactions with the geometry.
• We finally consider the possibility of an epoch of increasing curvature prior to the conventional decelerated stage of expansion and argue that this option is only reconcilable with the observed excesses provided the wavelengths crossing the comoving horizon at early times do not reenter in an epoch dominated by radiation.
In connection with these three complementary options the results obtained in this analysis are, in short, the following.
• A late-time modification of the comoving horizon may indeed alter the spectral energy density of the relic gravitons in the nHz range but the observed amplitudes and slopes are unfortunately neither compatible nor comparable with this minimal explanation.
• Conversely, if a dynamical refractive index evolves during a conventional inflationary phase of nearly minimal duration the corresponding h 2 0 Ω gw (ν, τ 0 ) compares well with the nHz hump even if the observational data coming from competing experiments pin down slightly different regions of the parameter space.
• In the context of the bouncing scenarios, it is finally possible to exclude a nHz excess associated with a single slope of the spectral energy density between 100 aHz and the GHz range.However, if the spectral energy density has a break in the nHz region the direct limits of the wide-band detectors constrain the high-frequency amplitude and, in this case, we could account for the hump provided the bunch of PTA wavelengths reenter the comoving horizon during a decelerated stage not yet dominated by radiation.
The three possibilities discussed here are therefore not mutually exclusive: in the bouncing case, the simplest way to match the amplitudes and the slopes of the spectra is to modify the comoving horizon also at late times.The potential nHz excesses make even more relevant the high-frequency determinations of the spectral energy density.In particular the direct bounds in the audio band are essential even if the largest signal of the relic gravitons is expected in the MHz and GHz domains where smaller detectors may play a crucial rôle, as repeatedly suggested in the past.While the observational aspects of the problem cannot be addressed with the theoretical approach reported here, it is nonetheless true that the physical interpretation of the results probably demands a conceptual framework (such as the one pursued in this analysis) that could clarify (and even exclude) the primeval origin of the nHz excesses.
As before, in Eq. (B.1) the parity-breaking terms associated with quadratic combinations involving either the dual Riemann or the dual Weyl tensors have been neglected; both terms would appear in the effective action and can polarize the backgrounds of relic gravitons [64] (see also [91]).Among the three function c i (τ ) (with i = 1, 2, 3) we have that c 1 (τ ) and c 2 (τ ) are related to the expanding dimensions while c 3 (τ ) may appear in the case of compact extra-dimensions [92].If we factor c 1 (τ ) in Eq. (B.1) the resulting expression will be given by: where n(τ ) and n(τ ) denote, respectively, the refractive indices associated with the expanding and with the compact dimensions [92], i.e. n(τ ) = c 1 (τ )/c 2 (τ ) and n(τ ) = c 1 (τ )/c 3 (τ ).Equation (B.2) simplifies after a rescaling of the background dependence and its final form becomes 14 : From Eq. (B.5) it follows that the crossing of a given wavelength occurs when k 2 = (∂ 2 η b)/b.This expression generalizes therefore the notion of the comoving horizon during the refractive phase.More specifically we may recall the connection between the derivations in the various time parametrizations introduced so far, namely where η, τ and t denote, once more, the η-time, the conformal time and the cosmic time coordinates.By using Eq.(B.6) in the condition k 2 = (∂ 2 η b)/b we obtain, after simple algebra, a condition similar to Eq. (3.1): (B.7) 14 Once the new parametrization has been introduced we can also rescale the background dependence so that b(η) = c1(η)/n(η) and rc(η) = n(η)/n(η).In the absence of a contribution from the internal dimensions (i.e.mc → 0) Eq. (B.3) reproduces exactly Eq. (A.3) for when n → 1 and c1(τ ) = a 2 (τ ).
In the limit n → 1 we get b → a, F → H and ϵ → ϵ.This is why it is appropriate to consider (b F ) −1 as the generalization of the comoving horizon during inflation.Different choices in Eq. (B.1) are in fact artificial since they are ultimately equivalent to the one of Eq. (B.2).For instance instead of factoring c 1 (τ ) we may factor c 2 (τ ).If we rescale c 2 (τ ) we simply get the analog of Eq. (B. µ−1 (−kη), (B.9) where µ = ζ + 1/2 and H µ (−kη) is the Hankel function of first kind [93] and |M| = π/2.From Eq. (B.9) the tensor power spectrum becomes: For α → 0 we have that r T → 16ϵ and the standard consistency relation is recovered; note that n T is given exactly by the expression already mentioned in Eq. (3.11); therefore in the limit α → 0 we have that n T → −2ϵ.The notations employed here imply that the definition of the power spectrum P T (k, η) (and equally for P T (k, τ )) directly follow from the quantum mechanical normalization of the corresponding field operators.In particular we have that As anticipated, the tensor power spectrum of Eq. (B.14) is defined as in Eq. (B.10).For a full quantum mechanical discussion of this class of problems problem see, for instance, Ref. [90].

ν
Figure 1:The ratio between the comoving horizon and the wavelength λ P T A = O(0.30)pc is illustrated; as already mentioned λ P T A corresponds to the comoving frequency ν P T A = 31.68nHz.Common logarithms are employed on both axes and the two blobs indicate the exit and to the reentry of the corresponding wavelength.The full line is associated with ν P T A while the dot-dashed and the dashed lines illustrate, respectively, the situations where the frequency is slightly larger (i.e. 10 ν P T A ) or slightly smaller (i.e.0.1 ν P T A ) than the benchmark value provided by ν ref .The timeline reported in this cartoon is characteristic of the concordance paradigm where all the wavelengths (including λ P T A ) do their first crossing during inflation (i.e. for a < a 1 ) and reenter in the radiation-dominated stage.

νFigure 2 :
Figure2: As in the case of Fig.1the ratio between the comoving horizon and λ P T A is illustrated with the difference that, in the present cartoon, the post-inflationary evolution is not immediately dominated by radiation; the curve δ → 1 (dashed line) is also reported for comparison.The postinflationary stage with δ ̸ = 1 lasts until the crossing time of λ P T A .On the one hand this choice maximizes the deviations from the standard form of the spectral energy density, on the other hand the modified post-inflationary evolution cannot last much longer; this is because the crossing time of λ P T A is just a fraction of τ bbn (see Eq. (1.3) and discussion therein).

Figure 3 :
Figure 3: The two straight lines illustrate Eq.(2.18)  for N (r T ) = 10 −17 (full line) and for N (r T ) = 10 −16 (dashed line).The two filled rectangles represent the regions probed by the Parkes PTA and by NANOgrav in the plane (log q 0 , β).Since the two diagonal lines do not overlap with the shaded areas appearing in the lower portion of the plot, the amplitudes and the slopes of the theoretical signal cannot be simultaneously matched with the corresponding observational determinations.Common logarithms have been employed on the horizontal axis since the range of the observational determinations of q 0 spans nearly one order of magnitude.

Figure 4 :
Figure4: The evolution of the comoving horizon is again illustrated in units of λ P T A but the postinflationary evolution is now characterized by two different expansion rates prior to the radiationdominated phase.In the plot at the left we have that the expansion rate is initially slower than radiation (i.e.δ 1 < 1) and then faster than radiation (i.e.δ 2 > 1).In the right plot the timeline is inverted and δ 1 > 1 while δ 2 < 1.In the left plot λ P T A crosses the comoving horizon when the plasma expands faster than radiation and this is why, for ν = O(ν ref ) the spectral energy density decreases in frequency.Conversely, in the right plot, the expansion rate is slower than radiation when λ P T A crosses the comoving horizon and, in this case, h 2 0 Ω gw (ν, τ 0 ) inherits a growing frequency spectrum.

ν
Figure5: The comoving horizon is illustrated in units of λ P T A when its early evolution is modified during the inflationary stage; as usual common logarithms are employed on both axes.The dynamics of the refractive index of the relic gravitons is responsible for the modified evolution illustrated in this cartoon.If the crossing of λ P T A occurs during the refractive phase (or at the end of it) the spectral energy density inherits a blue or violet spectrum that may eventually explain, under very specific conditions, the PTA excess.

* 3 Figure 6 :
Figure 6: The regions pinned down by the Parkes PTA are illustrated in the left plot in terms of the parameters appearing in Eqs.(3.8) and(3.12).The shaded areas follow by imposing the range of Eq.(3.16)  for different values of N t (denoting, as repeatedly mentioned, the total number of e-folds).We also remind that N * controls the length of the phase where the refractive index is effectively dynamical.In both plots we traded α for β at a fixed value of ϵ k = 0.0018 (see Eq. (3.15) and discussion therein); the value of ϵ k is related to r T = 0.03 by the consistency relations.While in the left plot we illustrated q 0 (β, N * ) for different values of N t , in the right plot β has been fixed to −2/3; in this case the allowed range of q 0 (N t , N * ) is slightly different from the one of Eq. (3.16) and it is given by 1.82 < q 0 (N t , N * ) < 2.29.Finally, the shaded regions in both plots are consistent with the higher-frequency bounds coming from the audio band.

Figure 7 :
Figure 7: The shaded region follows by requiring −0.65 < β < −0.25 (as implied by Eq. (3.16)) and by imposing the relation of Eq. (3.15).The thick line outside the shaded region corresponds to β = −2/3.If the consistency relations are imposed β is determined in the (α, r T ) plane.However, as the plot shows, the smallness of r T implies that ϵ k can be neglected for the determination of the spectral index in the region ν < ν * .

55 Nt = 53 Nt = 51 - 3 Figure 8 :
Figure 8: In the left plot we illustrate the region pinned down by the NANOgrav PTA in the plane (β, N * ) for different values of N t .Both plots can be compared with the ones of Fig. 6 where we discussed the case of the Parkes PTA.The shaded regions in the left plot follow by imposing Eq.(3.17) for different values of N t .The right plot is obtained by fixing β → −2/3 and the shaded region corresponds to the range 1.8 < q 0 (N t , N * ) < 3.1; this is because for β → −2/3 the NANOgrav collaboration suggests q 0 = 2.4 0.7 −0.6 .As in Fig.6the high-frequency bounds from the audio band have been also imposed.

N * = 14 N * = 16 N 61 PTA 61 Figure 9 :
Figure 9: We illustrate the common logarithm of the spectral energy density in critical units as a function of the common logarithm of the comoving frequency.In both plots N t = 61 but the values of β and r T do not coincide and they are indicated above each of the two cartoons.The arrows indicate the PTA signal for the spectral indices corresponding to the ones selected in each of the plots.The high-frequency region labeled by LVK refers to the Ligo-Virgo-Kagra bound that applies in the audio band.The increasing branch and the flat plateau corresponds, respectively, to the analytic estimate of Eqs.(3.5) and (3.8).

Figure 11 :
Figure 11: We illustrate the intermediate spectral index n T (γ, δ) appearing in Eqs.(4.3)-(4.4).In both plots γ indicates the expansion (or contraction) rates before the bounce; δ denotes instead the expansion rate at late times (after the bounce) when the intermediate wavelengths of the spectrum reenter the Hubble radius.In the left plot 0 ≤ γ < 1 and the evolution before the bouncing regime is characterized by an accelerated contraction.In the right plot γ < 0 implying a phase of accelerated expansion (with growing Hubble rate) before the bouncing regime.

Figure 12 :
Figure12: The values of β(γ, δ) are illustrated together with the observational constraints derived from the Parkes PTA data.The shaded region denotes, in both plots, the range −0.65 < β(γ, δ) < −0.25 which is pinned down by the Parkes PTA[13].In the left plot 0 ≤ γ < 1 and this values describe a stage of accelerated contraction, as it happens for the bounces of the scale factor.In the right plot we discuss the case γ < 0 (i.e.accelerated expansion with growing expansion rate) that occurs for the bounces of the scale factor.

Figure 15 :
Figure15: After imposing the big-bang nucleosynthesis constraints (see Eq. (4.11)) we illustrate the contours of constant q 0 (β, ν ref /ν max ).The intervals of β match the ones of the corresponding observations and the shaded regions correspond, as indicated in each of the plots, to the ranges of q 0 determined by the Parkes[13] and by the NANOgrav[14] PTA.

r 2 c
(η)m 2 c h i j h i j .(B.3) Equation (B.3) follows from Eq. (B.2) by first changing the time parametrization from τ (the conformal time coordinate) to η; the relation between the two time parametrizations is simply given byn(η)dη = dτ .Let us therefore consider the simplest situation where the refractive index increases during inflation as suggested in Eq. (3.4); in this case for a < a * we would have n(a) = n * (a/a * ) α (with α > 0) so that the relation between the conformal time coordinate τ and the η-time can be swiftly worked since dη = dτ /n(a).From the definition of η we therefore have: in Eq. (A.2), H = ȧ/a and the overdot denotes a derivation with respect to the cosmic time coordinate.The second equality in Eq. (B.4) follows after integration by parts since ε ≪ 1 and α = 0. Equation (B.4) also implies that a H n = −1/[(1 − ϵ + α)η]; note once more that when n → 1 we also have α → 0 and the standard relation a H = −1/[(1 − ϵ)τ ] is immediately recovered.In the η-time parametrization the evolution of the mode functions simplifies and it is given by