Gravitational wave astronomy and the expansion history of the universe

Contents

1. Introduction

1.1. Ten years of gravitational wave astronomy

Gravitational waves have been predicted by Einstein in 1916¹ as a direct consequence of general relativity.² Later on this problem has been revisited by Einstein and Rosen with somehow contradicting conclusions³ suggesting that gravitational waves could be unphysical. While the legacy of Ref. 3 brought eventually some late skepticism on the true physical nature of gravitational radiation (see, for instance, Ref. 4) the gauge-invariant nature of gravitational waves has been well established in the 1970s.⁵ In spite of the pioneering attempts of Weber^6,7 and of the subsequent resonant detectors of gravitational radiation in the early 1970s, the first direct evidence of gravitational radiation dates back to the early 1980s when the orbital decay of a binary neutron star system has been originally observed.⁸ Roughly speaking almost one century after the first speculations, the gravitational waves have been detected by the wide-band interferometers.^9,10,11 The signals observed so far mainly come from astrophysical processes occurring at late time in the life of the Universe and they are the result of accelerated mass distributions with nonvanishing quadrupole moment. One of the most exciting directions is however related to the possible existence of diffuse backgrounds of gravitational radiation produced thanks to the early variation of the space–time curvature. This collection of random waves encodes a snapshot of the early expansion history of the Universe prior to the formation of light nuclei. The purpose of this topical review is to summarize what can be said on the early expansion history of the Universe from the analyses of the stochastic backgrounds of relic gravitational waves.

1.2. Gravitational waves in curved backgrounds

In the 1960s and 1970s it was believed that the tensor modes of the geometry could not be excited in curved background geometries. Although the chain of arguments leading to such a conjecture would be per se interesting, this misleading perspective implied that both electromagnetic and gravitational waves could be considered invariant for a Weyl rescaling of the four-dimensional background geometry; from a practical viewpoint Weyl invariance implies that both electromagnetic and gravitational waves should obey the same equations in a Minkowski background and in curved geometries eventually obtained by Weyl rescaling from a flat space–time.^12,13 This viewpoint persisted until the mid 1970s when it was challenged by a series of papers^14,15 suggesting that gravitational waves can be indeed excited in curved backgrounds and, more specifically, in Friedmann–Robertson–Walker (FRW) cosmologies.^16,17

Almost 50 years after these pioneering analyses the relic signals represent today a well-defined (and probably unique) candidate source for typical frequencies exceeding the kHz region where wide-band detectors are currently operating. Following the formulation of the inflationary scenarios^18,19,20,21 it became gradually clear that the conventional lore would predict a minute spectral energy density in the MHz region.^22,23,24 This is ultimately the reason why the most stringent tests of the conventional lore could come, in the near future, from the largest scales²⁵ where the limits on the tensor to scalar ratio $r_T$ are in fact direct probes of the spectral energy density in the aHz region. Throughout the discussions of this paper the standard prefixes of the international system of units are systematically employed; so for instance $1\mathrm{\text{kHz}}=10^3\mathrm{\text{Hz}}$, $1\mathrm{\text{aHz}}=10^{-18}\mathrm{\text{Hz}}$ and similarly for all the other relevant frequency domains mentioned hereunder.

1.3. The expansion history

During the last 30 years cosmology astrophysics and particle physics experienced a progressive unification toward two complementary paradigms accounting for the observations at small and large distance scales. The Standard Model of particle interactions describes the strong and electroweak physics or, as we could say for short, the microworld; although there are various hints on its possible incompleteness (typically related to the existence of dark matter), so far the Standard Model has not been falsified. The so-called concordance paradigm (based on general relativity) is customarily employed to analyze the macroworld of cosmological and astrophysical observations involving, in particular, the data associated with the temperature and polarization anisotropies of the Cosmic Microwave Background (CMB), the large-scale structure data and the supernova observations. The concordance paradigm is sometimes dubbed $\mathrm{\Lambda }$CDM where $\mathrm{\Lambda }$ accounts for the dark energy component and CDM stands for the cold dark matter. It is fair to say that, at the moment, the Standard Model of particle interactions and the $\mathrm{\Lambda }$CDM scenario seem mutually consistent but conceptually incomplete.

In the concordance paradigm the source of large-scale inhomogeneities is represented by the adiabatic and Gaussian fluctuations produced during a stage of conventional inflationary expansion. The subsequent evolutionary history of the plasma assumes a long period of expansion dominated by radiation until the epoch of matter–radiation equality and this timeline is broadly compatible with the idea that all the particle species were in thermal equilibrium above typical temperatures of the order of 200GeV but there is no direct evidence either in favor of this hypothesis or against it. In the past the radiation dominance of the primeval plasma before big bang nucleosynthesis (BBN) has been taken as a general truism also because it was practically impossible to check directly the early timeline of the expansion rate by simply looking at electromagnetic effects. This was the viewpoint conveyed in the pioneering analyses of the hot big bang hypothesis formulated by Gamow et al.^26,27,28 and subsequently confirmed with the discovery of the CMB²⁹ by Penzias and Wilson also thanks to the neat theoretical interpretation formulated by Peebles.^30,31 As we know the plasma became transparent to radiation around the time of photon decoupling. After that moment the slightly perturbed geodesics of the photons could be used to reconstruct the temperature and polarization anisotropies of the CMB³² but the electromagnetic signals coming from the earlier expansion history were quickly reabsorbed by the plasma and are today completely inaccessible to any direct detection.

The sensitivities of operating detectors^33,34,35,36 are notoriously insufficient to measure the diffuse backgrounds of relic gravitons but in the future new detectors might cover different frequencies³⁷ even beyond the so-called audio band ranging between few Hz and 10kHz. The gravitational waves produced thanks to the variation of the space–time curvature should then become an object of future empirical investigations even at high frequencies while at intermediate frequencies (in the nHz range) the backgrounds of relic gravitons could be observed by the pulsar timing arrays (PTAs)^38,39,40,41 that are now primarily focused on the diffuse astrophysical signals. We actually know that every variation of the expansion rate produces shots of gravitons with given averaged multiplicities and specific statistical properties. If these spectra will ever be detected the timeline of the expansion rate might be directly tested without the need of postulating a particular post-inflationary paradigm before the curvature scale of BBN whose striking success is the last certain signature of radiation dominance for typical temperatures smaller than $\mathcal{𝒪}(10)$MeV. When considering these possibilities at face value there are at least two conceptually different issues that must be addressed.

•	The first problem concerns the early expansion history of the current Hubble patch and its physical properties: Is the conventional timeline of the $\mathrm{\Lambda }$CDM scenario really compelling or just plausible?
•	The second class of questions involves the way relic gravitons could be used as a diagnostic of the early expansion history: How sensitive is the spectral energy density of the relic gravitons on the early expansion rates deviating from the $\mathrm{\Lambda }$CDM timeline?

To address the first group of subjects we should first acknowledge that the causal structure of FRW models provides already a number of relevant constraints on the expansion history. However, even admitting that, at early times, the particle horizon should disappear or diverge (as it happens in the case of conventional inflationary scenarios) to be replaced by an event horizon, the subsequent evolution of the space–time curvature remains undetermined. To appreciate this relevant point we should actually observe that the total number of e-folds does depend on the post-inflationary rate of expansion. For instance when we say that 60 e-folds of accelerated expansion are necessary to suppress the spatial curvature we are actually referring to a post-inflationary evolution dominated by radiation. The same tacit assumption is systematically employed to confront the temperature and the polarization anisotropies of the CMB with the conventional inflationary scenarios.^42,43,44

1.4. The relic gravitons and the expansion history

One of the purposes of this paper is to argue that the spectra of relic gravitons provide the only direct probe of the post-inflationary evolution prior to the formation of light nuclei. This is why a detailed analysis of such a signal is mandatory even in the absence of sensitive detectors that might be available only in the far future. Various secondary effects may produce different backgrounds of gravitational radiation during a fixed post-inflationary evolution like the one endorsed in the context of the $\mathrm{\Lambda }$CDM scenario. These effects, however, always assume a specific knowledge that is still missing. Conversely the relic gravitons do represent the only conceivable direct diagnostic of the post-inflationary expansion history and this is the general perspective developed here. Since the spectrum of the relic gravitons extends from the aHz region up to the THz domain we can partition this broad frequency domain into three complementary ranges where different stages of the early expansion rate are correspondingly probed:

•	the low-frequency region (between few aHz and the fHz) is directly sensitive to the expansion rate during inflation; in this region the upper limits on the tensor to scalar ratio deduced from the temperature and polarization anisotropies of the CMB are in fact bounds on the early expansion rate; the CMB can be in fact considered as the largest electromagnetic detector of long-wavelength gravitational waves;
•	at intermediate frequencies various potential constraints are associated with the PTAs (typically operating in the nHz domain); from the viewpoint of the expansion history this region may set constraints both on the post-inflationary evolution and on the modifications introduced during the inflationary stage;
•	finally in the high-frequency domain the constraints from the operating wide-band detectors between few Hz and 10kHz (as well as from other electromagnetic detectors operating in the MHz or GHz region) will be essential for the analysis of potential peaks in the spectrum of relic gravitons.

The first speculations suggesting that the relic gravitons could be used as a direct probe of the post-inflationary expansion history goes back to the late 1990s and this will be the general inspiration of this paper. In particular in Ref. 45 it has been suggested that different post-inflationary stages modify the slopes of the spectral energy density of the relic gravitons for frequencies larger than the mHz. It was found, quite surprisingly, that when the expansion rate is slower than radiation the spectral energy density exhibits a high-frequency spike.^46,47 The original observation of Ref. 45 was that the post-inflationary evolution may be modified and this would be especially true if we have to accommodate a late-time dominance of the dark energy. In this case a post-inflationary evolution dominated by radiation would be less likely than a long stiff stage expanding slower than radiation.⁴⁵ One of the first frameworks where these observations have been applied are the quintessential inflationary models.⁴⁸ In this context the late-time dominance of dark energy occurs via a quintessence field that ultimately coincides with the inflaton. Later on different scenarios based on different premised have been proposed⁴⁹ with the aim of accommodating an intermediate stage expanding at a rate different from radiation. For the purpose of this paper, however, we do not want to commit ourselves to a specific scenario or to a specific class of models. Indeed, as suggested in Ref. 45, the spectra of the relic gravitons chiefly depend on the evolution of the space–time curvature and not on the particular features involving the different sources.

1.5. The layout of this paper

The interplay between the timeline of the expansion rate and the spectra of the relic gravitons promises a direct connection between cosmology, quantum field theory and the effective description of gravitational interactions. On a more practical ground, in this topical investigation astrophysics and gravitational wave astronomy are seen as a tool for high-energy physics. In the past the common wisdom suggested instead that high-energy physics was probably the sole tool to infer properties of the primeval plasma prior to BBN. This conventional viewpoint did rest on the assumption that the post-inflationary expansion rate had to be fixed and almost perpetually dominated by radiation down to the scale of matter–radiation equality. In our context the timeline of the post-inflationary expansion rate is only a working hypothesis subjected to the direct tests associated with the diffuse backgrounds of gravitational radiation. Given the wealth of the connections between the various aspects of the problem it is impossible to analyze in detail all the relevant themes and this is why various collateral topics are swiftly mentioned but the interested readers may usefully consult a recently published book that dwells on the physics of the relic gravitons⁴⁹ where most of the considerations omitted here are systematically addressed. The layout of this paper is, in short, the following. Before elaborating on the unknowns, Sec. 2 is focused on what it is understood about the early expansion history with the goal of distinguishing the facts from the tacit assumptions. Section 3 deals more directly with the interplay between the relic gravitons and the expansion history and since the various ranges of the spectra are directly sensitive to the evolutionary stages of the background geometry, it seems useful to examine separately the interplay between the relic gravitons and the expansion histories in the low (see Sec. 4), intermediate (see Sec. 5) and high frequency (see Sec. 6). In Sec. 4, we point out that the inflationary observables are either suppressed or enhanced depending upon the post-inflationary evolution that affects the total number of e-folds. In Sec. 5, we present a discussion on the mutual interplay between the modified expansion histories and the PTAs; toward the end of Sec. 5 we also argue that a the post-inflationary evolution may also produce signals between few $\mu$Hz and the Hz where usually different sources are claimed to be relevant for the (futuristic) space-borne detectors. Finally in Sec. 6, we specifically address the direct bounds on the post-inflationary expansion rate coming from the high frequency and ultra-high frequency regions where absolute bounds on the maximal frequency of the spectra can be derived. Some ideas related to the use of quantum sensing for the detection of the relic gravitons will also be analyzed. The obtained limits on the maximal frequency are deeply rooted in the quantumness of the produced gravitons whose multiparticle final sates are macroscopic but always nonclassical. As the unitary evolution preserves their coherence, the quantumness of the gravitons can be associated with an entanglement entropy that is related with the loss of the complete information on the underlying quantum field. In the appendices we elaborated on some of the technical aspects that are often recalled in the main discussions. In particular, Appendix A illustrates a number of relevant complements on the evolution of curvature inhomogeneities that are specifically needed in the discussion while Appendix B treats the forms of the action of the relic gravitons in different frames.

2. The Timeline of the Expansion Rate: Facts and Tacit Assumptions

In the last 50 years the interplay between high-energy physics, astrophysics and cosmology has been guided by the tacit assumption that prior to matter–radiation equality the primeval plasma was always dominated by radiation^50,51 and this general truism is also reflected in various cartoons that are customarily employed to represent the timeline of the expansion rate where different moments of the life of the Universe are illustrated with the supposed matter content of the plasma. This viewpoint has been also propounded by Weinberg in one of the first popular accounts of the subject.⁵² After the formulation of the inflationary paradigm in its different variants (see e.g. Refs. 18, 19, 20, 21) the hypothesis of a post-inflationary radiation dominance remained practically unmodified and even today it is customary to assume that after an explosive stage of reheating the Universe should become, almost suddenly, dominated by radiation (see, for instance, Refs. 53, 54, 55). Among the various conclusions that emanate from the assumption of an evolution dominated by radiation, the most notable one is probably that the plasma as a whole is described by a single temperature for most of its history. An equally relevant statement is that the inflationary expansion must (or should) last for at least 60 e-folds.^53,54,55 Since this tacit assumption of radiation dominance is not directly tested (at least for temperatures larger than few MeV) more general possibilities will be discussed.

2.1. What do we know about the early expansion history?

2.1.1. Particle horizon and causally disconnected regions

A relevant constraint on the early expansion history comes from the causal structure of FRW models whose line element in its canonical form is given by

where $g_{\mu \nu }$ denotes the metric tensor and $a(t)$ is the scale factor. In the parametrization of Eq. (2.1), $\kappa =0$ corresponds to a spatially flat Universe; if $\kappa >0$ the Universe is spatially closed and, finally, $\kappa <0$ describes an open spatial section. In Eq. (2.1), the time t indicates the cosmic time coordinate but depending upon the physical problem at hand, different time parametrizations can be also adopted. A particularly useful one is the so-called conformal time parametrization that turns out to be particularly useful in the analysis of the inhomogeneities (see, in this respect, Appendix A and B). In the conformal time coordinate $\tau$ the line element of Eq. (2.1) becomes Since light rays follow null geodesics in Eq. (2.1) we may suppose that a signal is emitted at the time $t_{\mathrm{\text{e}}}$ (at a radial position $r_{\mathrm{\text{e}}}$) and received at the time $t_{\mathrm{\text{r}}}$ (at a radial position $r_{\mathrm{\text{r}}}=0$). Then from Eq. (2.1) with $d{s}^2=0$ and $d\vartheta =d\phi =0$ we will have, for a null radial geodesic The position of the emitter is fixed in the comoving coordinate system. We can then say that the signal was emitted at a physical distance $d(t)$ from the origin : If we now introduce the concept of $t_{min}$ (corresponding to the maximal past extension of the time coordinate on the FRW space–time), in the limit $t_e\to t_{min}$ we can introduce the particle horizon at time t as In short, for any given observer the particle horizon divides the regions of the space–time already observed from the ones that have not been observed yet. In the hot big bang scenario the background expands but it is simultaneously decelerated and this means that the first derivative of the scale factor $a(t)$ with respect to the (cosmic) time coordinate t is positive while its second derivative gets negative : where the overdot denotes here a derivation with respect to the cosmic time coordinate t and $H=ȧ∕a$ indicates the usual Hubble rate. If the scale factor is parametrized in terms of a power law $a(t)\simeq a_1{(t∕t_1)}^{\alpha }$ the conditions (2.6) imply $0<\alpha <1$ and this means that the particle horizon exists and it is finite In the limit $t_{min}\to 0$ (when t remains finite) $d_p(t)$ coincides (up to an irrelevant multiplicative constant factor) with $H^{-1}(t)\simeq t$. The particle horizon of Eq. (2.7) measures the extension of causally connected regions at time t and its finiteness poses a problem if the Universe always expands in a decelerated manner. Since during a decelerated stage of expansion the extension of a causal patch is of the order of $d_p(t)\simeq t$, the Hubble radius at any time preceding the current epoch must contain a finite number of causally disconnected regions. If we indicate with $H_0^{-1}$ the Hubble rate at the present time, at equality $H_0^{-1}(a_{\mathrm{\text{eq}}}∕a_0)<H_0^{-1}$. In other words, the Hubble patch at equality is comparatively smaller than today but the typical size of causally connected regions at equality is $d_p(t_{\mathrm{\text{eq}}})\sim t_{\mathrm{\text{eq}}}$, as suggested by Eq. (2.7). If we then measure $H_0^{-1}(a_{\mathrm{\text{eq}}}∕a_0)$ in units of $t_{\mathrm{\text{eq}}}$ we obtain As $d_p(t_{\mathrm{\text{eq}}})$ measures the extension of causally connected domains at $t_{\mathrm{\text{eq}}}$, Eq. (2.8) suggests that, at matter–radiation equality, the region corresponding to the present Hubble patch contained about 50 causally disconnected regions or $50^3\simeq 10^5$ disconnected volumes. The reference time selected in Eq. (2.8) can be modified but the essence of the problem remains the same. Furthermore if the typical reference time is larger than $t_{\mathrm{\text{eq}}}$ the number of causally disconnected regions is comparatively smaller. Conversely, when $t<t_{\mathrm{\text{eq}}}$ the number of the disconnected regions increase and quickly approaches its Planckian limit.^a

2.1.2. Event horizon

The previous discussion clarifies why the existence of the particle horizon leads necessarily to causally disconnected volumes; this occurrence clashes, among other things, with the high degree of homogeneity and isotropy of the Universe as it follows, for instance, from the analysis of the temperature and polarization anisotropies of the CMB. How come that regions emitting a highly homogeneous and isotropic CMB were causally disconnected in the past? To solve the causality problems of the conventional big bang scenario the idea is then to complement the standard decelerated stage of Eq. (2.6) with an epoch where the scale factor accelerates

and the particle horizon diverges. If the scale factor is parametrized with a power law $a(t)\simeq a_1{(t∕t_1)}^{\alpha }$ the conditions (2.9) demand that $\alpha >1$ and, in this situation, the particle horizon does not exist: when $\alpha >1$ the integral of Eq. (2.7) is divergent in the limit $t_{min}\to 0$ and for any finite value of the cosmic time coordinate t. If the Universe expands as in Eq. (2.9) there exist however an event horizon where the second approximate equality holds for t finite and $t_{max}\to +\infty$. In Eq. (2.7) $t_{min}$ measured the maximal past extension of the time coordinate in the given FRW space–time; in the case of the event horizon $t_{max}$ measures instead the maximal future extension of the cosmic time coordinate. For this reason the event horizon measures the maximal distance over which we can admit, even in the future, a causal connection. If $d_e(t)$ is finite in the limit $t_{max}\to \infty$ (for finite t) we can conclude that the event horizon exist. When the phase of accelerated expansion is parametrized in terms of the (expanding) branch of four-dimensional de Sitter space–time, namely, $a(t)\simeq e^{H_{i}t}$ (with $H_i>0$) the particle and event horizons are, respectively,

$$d_p(t)=H_i^{-1}[e^{H_i(t-t_{min})}-1],$$(2.11)

$$d_e(t)=H_i^{-1}[1-e^{H_i(t-t_{max})}].$$(2.12)

The cosmic time coordinate is allowed to run from^b $t_{min}\to -\infty$ up to $t_{max}\to +\infty$. Consequently, for $t_{\mathrm{\text{min}}}\to -\infty$ (at fixed t) the particle horizon will diverge and the typical size of causally connected regions at time t scales as $L_i(t)\simeq H_i^{-1}a(t)∕a(t_{min})$. While for the standard decelerated expansion the particle horizon increases faster than the scale factor, the typical size of causally connected regions scales exactly as the scale factor. In the limit $t_{max}\to \infty$ the event horizon exist and it is given, from Eq. (2.12), by $d_e(t)\simeq H_i^{-1}$.

2.1.3. Total number of e-folds?

If the accelerated stage of expansion is sufficiently long, all the scales that were inside $H^{-1}$ at the onset of inflation are today comparable (or larger) than the Hubble radius. It is essential to appreciate that the quantitative meaning of the locution sufficiently long depends also on the post-inflationary evolution and not only on the inflationary dynamics itself. The duration of the accelerated stage of expansion is customarily parametrized in terms of the ratio between the scale factors at the end (i.e. $a_f$) and at the beginning (i.e. $a_i$) of inflation^c :

where N denotes the number of e-folds. Later on in this section we shall be introducing with $N_k$ (namely, the number of e-folds elapsed since a given scale crossed the Hubble radius) as well as other notions derived from Eq. (2.13); the notion of $N_k$ becomes particularly relevant for the analysis of Sec. 4. Equation (2.14) is clearly equivalent to

The number of e-folds required for the consistency of a given inflationary scenario does not only depend on the inflationary dynamics as it might seem to follow from Eqs. (2.13) and (2.14). In other words, while the physical features of the decelerated and of the accelerated expansions are per se relevant, what we want to stress here is that the indetermination of the post-inflationary evolution affects the specific value of the total number of e-folds. To clarify this point we consider the ratio between the intrinsic (spatial) and the extrinsic (Hubble) curvatures and recall that it is notoriously given by

The right-hand side of Eq. (2.15) gets suppressed when the Universe accelerates (see Eq. (2.9)) while $\kappa ∕{ȧ}^2$ increases during a stage where $ä<0$ (see Eq. (2.6)). The total suppression of the ratio given in Eq. (2.15) cannot be simply attributed to the inflationary stage of expansion unless we artificially assume that the post-inflationary evolution is known and corresponds, for instance, to the dominance of radiation. This is ultimately the reason why the total number of e-folds suffers a theoretical indetermination associated with the post-inflationary evolution. In Fig. 1, we illustrate with a cartoon the sensitivity of the number of e-folds to the post-inflationary evolution. In the rightmost part of the plot we have the current value of the Hubble radius and the thick line denotes the standard evolution where the post-inflationary expansion rate is always dominated by radiation. If the blueshifted value of the current Hubble radius must fit exactly inside the inflationary event horizon, the value of N becomes, as we shall see, $\mathcal{𝒪}(60)$. In Fig. 1, two qualitatively different possibilities are also mentioned: in the first case the post-inflationary expansion rate is faster than radiation (see the dashed line of Fig. 1) in the second case the post-inflationary expansion rate is slower than radiation (see the dot-dashed line of Fig. 1). When expansion rate is faster than radiation the number of e-folds required to fit the blueshifted Hubble radius inside the inflationary event horizon is comparatively smaller than N (i.e. $N_{\mathrm{\text{faster}}}<N$); the opposite is true in case the post-inflationary expansion rate is slower than radiation (i.e. $N_{\mathrm{\text{slower}}}>N$). The quantitative aspects of Fig. 1 are analyzed in Subsec. 2.3 after a discussion of the early expansion rate (see Subsec. 2.2) which is relevant also for the determination of the number of e-folds.

2.2. The early expansion rate

2.2.1. Conventional inflationary stages

The expansion history during the inflationary stage follows from the equations connecting the Hubble rate to the corresponding sources. The single-field inflationary models can be notoriously analyzed in terms of the following scalar–tensor action (see, for instance, Ref. 55)

where ${\ell }_P$ denotes the Planck length and the following notations will be used throughout the whole discussion : Equation (2.16) should be regarded as the first term of an effective description where the higher derivatives are suppressed by the negative powers of a large mass M associated with the fundamental theory that underlies the effective action. The leading corrections to Eq. (2.16) consist of all possible terms containing four space–time derivatives⁵⁶ and Eq. (2.16) itself can be studied in two complementary perspectives:

•	if we presume, by fiat, that the post-inflationary evolution is fixed and known then the only sensible question and the sole concern should be somehow to reconstruct the functional dependence of the potential;
•	conversely if post-inflationary evolution is unknown (or only partially known) it is less meaningful to aim at a reconstruction of the inflaton potential from the large-scale data since the number of e-folds ultimately depends on the post-inflationary evolution.

For different reasons both approaches are appealing but the former corresponds to the conventional lore while the latter perspective is pursued in this discussion: the ultimate goal would be to test the post-inflationary expansion rate rather than arbitrarily postulating a specific timeline. The post-inflationary evolution is not going to be fixed and this choice has a specific impact on the remaining part of the discussion. For this reason the properties of the expansion rate during inflation are technically more essential than the form of the potential although the obtained results can be related (at any step) to the more conventional approach. In the single-field case (and for the background geometry of Eqs. (2.1) and (2.2)) the evolution equations follow from Eq. (2.16) and they can be written as

where, as previously remarked, the overdot denotes a derivation with respect to the cosmic time coordinate t. If Eqs. (2.18) and (2.19) are combined we obtain $2Ḣ=-{\ell }_P^2{\dot{\phi }}^2+2\kappa ∕a^2$; this is, in practice, the explicit form of the Raychaudhuri equation⁴⁹ written in the case of a scalar field source. In a stage where the decrease of the Hubble rate is sufficiently slow (i.e. $Ḣ\ll -H^2$) Eqs. (2.18) and (2.19) can be approximated as where the contribution of the spatial curvature is also neglected since it is sharply suppressed during an accelerated stage of expansion. In connection with Eqs. (2.19) and (2.20) the last technical remark concerns the cosmic time parametrization that can be traded for the conformal time coordinate defined as $a(\tau )d\tau =dt$. In the conformal time coordinate the expansion rate and its derivative are defined as with the prime now denoting a derivation with respect to $\tau$. Both parametrizations of the time coordinate will be used interchangeably; in the analysis of the effective action of the tensor modes of the geometry (see Appendix B) the conformal time parametrization turns out to be more convenient, as we are going to see.

2.2.2. The early expansion rate

Although the post-inflationary expansion rate modifies the number of e-folds (and consequently all the inflationary observables), the early expansion rate can be estimated, at least approximately, without a detailed knowledge of the post-inflationary evolution. This happens since the early expansion rate ultimately follows from the analysis of the spectrum of curvature inhomogeneities associated with the CMB scales that left the Hubble radius during the first stages of inflation and reentered before matter radiation equality. Since the curvature inhomogeneities are conserved when they evolve for scales larger than the Hubble radius, the early expansion rate does not depend upon the total number of e-folds and the rationale for this statement is illustrated in Fig. 2 where the common logarithm of $aH$ is reported as a function of the common logarithm of the scale factor. While during inflation $aH\propto a$, in a radiation-dominated stage $aH\propto a^{-1}$; the two ellipses of Fig. 2 parametrize the unknowns of the intermediate evolution but a detailed knowledge of that regime is not strictly necessary to set initial conditions for the temperature and for the polarization anisotropies. The CMB observations involve in fact a bunch of wavenumbers $k=\mathcal{𝒪}(k_p)$ where $k_p=0.002{\mathrm{\text{Mpc}}}^{-1}$ is the conventional pivot scale that is used to normalize the large-scale power spectra. These typical scales are pictorially indicated in the lower part of Fig. 2 where the two filled squares denote the moment where $k=\mathcal{𝒪}(k_p)$ gets of the order of $aH$. While the first crossing time occurs during inflation, the second one takes place prior to matter–radiation equality (see the right part of the cartoon). The scales $k=\mathcal{𝒪}(k_p)$ become again of the order of $aH$ when the Universe is already dominated by a radiation plasma (i.e. before matter–radiation equality) and their evolution is not affected by the unknowns of the post-inflationary evolution that may however modify the spectra at smaller scales; in this case the reentry of the fluctuations might not take place when the plasma is dominated by radiation.

2.2.3. Adiabatic and nonadiabatic solutions

The argument of Fig. 2 holds only under the hypothesis that curvature inhomogeneities are conserved in the limit $k<aH$, i.e. for typical wavelengths larger than the Hubble radius. This is exactly what happens when the evolution of the curvature inhomogeneities on comoving orthogonal hypersurfaces (conventionally denoted by $\mathcal{ℛ}$) is analyzed in the limit $k<aH$ (or $k\tau <1$). A complementary possibility is to employ $\zeta$ which measures the curvature inhomogeneities on the hypersurfaces where the density contrast is constant (see, for instance, Ref. 54 and references therein). Although $\mathcal{ℛ}$ and $\zeta$ are different variables, the Hamiltonian constraint associated with the relativistic fluctuations of the geometry stipulates that

where $\mathrm{\Psi }$ is the gauge-invariant generalization of the Newtonian potential (the so-called Bardeen potential⁵⁷) while $(p_t+{\rho }_t)$ is the total enthalpy density of the sources. According to Eq. (2.22), $\mathcal{ℛ}$ and $\zeta$ must obey the same evolution for $k\ll aH$ : In Eq. (2.23), $\delta p_{\mathrm{\text{nad}}}$ indicates the nonadiabatic pressure fluctuation^d which may arise if different barotropic fluids are simultaneously present.^53,54,55 In the single-field case $\delta p_{\mathrm{\text{nad}}}=0$ and heeding observations the temperature and polarization anisotropies of the CMB are consistent with adiabatic and Gaussian initial conditions.^42,43,44 The true question to ask in connection with Eq. (2.23) is how many adiabatic solutions and how many nonadiabatic solutions are compatible, for instance, with the conservation of $\mathcal{ℛ}$ and $\zeta$. This question, usually approached within the separate Universe picture⁵⁸ (see also Ref. 59 for a reintroduction of some of arguments given in Ref. 58), has been addressed by Weinberg in a series of papers.^60,61,62 In short the argument suggested in Refs. 60, 61, 62 stipulates that the evolution equations of the relativistic fluctuations of the geometry have always a pair of physical solutions for which $\delta p_{\mathrm{\text{nad}}}\to 0$ and $\mathcal{ℛ}$ approaches a constant for $k\tau \to 0$. In other words, following the usual terminology, there will be always at least a pair of adiabatic solutions with $\delta p_{\mathrm{\text{nad}}}=0$ and $\mathcal{ℛ}$ constant in the limit $k\tau \to 0$: one solution with $\mathcal{ℛ}\ne 0$; the other with $\mathcal{ℛ}=0$. Although this analysis is most easily performed in the conformally Newtonian gauge^57,63 the result is in fact gauge-invariant since $\mathcal{ℛ}$ is itself gauge-invariant. All in all we have that the adiabatic modes are ultimately conserved and this is why the early expansion rate does not depend on the post-inflationary evolution. Therefore, if the only source of large-scale inhomogeneity are the scalar and tensor modes of the geometry excited during the inflationary stage an estimate of the early expansion rate follows from the amplitude of the curvature inhomogeneities assigned for typical scales $k=\mathcal{𝒪}(k_p)$ as where $n_s$ and $n_T$ are, respectively, the scalar and the tensor spectral indices while ${\mathcal{𝒜}}_{\mathcal{ℛ}}$ and ${\mathcal{𝒜}}_T$ are the corresponding amplitudes at the scale $k_p$. The ratio between ${\mathcal{𝒜}}_T$ and ${\mathcal{𝒜}}_{\mathcal{ℛ}}$ is the tensor to scalar ratio and it is customarily denoted by $r_T$. The values of $r_T$, ${\mathcal{𝒜}}_{\mathcal{ℛ}}$ and $n_T$, in various different combinations, determine the properties of the expansion rate during inflation. In the context of single-field inflationary models these quantities are related via the so-called consistency conditions, as we are now going to discuss in further detail.

2.2.4. The scale-dependence of the expansion rate

In the absence of nonadiabatic contributions the evolution of the curvature inhomogeneities obeys a source-free evolution equation that can be written in a decoupled form (see Appendix A and discussion therein); since the inflationary bound on the expansion rate follows from the large-scale evolution of $\mathcal{ℛ}$, it is practical to recall the equation obeyed by the corresponding Fourier amplitudes :

During an inflationary stage of expansion the evolution of $z_{\phi }$ is proportional to the scale factor a via the (time-dependent) expression of the slow-roll parameter $ϵ(\tau )$ : The role of $ϵ(\tau )$ is here the most relevant since it measures the progressive suppression of the Hubble rate during the inflationary stage of expansion. Recalling now Fig. 2, Eq. (2.25) may be approximately^e solved by iteration when $k\ll aH$ : The first term of Eq. (2.27) is the constant adiabatic mode which obeys ${\mathcal{ℛ}}^{\prime }\simeq 0$ while the second term vanishes asymptotically and corresponds to the second adiabatic solution with $\mathcal{ℛ}\to 0$ (see also Eq. (2.22) and discussion thereafter); finally the third term of Eq. (2.27) vanishes exactly in the large-scale limit, i.e. for $k\to 0$. The form of the large-scale solution given by Eq. (2.27) is actually a concrete example of the general argument suggested in Refs. 60, 61, 62. Let us now go back to Fig. 2 and focus on the wavenumbers $k=\mathcal{𝒪}(k_p)$ that are larger than $aH$ (or which is the same $k^2\gg |z_{\phi }^{\prime \prime }∕z_{\phi }|$); in this case ${\mathcal{ℛ}}_k(\tau )=q_k(\tau )∕z_{\phi }(\tau )$ where $q_k(\tau )=e^{-ik\tau }∕\sqrt{2k}$. Since the solution must be continuous and differentiable in ${\tau }_{\mathrm{\text{ex}}}$ we can compute the approximate form of the scalar power spectrum From Eq. (2.28) recalling that $|q_k(\tau ){|}^2={(2k)}^{-1}$ we obtain where the expression of $z_{\phi }$ has been made explicit. The overall normalization of the scalar power spectrum is determined by the expansion rate at the typical time ${\tau }_{\ast }\simeq 1∕k$; thanks to Eq. (2.20) we can use $3H^2{\overline{M}}_P^2=V$ and $2{\overline{M}}_P^2Ḣ=-{\dot{\phi }}^2$ to simplify Eq. (2.29) : where, as before, $ϵ(\tau )=-Ḣ∕H^2$. During a de Sitter stage of expansion the scale factor is approximately given by $a(\tau )={(-H\tau )}^{-1}$ so that Eq. (2.30) ultimately becomes For a typical time ${\tau }_{\ast }=1∕k$ we then obtain the power spectrum that should be directly compared with the first of the two parametrizations of Eq. (2.24) where $H_{\ast }=H_k$ and ${ϵ}_{\ast }={ϵ}_k$ are evaluated for ${\tau }_{\ast }=1∕k$. Thus, after comparing Eqs. (2.32) and (2.24) for $k=\mathcal{𝒪}(k_p)$ we obtain the wanted estimate of $H_k$ : The condition (2.33) determines the expansion rate in Planck units which is then given by The second equality of Eq. (2.34) follows from the consistency relations stipulating that $r_T(k)\simeq 16{ϵ}_k$; this condition is typical of single-field inflationary scenarios (see Appendix A and discussion therein) and in Eq. (2.34) we adopted the notation^f $r_T=r_T(k_p)$. If we assume that $r_T\le 0.03$ and ${\mathcal{𝒜}}_{\mathcal{ℛ}}=\mathcal{𝒪}(10^{-9})$, Eq. (2.34) implies $H_k∕M_P\ll 1$ and since H decreases very little during inflation, the expansion rate few e-folds before the end of inflation is also comparable with $H_k$, i.e. $H_k=\mathcal{𝒪}(H)$. The ratio $H∕H_k$ may be estimated from the condition that defines the crossing of a given scale, i.e. $a_k{H}_k=k$ But for typical wavenumbers $k=\mathcal{𝒪}(k_p)$ it turns out that $k_p\ll |a_f{H}_f|$; more specifically the approximate value of $k_p∕(a_f{H}_f)$ is estimated as where $H_f\simeq H$. If we now insert Eq. (2.36) into Eq. (2.35), we can conclude, as previously anticipated that $H_k\simeq H$ so that Equation (2.37) estimates the expansion rate but does not imply any specific duration of the inflationary phase. As we are going to see in the following subsection, the duration of the inflationary stage of expansion is ultimately related to the nature and to the rate of the post-inflationary evolution.

2.3. What do we know about the late expansion history?

As already discussed after Eq. (2.13), the duration of inflation does depend on the post-inflationary evolution and this means that different expansion histories affect the number of e-folds required to bring all the physical scales of the model in causal contact. Different possibilities are examined hereunder with the purpose of quantifying the theoretical indetermination on the total number of e-folds.

2.3.1. A radiation-dominated universe?

One of the standard (unproven) assumptions both of the hot big bang model and of the conventional post-inflationary evolution is that the plasma must always be dominated by radiation even before the scale of BBN where the deviations from radiation dominance are severely constrained (see Sec. 4 and discussion therein). The gist of this argument is that, in the early hot and dense plasma, it is appropriate to assume an equation of state corresponding to a gas of relativistic particles; this choice is compatible with all the current data but it is neither compelling nor unique. A radiation-dominated stage of expansion extending between $H=\mathcal{𝒪}(10^{-5})M_P$ and the equality time is one of the assumptions customarily adopted for the timeline of the expansion rate in the context $\mathrm{\Lambda }$CDM paradigm. In terms of the cartoon of Fig. 2 we would then have that $aH\propto a$ during inflation while in the radiation stage $aH\propto a^{-1}$. This means, in practice, that the energy density of the background scales approximately as $a^{-4}$ between the end of inflation and the equality time. The critical number of e-folds required to fit inside the current Hubble patch the redshifted value of $H^{-1}$ (i.e. the approximate size of the event horizon at the onset of inflation) follows from the condition :

where $a_0$ is the current value of the scale factor.^g Equation (2.38) can be made even more explicit by rewriting it in a slightly different manner : The terms appearing in the second equality of Eq. (2.39) can be directly evaluated when the post-inflationary evolution is dominated by radiation; for instance, by definition, $3H_{\mathrm{\text{eq}}}^2{\overline{M}}_P^2=2{\rho }_{M0}{(a_0∕a_{\mathrm{\text{eq}}})}^3$ where ${\rho }_{M0}$ denotes the present matter density and the factor 2 follows since, at equality, the matter and radiation energy density coincide; furthermore the redshift to equality can be estimated as $(a_0∕a_{\mathrm{\text{eq}}})={\mathrm{\Omega }}_{M0}∕{\mathrm{\Omega }}_{R0}$ where ${\mathrm{\Omega }}_{M0}$ and ${\mathrm{\Omega }}_{R0}$ are the critical fractions of matter and radiation in the concordance scenario. All in all we can eventually estimate Moreover, between the equality time and $a_r$ the evolution is dominated by radiation, thanks to Eq. (2.40), Eq. (2.39) becomes where $g_{\rho }$ denotes the number of effective relativistic degrees of freedom appearing in the energy density of the plasma while $g_s$ corresponds to the number of effective relativistic degrees of freedom of the entropy density. If the entropy density is conserved between the r-stage and the equality epoch we should have that $g_{s,r}{a}_r^3{T}_r^3=g_{s,\mathrm{\text{eq}}}{a}_{\mathrm{\text{eq}}}^3{T}_{\mathrm{\text{eq}}}^3$ and this observation affects the redshift between the two epochs^h : Because during inflation the Hubble rate is nearly constant (i.e. $H_f\simeq H_i\simeq H$), once Eqs. (2.41) and (2.42) are inserted into Eq. (2.39) the number of e-folds ${\overline{N}}_{max}=ln(a_f∕a_i)$ can be determined by requiring that Eqs. (2.38) and (2.39) are satisfied. The final result for ${\overline{N}}_{max}$ becomes where, for the sake of conciseness, we wrote $\mathcal{𝒞}(g_s,g_{\rho },{\tau }_r,{\tau }_{\mathrm{\text{eq}}})={(g_{\rho ,r}∕g_{\rho ,\mathrm{\text{eq}}})}^{1∕4}{(g_{s,\mathrm{\text{eq}}}∕g_{s,r})}^{1∕3}$. Once more Eq. (2.43) determines the critical number of e-folds necessary to fit the redshifted value of $H^{-1}$ inside $H_0^{-1}$, as postulated in Eq. (2.38). It should be stressed that ${\overline{N}}_{max}$ corresponds to the maximal number of e-folds currently accessible to large-scale observations. In other words the conditions (2.38) and (2.39) fix ${\overline{N}}_{max}$ by requiring that all the physical scales inside the inflationary (event) horizon are all contained inside the current Hubble patch $H_0^{-1}$. It is of course possible that the total number of e-folds exceeds ${\overline{N}}_{max}$ and this happens if we require The condition (2.44) implies that some of the scales originally contained inside the inflationary (event) horizon are today larger than the current value of the Hubble patch; in this case the causal connection is realized on a region possibly larger than $H_0^{-1}$. The overlines appearing both in $\overline{N}$ and ${\overline{N}}_{max}$ remind that the corresponding quantities have been deduced for a post-inflationary evolution dominated by radiation. From Eq. (2.43), the explicit value of ${\overline{N}}_{max}$ becomes and it is, as anticipated, $\mathcal{𝒪}(60)$. For the actual estimates relating Eqs. (2.43) and (2.45) the following three observations should be emphasized:

•	the inflationary expansion rate is estimated from the amplitude of the scalar power spectrum and, more specifically, from Eq. (2.37);
•	it is assumed that, in practice, there is no energy loss between the inflationary phase and the post-inflationary evolution (i.e. $H_r\simeq H$);
•	in the standard situation where $g_{s,r}=g_{\rho ,r}=106.75$ and $g_{s,\mathrm{\text{eq}}}=g_{\rho ,\mathrm{\text{eq}}}=3.94$ the value of $\mathcal{𝒞}(g_s,g_{\rho },{\tau }_r,{\tau }_{\mathrm{\text{eq}}})$ is given by $0.75$; the contribution of $\mathcal{𝒞}(g_s,g_{\rho },{\tau }_r,{\tau }_{\mathrm{\text{eq}}})$ to Eq. (2.45) is numerically not essential for the determination of ${\overline{N}}_{max}$.

The approximation $H_r\simeq H$ is customarily enforced by CMB experiments when setting bounds, for instance, on the total number of e-folds^42,43,44 and although energy is lost during reheating, in the case of single-field inflationary models this approximation is rather plausible since the combined action of the reheating and of the preheating dynamics leads to a process that is almost sudden.^53,54,55 In this sense, if $H_{\mathrm{\text{last}}}$ denotes the expansion rate during the last few e-folds of inflationary expansion, it is true that $H_r<H_{\mathrm{\text{last}}}$; however, even for a difference of few orders of magnitude the quantitative arguments illustrated here will not be crucially affected. We recall that, conventionally, the reheating is the period where the entropy observed in the present Universe is produced and it typically takes place when all the large-scale inhomogeneities of observational interest are outside the horizon. The different approaches to the reheating dynamics are not expected to affect the large-scale power spectra.⁶³

The number of inflationary e-folds introduced in Eqs. (2.13)–(2.15) depends on the post-inflationary evolution but it also scale-dependent. This happens because the actual observations always probe a typical scale so that this dependence also enters the number of e-folds and the expansion rate. In what follows $N_k$ and $H_k$ are associated, respectively, with the number of e-folds and with the expansion rate at the crossing of the CMB scales $k=\mathcal{𝒪}(k_p)$. Even though $H_k$ and H are conceptually different, $H_k∕H=\mathcal{𝒪}(1)$ since the curvature scale decreases very slowly during the inflationary stage. Most of the previous estimates can then be repeated in the case of ${\overline{N}}_k=ln(a_f∕a_k)$. As in the case of ${\overline{N}}_{max}$, also ${\overline{N}}_k$ is estimated in the present section for a post-inflationary thermal history dominated by a radiation background and this is why the overline is included; the values of ${\overline{N}}_k$ are implicitly determined from

and when the given wavenumber is of the order of the comoving expansion rate $k\simeq a_k{H}_k$. The latter condition fixes the value of ${\overline{N}}_k$ not only in terms of $H_f$ (the expansion rate at the end of inflation) but also as a function of the subsequent expansion history, exactly as in the case of ${\overline{N}}_{max}$. By then repeating all the different steps in the case of Eq. (2.46) we deduce If the determinations of ${\overline{N}}_k$ and ${\overline{N}}_{max}$ are compared in the case of a post-inflationary evolution dominated by radiation we obtain As already mentioned in Eq. (2.37), $H_k=\mathcal{𝒪}(H_f)$ so that ${\overline{N}}_k$ and ${\overline{N}}_{max}$ are of the same order as long as $k\simeq a_0{H}_0$. The explicit value of ${\overline{N}}_k$ can then be written asⁱ

2.3.2. An extra phase preceding big bang nucleosynthesis

In the previous subsection, we considered a timeline dominated by radiation between the end of inflation and the equality epoch. We are now going to suppose that, prior to radiation dominance, the expansion rate is modified for a sufficiently long period where the expansion rate can be either faster or slower than radiation. Probably the simplest example along this perspective consists in adding a further stage of expansion between the end of inflation and the onset of the radiation-dominated phase. The ellipses of Fig. 2 are now replaced by the cartoon of Fig. 3 and the following comments are in order:

•	as before during inflation we have that $Ha\propto a$ while in a radiation stage we would get $aH\propto a^{-1}$: the simplest timeline is then the one illustrated with the full thick line;
•	prior to the onset of the radiation stage and after inflation we have instead that $aH\propto a^{-1∕\delta }$ where now $\delta$ parametrizes the expansion rate in the intermediate regime;
•	if $\delta >1$ the expansion rate is faster than radiation; conversely when $\delta <1$ the expansion rate is slower than radiation (see, in this respect, the dashed timelines of Fig. 3).

According to Fig. 3, the condition imposed by Eq. (2.39) becomes different and its modification depends on $\delta$. Indeed, if the estimate of Eq. (2.39) is repeated, the value of ${\overline{N}}_{max}$ gets shifted^45,46,47 (see also Refs. 64 and 65)

where we now denote with $N_{max}$ the maximal number of e-folds for a generic post-inflationary evolution while ${\overline{N}}_{max}$ corresponds to the case of a timeline dominated by radiation right after the end of the inflationary expansion. This is why, as anticipated, in the case $\delta \to 1$ the timeline of Fig. 3 reproduces a (single) radiation-dominated stage of expansion and $N_{max}\to {\overline{N}}_{max}$ (see Eqs. (2.43)–(2.50) and discussion therein). Because $H_r<H<H_i$ in Eq. (2.50), for arbitrary values of $\delta$ the following two remarks are in order:

•	when the background expands faster than radiation (i.e. $\delta >1$) the value of $N_{max}$ gets smaller than in the case of radiation dominance (i.e. $N_{max}<{\overline{N}}_{max}$);
•	conversely when the expansion rate is slower than radiation (i.e. $\delta <1$) we have that $N_{max}>{\overline{N}}_{max}$.

The orders of magnitude involved in Eq. (2.50) are estimated by considering that the typical expansion scale of BBN is approximately $H_{\mathrm{\text{bbn}}}=\mathcal{𝒪}(10^{-44})M_P$ whereas the inflationary expansion rate follows from Eq. (2.37) (i.e. $H\simeq \sqrt{\pi {\mathcal{𝒜}}_{\mathcal{ℛ}}{r}_T}∕4$). This means that the relation between $N_{max}$ and ${\overline{N}}_{max}$ is approximately given by

Let us now suppose, for instance, that $\delta >1$. If the sources for the evolution of the geometry are parametrized in terms of perfect fluids with barotropic equation of state $\delta =2∕(3w+1)$ so that ${\delta }_{max}\to 2$ corresponds^j to $w\to w_{min}=0$. In this case, from Eqs. (2.50) and (2.51), $N_{max}={\overline{N}}_{max}-\mathcal{𝒪}(15)$. In case $\delta <1$ the post-inflationary expansion rate between H and $H_r$ is instead slower than radiation so that we would have $N_{max}>{\overline{N}}_{max}$. Again, assuming the background is driven by perfect barotropic fluids, ${\delta }_{min}=1∕2$ and it corresponds to a plasma $w_{max}=1$ where the sound speed and the speed of light coincide. Therefore for $\delta \to {\delta }_{min}=1∕2$ Eqs. (2.50) and (2.51) imply that $N_{max}={\overline{N}}_{max}+\mathcal{𝒪}(15)$. In summary the critical number of e-folds required to fit the redshifted event horizon inside the current value of the Hubble radius does depend on the post-inflationary expansion rate; thanks to the results of Eqs. (2.50) and (2.51) we can then estimate the theoretical error associated with the unknown post-inflationary expansion rate as

where, we remind, the value of ${\overline{N}}_{max}$ is determined in the case $\delta \to 1$ corresponding to a radiation-dominated stage of expansion. The same kind of evaluation leading to Eq. (2.52) can be repeated for different classes of sources driving the background evolution. For instance the post-inflationary expansion rate might correspond to a stage dominated by an oscillating scalar field with an approximate potential ${\phi }^{2q}$ near the origin⁶⁶ (see also Refs. 67, 68,69). In this case $\delta =(q+1)∕(2q-1)$ and the condition $\delta \ge 1$ implies that $q\le 2$; this means, once more, that ${\delta }_{max}=2$ while ${\delta }_{min}\to 1∕2$ corresponding either to the asymptote $q\gg 1$ or to the absence of the potential. Thus the case ${\delta }_{min}\to 1∕2$ may be realized in a number of physically different situations.⁷⁰

All in all, if the total number of e-folds is $\mathcal{𝒪}(60)$ in the case of a radiation-dominated Universe, Eq. (2.52) suggests that the potential indetermination due to a modified expansion rate ranges^k between 45 and 75. The same indetermination affecting $N_{max}$ also enters the value of $N_k$. Indeed even in the presence of an intermediate stage preceding the conventional radiation-dominated epoch Eqs. (2.47) and (2.48) remain fully valid. The value of $N_k$ is relevant for various phenomenological aspects of the problem since it affects the inflationary observables that are specifically discussed later on^l in Sec. 4. We finally recall that Eq. (2.36) has been correctly deduced in the case of radiation dominance (i.e. $\delta \to 1$ in the language of this subsection) and the same indetermination affecting the number of e-folds may also modify the value of the pivot scale in units of the inflationary expansion rate. In the presence of the $\delta$-phase illustrated in Fig. 3 we have that Eq. (2.36) gets modified as

Depending on the value of $H_r$, when the expansion rate is faster than radiation the value of $k_p∕(a_f{H}_f)$ may get smaller than $10^{-26}$. The opposite is true when the background expands at a rate slower than radiation since, in this second instance, $k_p∕(a_f{H}_f)$ gets larger than $10^{-26}$. In both situations, however, it is fully justified to assume $H_f\simeq H_k\simeq H$, as already established in Eq. (2.35).

2.3.3. Multiple stages preceding big bang nucleosynthesis

A natural extension of the results obtained in Eqs. (2.50)–(2.52) involves the presence of multiple post-inflationary stages parametrized by different values of the expansion rate conventionally denoted by ${\delta }_i$ with $i=1,\dots ,n$. It is actually plausible to generalize the previous considerations by replacing the single $\delta$ stage with n intermediate phases of expansion preceding the epoch of radiation dominance, as illustrated in Fig. 4. The cartoon of Fig. 3 is then substituted by the timeline of Fig. 4 where the initial stage of the post-inflationary evolution begins after the end of inflation (i.e. $H\simeq H_f=H_1$) while the nth stage conventionally coincides with the standard radiation-dominated evolution i.e. $a_r=a_n$ and ${\delta }_n=1$. As already explained before, we should always require $H_r>10^{-44}{M}_P$ implying that the BBN takes place when radiation is already dominant. During the ith stage of the sequence $aH\propto a^{-1∕{\delta }_i}$ and the expression of $N_{max}$ given in Eq. (2.50) can be generalized to the timeline of Fig. 4 :

The various ${\xi }_i$ appearing in Eq. (2.54) measure the duration of each post-inflationary stage of expansion and since the rate is always decreasing we may conclude that Since, by construction, $a_n=a_r$ we also have that $H_n=H_r$; this means that ${\xi }_{n-1}=H_n∕H_{n-1}=H_r∕H_{n-1}$. It finally follows from Fig. 4 that the product of all the ${\xi }_i$ coincides with $H_r∕H$, namely, This also means that if all the ${\delta }_i$ are equal the result of Eq. (2.54) coincides with the one of Eq. (2.50) obtained for a single $\delta$-phase. If the post-inflationary plasma is only dominated by radiation then in Eq. (2.54) all the ${\delta }_i$ go to 1 and the whole contribution disappears. Conversely when some of the ${\delta }_i$ are smaller than 1 both $N_{max}$ and $N_k$ increase. For ${\delta }_i>1$ we may have the opposite effect suggesting an overall reduction of $N_{max}$ and $N_k$. Both effects are relevant in low-frequency region of the relic graviton spectrum, as we are going to see more specifically in Sec. 4. Indeed the result of Eq. (2.48) remains valid also for the timeline of Fig. 4; this means that not only $N_{max}$ but also $N_k$ gets reduced or enhanced depending on the values of the various ${\delta }_i$, as it follows from the explicit expression of $N_k$ : Besides the case of Fig. 4 we can also take into account a further possibility that is illustrated in Fig. 5. Prior to a conventional stage of inflationary expansion there could be a stage where the expansion rate is different. This may happen for various reasons and, in the most conservative perspective, it could be that the evolution of the relic gravitons develops a refractive index even thought the dynamics of the background is always inflationary^m (see, in this respect, Appendix B and Sec. 5).

In the previous cartoons of this section, we illustrated the effective rate of expansion in Planck units even though, in various cases, it is also useful to reason in terms of the inverse of $aH$. For this reason in Fig. 6 we now plot ${(aH)}^{-1}$. Sometimes in the literature ${(aH)}^{-1}$ is referred to as the horizon or simply the Hubble radius. According to this terminology the different wavelengths of the gravitational waves and of the scalar modes of the geometry cross the Hubble radius at different times. The first crossing typically occurs during inflation (see the left part of Fig. 6); after the first crossing the wavelength gets larger than the Hubble radius. This moment is then referred to as the exit of the given wavelength. The second crossing (see the right part of Fig. 6) occurs in the decelerated stage of expansion and it is conventionally referred to as the reentry of the given wavelength since after this typical time the wavelength gets again smaller than the Hubble radius. The filled squares in Fig. 6 define the exit of a given (comoving) wavelength while the dots in the right portion of the plot denote reentry of the selected scale. According to Fig. 6, the wavelengths smaller than ${\lambda }_r$ reenter before radiation dominance while the wavelengths $\lambda >{\lambda }_r$ reenter between the onset of radiation dominance and the epoch of matter–radiation equality. For $\lambda <{\lambda }_r$ the wavelength $\mathcal{𝒪}({\lambda }_{min})$ corresponds to comoving frequencies close to the maximal (i.e. $\nu =\mathcal{𝒪}({\nu }_{max})$). The scales ${\lambda }_r<\lambda <{\lambda }_{\mathrm{\text{eq}}}$ were still larger than the comoving horizon prior to matter–radiation equality and exited about $N_k$ e-folds before the end of inflation; the corresponding wavenumbers range therefore between $0.05{\mathrm{\text{Mpc}}}^{-1}$ and $0.002{\mathrm{\text{Mpc}}}^{-1}$.

3. The Relic Gravitons and the Expansion History

During the last 50 years a recurrent viewpoint has been that, ultimately, high-energy physics is a tool for cosmology and astrophysics. This argument rests on the observation that the plasma became transparent to electromagnetic radiation only rather late (i.e. after the last scattering of photons). Therefore there cannot be direct signals coming, for instance, from an expanding stage with a typical temperature of the order of few TeV. However, since these energy scales are reachable by colliders, particle physics is the only tool that we might have to scrutinize the early Universe. This perspective (implicitly assuming the dominance of radiation and the existence of a prolonged stage of local thermal equilibrium) should be probably revamped in the light of the direct detection of gravitational radiation. Indeed we do know that every variation of the space–time curvature produces shots of relic gravitons with given multiplicities and specific spectra.⁴⁹ Since the sensitivities of gravitational wave detectors greatly improved in the last 30 years, it is plausible to assume that the direct observations might hit the thresholds of the cosmological signals during the next score year or so. Under this hypothesis the timeline of the expansion rate illustrated in the previous section may be one day testable in practice as it is already scrutinized in principle. Along this revamped perspective gravitational wave astronomy could become a tool for high-energy physics by conveying a more specific knowledge of energy scales that might not be accessible to colliders in the future.

The relic gravitational waves produced by the early variation of the space–time curvature^14,15,16,17 lead to a late-time background of diffuse radiation. In the simplest situation the relic gravitons are produced in pairs of opposite three-momenta from the inflationary vacuum and this is why they appear as a collection of standing (random) waves which are the tensor analog of the so-called Sakharov oscillations⁷¹; this phenomenon has been also independently discussed in the classic paper of Peebles and Yu⁷² (see also Ref. 73). The late-time properties of the signal not only rest on the features of the inflationary vacuum but also on the post-inflationary evolution. It is well established that in the concordance paradigm the spectral energy density at late times is quasi-flat^22,23,24 and it gets larger at smaller frequencies of the order of the aHz.²⁵ This happens because, in the concordance scenario, the spectral energy density scales as ${\nu }^{-2}$ between few aHz and 100aHz in the region where the current CMB observations are now setting stringent limits on the contribution of the relic gravitons to the temperature and polarization anisotropies.^42,43,44 Along this perspective the low-frequency constraints translate into direct bounds on the tensor to scalar ratio $r_T$ and seem to suggest that at higher frequencies (i.e. in the audio band and beyond) the spectral energy density in critical units should be $\mathcal{𝒪}(10^{-17})$ or even smaller. This result has been realized, at a different level of accuracy, in various papers starting from Refs. 22, 23, 24 (see also Refs. 74 and 75). The minuteness of the spectral energy density follows from the presumption that radiation dominates (almost) right after the end of inflation and it is otherwise invalid. As we argued in Sec. 2, the post-inflationary evolution prior to BBN is not probed by any direct observation and may deviate from the radiation-dominated timeline; if this is the case, the high-frequency spectrum of the relic gravitons can be much larger.^45,46,47 In what follows, we are going to discuss first the statistical properties of the gravitons produced by the variation of the space–time curvature; in the second part of the section the discussion is focused on the slopes of the spectral energy density of the relic gravitons and on their connection with the expansion rate of the Universe.

3.1. Random backgrounds and quantum correlations

The random backgrounds associated with the relic gravitons are homogeneous but not stationary and this property is ultimately related with their quantum mechanical origin. Conversely the homogeneity of the background does not directly follow from the properties of the quantum mechanical correlations. In what follows, we shall try to clarify the analogies and the differences between these two aspects of the problem by swiftly summarizing the main conclusions of a recent analysis⁷⁶ that follows previous attempts along similar directions.⁷⁷

3.1.1. The energy density of random backgrounds

We start by considering a tensor random field $h_{ij}(\overrightarrow{x},\tau )$ and its Fourier transformⁿ :

where $\overrightarrow{k}$ is the comoving three-momentum. The Fourier amplitude $h_{ij}(\overrightarrow{k},\tau )$ can be decomposed in terms of the tensor polarizations as where the sum over $\lambda$ runs over $\oplus$ and $\otimes$. If we introduce a triplet of mutually orthogonal unit vectors $\widehat{m}$, $\widehat{n}$ and $\widehat{k}$ (where $\widehat{m}\times \widehat{n}=\widehat{k}$) the two tensor polarizations are If background is isotropic and unpolarized the corresponding ensemble averages of the Fourier amplitudes (and of their first derivatives) can be expressed as

$$〈h_{ij}(\overrightarrow{k},\tau )h_{mn}(\overrightarrow{k^{\prime }},\tau )〉=\frac{2{\pi }^2}{k^3}{P}_T(k,\tau ){\delta }^{(3)}(\overrightarrow{k}+\overrightarrow{k^{\prime }}){\mathcal{𝒮}}_{ijmn}(\widehat{k}),$$(3.4)

$$〈{\partial }_{\tau }{h}_{ij}(\overrightarrow{k},\tau ){\partial }_{\tau }{h}_{mn}(\overrightarrow{k^{\prime }},\tau )〉=\frac{2{\pi }^2}{k^3}{Q}_T(k,\tau ){\delta }^{(3)}(\overrightarrow{k}+\overrightarrow{k^{\prime }}){\mathcal{𝒮}}_{ijmn}(\widehat{k}),$$(3.5)

where the two tensor power spectra $P_T(k,\tau )$ and $Q_T(k,\tau )$ fully describe the tensor random field; $〈\cdots 〉$ denotes an average over an ergodic ensemble of random functions. In Eq. (3.5), the tensor ${\mathcal{𝒮}}_{ijmn}(\widehat{k})$ arises from the sum over the two tensor polarizations : where $p_{ij}=({\delta }_{ij}-{\widehat{k}}_i{\widehat{k}}_j)$. From the $(00)$ component of the energy–momentum pseudo-tensor discussed in Appendix B (see in particular Eq. (B.16)) the energy density of the relic gravitons becomes If we now insert Eq. (3.1) inside Eq. (3.7) and average the obtained result according to Eqs. (3.4) and (3.5) we obtain In Eq. (3.8), ${\overline{\rho }}_{gw}=〈{\rho }_{gw}〉$ represents the ensemble average of the energy density; the second equality in Eq. (3.8) directly follows from Eq. (3.1) after taking the ensemble average of each term according to Eqs. (3.4) and (3.5). From Eq. (3.8), we can always introduce the spectral energy density in critical units : where, as before, ${\rho }_{\mathrm{\text{crit}}}=3{\overline{M}}_P^2{H}^2$. The value of ${\mathrm{\Omega }}_{gw}(k,\tau )$ depends both on $P_T(k,\tau )$ and $Q_T(k,\tau )$. Sometimes ${\mathrm{\Omega }}_{gw}(k,\tau )$ is swiftly referred to as the energy density (in critical units) of the random background but this terminology is incorrect: the energy density does not depend on the frequency (or on the momentum); ${\mathrm{\Omega }}_{gw}(k,\tau )$ represents the energy density (in critical units) and per logarithmic interval of momentum (or frequency) since, in our units, $\omega =k=2\pi \nu$. Since $Q_T(k,\tau )\to k^2{P}_T(k,\tau )$ for $k\gg aH$, the spectral energy density for typical wavelengths shorter that the Hubble radius can also be expressed as^o

3.1.2. Homogeneity in space

The results of Eqs. (3.9) and (3.10) follow by considering the basic features of traceless and solenoidal tensor random fields supplemented by the notion of stochastic average introduced in Eqs. (3.4) and (3.5). A relevant result following from the previous considerations is that the two-point function of the tensor modes is homogeneous in space. By this we mean that the two-point function only depends on the distance between two spatial locations. If we compute the correlation functions of $h_{ij}(\overrightarrow{x},\tau )$ and of its derivative at equal times (but for two different spatial locations) we obtain

where $j_0(kr)$ is the spherical Bessel function of zeroth order.^78,79 We remark that Eqs. (3.11) and (3.12) follow directly from the definition of Eq. (3.1) and from the averages of Eqs. (3.4) and (3.5). We note that the homogeneity in space implies that both correlators are evaluated at the same values of the conformal time coordinate $\tau$. The results of Eqs. (3.11) and (3.12) demonstrate that the tensor random fields, heuristically defined by Eqs. (3.1) and (3.4), (3.5), can be described by stochastic processes that are homogeneous in space. This means also that the two-point function computed at equal times (but for different locations) is invariant under spatial translations.

3.1.3. Homogeneity in time (stationarity)

It would now seem that the same kind of invariance should also hold when the spatial location is fixed but the time coordinates are shifted. In this case the two-point function of the tensor fluctuations would also be stationary, i.e. invariant under time translations. The stationarity is actually more restrictive than homogeneity if the random background is defined by Eqs. (3.1) and (3.4), (3.5). Indeed, as we are going to see, the stationarity ultimately restricts the time-dependence of the power spectra $P_T(k,\tau )$ and $Q_T(k,\tau )$. If we then avoid the complication of the spatial dependence and directly discuss a single tensor polarization $h(\tau )$, instead of an ensemble or random fields we deal an ensemble of real random functions $h(\tau )$. We then introduce the autocorrelation function ${\mathrm{\Gamma }}_h(\mathrm{\Delta }\tau )$ defined in the context of the generalized harmonic analysis and associated with the finiteness of the integral^80,81

Wiener considered the class of functions (all measurable in a Lebesgue sense) for which the integral (3.13) exists and demonstrated that the spectral density exists.⁸² In the case of a stationary and ergodic ensemble of random functions, the autocorrelation of Eq. (3.13) can be replaced by where $〈\cdots 〉$ now denotes an ensemble average and the results of Eqs. (3.13) and (3.14) must ultimately coincide under the hypotheses of ergodicity. The property expressed by Eq. (3.14) is characteristic of a stationary process whose autocorrelation function is invariant under a shift of the time coordinate. This is why the Fourier transform of the autocorrelation function is associated with a well-defined spectral amplitude.^82,83 Recalling Eq. (3.14) we can Fourier transform $h(\tau )$, obtain $h(\nu )$ as a function of the frequency $\nu$ and eventually evaluate the corresponding ensemble average; the result is From Eqs. (3.14) and (3.15), the autocorrelation function and the spectral amplitudes are then related as According to Eqs. (3.15) and (3.16), the spectral amplitude and the autocorrelation function of the process form a Fourier transform pair; this statement is often referred to as Wiener–Khintchine theorem (see e.g. Ref. 81) and was originally developed in the framework of the generalized harmonic analysis that establishes a rigorous connection between Eqs. (3.13) and (3.14).^82,83 The possibility of defining a spectral amplitude relies then on the stationary nature of the underlying random process. The spectral amplitude is actually measured in units of inverse frequencies and can also be assigned in the case of a generic spatial dependence. Both stationarity and homogeneity play an important role when analyzing the correlation between gravitational wave detectors of arbitrary geometry.^p84,85,86,87

3.2. Random backgrounds and quantum mechanics

For a quantum description of the relic gravitons the first step is to recall the second-order action for the tensor inhomogeneities deduced in Appendix B (see, in particular, Eq. (B.15)). The canonical momentum deduced from Eq. (B.15) is in fact given by ${\pi }_{ij}=a^2{\partial }_{\tau }{h}_{ij}∕(8{\ell }_P^2)$ and the resulting classical Hamiltonian is

By promoting the classical fields to the status of quantum operators (i.e. $h_{ij}(\overrightarrow{x},\tau )\to {ĥ}_{ij}(\overrightarrow{x},\tau )$ and ${\pi }_{ij}(\overrightarrow{x},\tau )\to {\widehat{\pi }}_{ij}(\overrightarrow{x},\tau )$) the quantum Hamiltonian ${Ĥ}_g(\tau )$ becomes where ${ĥ}_{ij}^{†}(\overrightarrow{x},\tau )={ĥ}_{ij}(\overrightarrow{x},\tau )$ and ${\widehat{\pi }}_{ij}^{†}(\overrightarrow{x},\tau )={\widehat{\pi }}_{ij}(\overrightarrow{x},\tau )$ are both Hermitian; the dagger denotes, as usual, the Hermitian conjugation. From Eq. (3.18), the evolution equations of the field operators in the Heisenberg description are and their explicit form in the Heisenberg representation is

$${ĥ}_{ij}(\overrightarrow{x},\tau )=\frac{\sqrt{2}{\ell }_P}{{(2\pi )}^{3∕2}}\sum _{\alpha =\oplus ,\otimes }\int d^{3}k{e}_{ij}^{(\alpha )}(\widehat{k})[F_{k,\alpha }(\tau ){\widehat{b}}_{\overrightarrow{k},\alpha }{e}^{-i\overrightarrow{k}\cdot \overrightarrow{x}}+\mathrm{\text{H.c.}}],$$(3.20)

$${\widehat{\pi }}_{ij}(\overrightarrow{x},\tau )=\frac{a^2(\tau )}{4\sqrt{2}{\ell }_P{(2\pi )}^{3∕2}}\sum _{\beta =\oplus ,\otimes }\int d^{3}k{e}_{ij}^{(\beta )}(\widehat{k})[G_{k,\beta }(\tau ){\widehat{b}}_{\overrightarrow{k},\beta }{e}^{-i\overrightarrow{k}\cdot \overrightarrow{x}}+\mathrm{\text{H.c.}}].$$(3.21)

In Eqs. (3.20) and (3.21), the second term inside the square bracket denotes the Hermitian conjugate of the preceding one. As before the sum runs over the two tensor polarizations defined in Eq. (3.3); because of Eq. (3.19) the mode functions $F_k(\tau )$ and $G_k(\tau )$ obey : where, as usual, $\mathcal{\mathscr{H}}=a^{\prime }∕a$. The Fourier transforms of the Hermitian field operators of Eqs. (3.20) and (3.21) are

$${ĥ}_{ij}(\overrightarrow{q},\tau )=\sqrt{2}{\ell }_P\sum _{\alpha }[e_{ij}^{(\alpha )}(\widehat{q}){\widehat{b}}_{\overrightarrow{q},\alpha }{F}_{q,\alpha }(\tau )+e_{ij}^{(\alpha )}(-\widehat{q}){\widehat{b}}_{-\overrightarrow{q},\alpha }^{†}{F}_{q,\alpha }^{\ast }(\tau )],$$(3.23)

$${\widehat{\pi }}_{mn}(\overrightarrow{p},\tau )=\frac{a^2}{4\sqrt{2}{\ell }_P}\sum _{\beta }[e_{mn}^{(\beta )}(\widehat{p}){\widehat{b}}_{\overrightarrow{p},\beta }{G}_{p,\beta }(\tau )+e_{mn}^{(\beta )}(-\widehat{p}){\widehat{b}}_{-\overrightarrow{p},\beta }^{†}{G}_{p,\beta }^{\ast }(\tau )].$$(3.24)

The field operators of Eqs. (3.23) and (3.24) obey the canonical commutation relations provided the mode functions obey the Wronskian normalization The condition expressed by Eq. (3.26) is essential to obtain the correct form of the commutation relations that must be preserved throughout all the stages of the dynamical evolution. The mode functions can also be rescaled as $F_k(\tau )=a{f}_k(\tau )$ and $G_k(\tau )=a{g}_k(\tau )$; in this case Eq. (3.26) becomes $f_k(\tau )g_k^{\ast }(\tau )-f_k^{\ast }(\tau )g_k(\tau )=i$.

3.2.1. Quantum mechanics and nonstationary processes

To analyze the stationarity of the process we need to introduce the autocorrelation functions depending on two different times ${\tau }_1$ and ${\tau }_2$:

The explicit form of Eqs. (3.27) and (3.28) follows from the actual expressions of the field operators in Fourier space (see Eqs. (3.23) and (3.24)) and the final result becomes

$${\mathrm{\Gamma }}_{ijmn}(\overrightarrow{k},\overrightarrow{p},{\tau }_1,{\tau }_2)={\mathcal{𝒮}}_{ijmn}(\widehat{k}){\delta }^{(3)}(\overrightarrow{k}+\overrightarrow{p}){\mathrm{\Delta }}_k({\tau }_1,{\tau }_2),$$(3.29)

$${\overline{\mathrm{\Gamma }}}_{ijmn}(\overrightarrow{k},\overrightarrow{p},{\tau }_1,{\tau }_2)={\mathcal{𝒮}}_{ijmn}(\widehat{k}){\delta }^{(3)}(\overrightarrow{k}+\overrightarrow{p}){\overline{\mathrm{\Delta }}}_k({\tau }_1,{\tau }_2),$$(3.30)

where ${\mathrm{\Delta }}_k({\tau }_1,{\tau }_2)$ and ${\overline{\mathrm{\Delta }}}_k({\tau }_1,{\tau }_2)$ are now given by

$${\mathrm{\Delta }}_k({\tau }_1,{\tau }_2)=4{\ell }_P^2[F_k({\tau }_1)F_k^{\ast }({\tau }_2)+F_k({\tau }_2)F_k^{\ast }({\tau }_1)],$$(3.31)

$${\overline{\mathrm{\Delta }}}_k({\tau }_1,{\tau }_2)=4{\ell }_P^2[G_k({\tau }_1)G_k^{\ast }({\tau }_2)+G_k({\tau }_2)G_k^{\ast }({\tau }_1)].$$(3.32)

If we now consider the limit ${\tau }_1\to {\tau }_2$, we have, from Eqs. (3.29), (3.30) and (3.31), (3.32), that

$$〈{ĥ}_{ij}(\overrightarrow{k},\tau ){ĥ}_{mn}(\overrightarrow{p},\tau )〉=\frac{2{\pi }^2}{k^3}{\mathcal{𝒮}}_{ijkmn}(\widehat{k})P_T(k,\tau ){\delta }^{(3)}(\overrightarrow{k}+\overrightarrow{p}),$$(3.33)

$$〈{\partial }_{\tau }{ĥ}_{ij}(\overrightarrow{k},\tau ){\partial }_{\tau }{ĥ}_{mn}(\overrightarrow{p},\tau )〉=\frac{2{\pi }^2}{k^3}{\mathcal{𝒮}}_{ijkmn}(\widehat{k})Q_T(k,\tau ){\delta }^{(3)}(\overrightarrow{k}+\overrightarrow{p}).$$(3.34)

Equations (3.33) and (3.34) reproduce exactly Eqs. (3.11) and (3.12) with the difference that random fields are replaced by the field operators and the ensemble averages are now quantum mechanical expectation values. Furthermore the two power spectra $P_T(k,\tau )$ and $Q_T(k,\tau )$ depend on the specific form of the mode functions Since Eqs. (3.11) and (3.12) hold also in the quantum case, the production of relic gravitons can be certainly mimicked with a homogeneous process. To analyze of the stationarity we need the explicit form of the autocorrelation functions, at late times.⁷⁶ For this purpose, the phases must be correctly determined from the continuity of the mode functions and from the Wronskian normalization condition of Eq. (3.26). These requirements ultimately lead to the standing oscillations^q both in the mode functions and in the autocorrelation functions. This means that the autocorrelation functions at late times do not only depend on the time-difference $|{\tau }_1-{\tau }_2|$ (as it would happen in the case of a stationary process) but also on $({\tau }_1+{\tau }_2)$. The nonstationary features are simpler to illustrate in the radiation epoch since the explicit expressions are less cumbersome; during the radiation stage the mode functions can be expressed in terms of a $2\times 2$ unitary matrix where, as already mentioned after Eq. (3.26), the mode functions have been rescaled as $F_k(\tau )=f_k(\tau )∕a(\tau )$ and $G_k(\tau )=g_k(\tau )∕a(\tau )$; in Eq. (3.36) ${\overline{f}}_k$ and ${\overline{g}}_k$ indicate the initial conditions of the mode functions as they emerge from the inflationary stage of expansion. To describe the smooth evolution of the background and of the mode functions in Eq. (3.36) it is appropriate to introduce the variable $u(\tau )$ where ${\tau }_r$ conventionally marks the onset of the radiation-dominated stage and $ϵ$ is the standard slow-roll parameter; by definition $u_r=u(-{\tau }_r)=k(1-ϵ){\tau }_r$. The inflationary initial conditions determine the amplitude of the tensor power spectrum during inflation for typical wavelengths larger than the Hubble radius and also the normalization of the autocorrelation function. In particular, we can introduce where we defined, for the sake of conciseness, ${\overline{P}}_T^{(r)}={\overline{P}}_T(k,{\tau }_r)$ and $n_T=-2ϵ=-r_T∕8$. In Eq. (3.38) V denotes the inflationary potential which is related to the expansion rate in the slow-roll approximation $3H^2{\overline{M}}_P^2\simeq V$; as already mentioned prior to equation and in the related footnote $H_r$ is, roughly speaking, the expansion rat at the end of inflation. It is relevant to point out that in the limit $n_T\to 0$ the last equality in Eq. (3.38) is ill defined even though it is still true that ${\overline{P}}_T^{(r)}=(16∕\pi ){(H_r∕M_P)}^2$. In summary, during the radiation stage the autocorrelation functions of Eqs. (3.31) and (3.32) become

$${\mathrm{\Delta }}_k({\tau }_1,{\tau }_2)=\frac{{\pi }^2{\overline{P}}_T^{(r)}}{k^3}\frac{[cos(u_1-u_2)-cos(u_1+u_2)]}{u_1{u}_2},$$(3.39)

$${\overline{\mathrm{\Delta }}}_k({\tau }_1,{\tau }_2)=\frac{{\pi }^2{\overline{P}}_T^{(r)}}{k}\left[\frac{cos(u_1-u_2)}{u_1{u}_2}\left(1+\frac{1}{u_1{u}_2}\right)+\frac{sin(u_1-u_2)}{u_1{u}_2}\left(\frac{1}{u_2}-\frac{1}{u_1}\right)\right.$$

$$\left.+\frac{cos(u_1+u_2)}{u_1{u}_2}\left(1-\frac{1}{u_1{u}_2}\right)-\frac{sin(u_1+u_2)}{u_1{u}_2}\left(\frac{1}{u_2}+\frac{1}{u_1}\right)\right],$$(3.40)

where, by definition, $u_1=u({\tau }_1)$ and $u_2=u({\tau }_2)$; at late times $u_1\simeq k{\tau }_1$ and $u_2=k{\tau }_2$. The autocorrelation functions of Eqs. (3.39) and (3.40) do not only depend on the time difference (as implied in the case of stationary processes); on the contrary both ${\mathrm{\Delta }}_k({\tau }_1,{\tau }_2)$ and ${\overline{\mathrm{\Delta }}}_k({\tau }_1,{\tau }_2)$ include distinct functions of $({\tau }_1-{\tau }_2)$ and of $({\tau }_1+{\tau }_2)$. In Eqs. (3.39) and (3.40), we can also see a number of corrections going as inverse powers of $u_1$ and $u_2$; some of these corrections are suppressed when the wavelengths of the gravitons are much smaller than the Hubble radius during the radiation stage. In the matter stage the discussion is technically more involved but the final result is the same.⁷⁶ This means that the backgrounds of relic gravitons are homogeneous in space but they cannot be reduced to a stationary stochastic process. This conclusion has various implications that are not analyzed here (see, however, Ref. 76). It is appropriate to remark, however, that the use of the spectral amplitude should be limited to the signals that are describable in terms of homogeneous and stationary processes.^r

3.2.2. The averaged multiplicity

In a quantum mechanical perspective the amplification of the field amplitudes corresponds to the creation of gravitons either from the vacuum or from some other initial state. Since the production of particles of various spin in cosmological backgrounds is a unitary process^12,13 (see also Refs. 89, 90,91) which is closely analog to the ones arising in the context of the quantum theory of parametric amplification,^{92,93,94,95,96,97,98} the relation between the creation and the annihilation operators in the asymptotic states is given by

The time-dependent (complex) functions ${\alpha }_{p,\lambda }(\tau )$ and ${\beta }_{p,\lambda }(\tau )$ appearing in Eq. (3.42) satisfy $|{\alpha }_{p,\lambda }(\tau ){|}^2-|{\beta }_{p,\lambda }(\tau ){|}^2=1$; because of the unitary evolution $[{â}_{\overrightarrow{p},\lambda },{â}_{\overrightarrow{k},{\lambda }^{\prime }}^{†}]={\delta }^{(3)}(\overrightarrow{k}-\overrightarrow{p}){\delta }_{\lambda {\lambda }^{\prime }}$ and $[{\widehat{b}}_{\overrightarrow{p},\lambda },{\widehat{b}}_{\overrightarrow{k},{\lambda }^{\prime }}^{†}]={\delta }^{(3)}(\overrightarrow{k}-\overrightarrow{p}){\delta }_{\lambda {\lambda }^{\prime }}$. Since the total three-momentum must be conserved, Eqs. (3.41) and (3.42) describe the production of graviton pairs with opposite momenta and the averaged multiplicity is obtained by computing the mean number of gravitons for with momentum $\overrightarrow{k}$ and $-\overrightarrow{k}$, i.e. where $\overline{n}(k,\tau )=|v_k(\tau ){|}^2$ denotes the multiplicity of the pairs of relic gravitons and the further factor of 2 counts the polarizations. From Eq. (3.43), the spectral energy density in critical units is expressed in terms of the averaged multiplicity of the produced gravitons with opposite three-momenta as The result of Eq. (3.44) does not show any specific dependence on $\hslash$ and c just because we are working here in the natural system of units where $\hslash =c={\kappa }_B=1$. The $\hslash$ dependence can be however restored by recalling that the energy of a single graviton is given by $\hslash \omega$ (which we simply wrote as k in units $\hslash =c=1$); moreover another $\hslash$ is present in the definition of Planck mass squared. Thus, as suggested in Ref. 99${\mathrm{\Omega }}_{gw}(\nu ,{\tau }_0)\propto {\hslash }^2$ and this means that the relic gravitons have a truly quantum mechanical origin since their energy density goes to zero in the limit $\hslash \to 0$.

3.2.3. Upper bound on the maximal frequency of the spectrum

Since we are here normalizing the scale factor as $a_0=1$, the physical and the comoving frequencies coincide at the present time and from Eq. (3.44) the spectral energy density in critical units becomes

where ${\tau }_0$ denotes the current value of the conformal time coordinate. Equation (3.45) suggests that the maximal frequency of the spectrum corresponds to the production of a single pair of gravitons with opposite three-momenta. For this reason we can always refer Eq. (3.45) to a putative maximal frequency of the spectrum (be it ${\nu }_{max}$) and obtain^100,101 While in Eqs. (3.45) and (3.46) the single graviton limit¹¹⁰ is achieved when $\overline{n}({\nu }_{max},{\tau }_0)\to 1$, in a classical perspective the maximal frequencies correspond to the bunch of wavenumbers that experience the minimal amplification. All the wavelengths reentering the Hubble radius between the end of inflation and BBN must comply with the bound^{s102,103,104,105,106} where ${\mathrm{\Omega }}_{\gamma 0}$ is the (present) critical fraction of CMB photons. Equation (3.47) sets an indirect constraint on the extra-relativistic species possibly present at the time of nucleosynthesis. Since Eq. (3.47) is also relevant in the context of neutrino physics the limit is often expressed in terms of $\mathrm{\Delta }{N}_{\nu }$ (i.e. the contribution of supplementary neutrino species). The actual bounds on $\mathrm{\Delta }{N}_{\nu }$ range from $\mathrm{\Delta }{N}_{\nu }\le 0.2$ to $\mathrm{\Delta }{N}_{\nu }\le 1$ so that the integrated spectral density in Eq. (3.47) must range, at most, between $10^{-6}$ and $10^{-5}$. The averaged multiplicity for $\nu \ll {\nu }_{max}$ corresponds to the mean number of produced pairs and it is approximately given by where $n_T$ denotes, in practice, to the spectral slope of ${\mathrm{\Omega }}_{gw}(\nu ,{\tau }_0)$ in a given frequency interval. In the conventional lore where the consistency relations are enforced $n_T=n_T^{(\mathrm{\text{low}})}\simeq -r_T∕8+\mathcal{𝒪}(r_T^2)$. There are, however, different physical situations where $n_T>0$^45,46,47 or even $n_T>1$ (see, for instance, Ref. 49). In all these situations ${\mathrm{\Omega }}_{gw}(\nu ,{\tau }_0)$ increase at high frequencies while the averaged multiplicity for $\nu \ll {\nu }_{max}$ is comparatively less suppressed than in the standard lore where $n_T\simeq n_T^{(\mathrm{\text{low}})}\to 0$. The pair production process implies that for $\nu >{\nu }_{max}$ the averaged multiplicity is suppressed^12,13,90 : The degree of suppression of Eq. (3.49) depends on $\gamma$, i.e. a numerical factor $\mathcal{𝒪}(1)$ controlled by the smoothness of the transition between the inflationary and the post-inflationary phase; the value of $\gamma$ can be numerically estimated if the evolution of the mode functions is carefully integrated frequency by frequency.^107,108 The mean number of pairs produced from the vacuum can be written in the following suggestive form : where $n_T$ coincides, in practice, with the high-frequency spectral index $n_T^{(\mathrm{\text{high}})}$. Equation (3.50) interpolates between the power-law behavior of Eq. (3.48) and the exponential suppression of Eq. (3.49) and suggests that the spectrum of relic gravitons should be represented by a distorted thermal spectrum as argued long ago.⁸⁹ There are furthermore situations where $n_T^{(\mathrm{\text{high}})}\to 3$ and this happens in some bouncing scenarios (see e.g. Ref. 49 and discussions therein). It is therefore possible that a thermal spectrum of relic gravitons is produced purely from geometric effects, as originally suggested in Ref. 89. We should also stress, incidentally, that the existence of the exponential suppression for $\nu >{\nu }_{max}$ guarantees the convergence of the integral (3.47) also in the case when the integration is performed up to $\nu \to \infty$. If, for some reason, the diffuse background of relic gravitons has been formed after, BBN ${\mathrm{\Omega }}_{gw}(\nu ,{\tau }_0)$ must always be smaller than the current fraction of relativistic species to avoid an observable impact on the CMB and matter power spectra.^105,109 Inserting therefore Eq. (3.50) into Eq. (3.46) we can directly use Eq. (3.47) (or its analog at even later times) to obtain a general bound on ${\nu }_{max}$ : where the inequality follows in the limit $\overline{n}({\nu }_{max},{\tau }_0)\to \mathcal{𝒪}(1)$, i.e. in the case where a single pair of gravitons is produced.^101,100 The argument leading to the bound of Eq. (3.51) follows directly from the quantum mechanical result of Eq. (3.44) that vanishes in the limit $\hslash \to 0$; in short we can therefore say that, according to a purely quantum perspective, ${\nu }_{max}$ coincides with the frequency where only a graviton pair is produced. Since $h_0^2{\mathrm{\Omega }}_{gw}(\nu ,{\tau }_0)$ scales as ${\nu }^4$ when $\overline{n}(\nu ,{\tau }_0)$ is (approximately) frequency-independent single gravitons (or bunches of coherent gravitons) could be preferentially detected in the high-frequency range: this observation has been eventually pointed out in Ref. 110 and a similar argument is in fact due to Dyson¹¹¹ who suggested that only at high frequencies it might be eventually possible to detect single gravitons. The existence of ${\nu }_{max}$ rules out the possibility of considering gravitons of arbitrary large frequencies as sometimes assumed by those who prefer to ignore the physical implications of high-frequency gravitons. On the contrary, for frequencies $\nu =\mathcal{𝒪}({\nu }_{max})$, as we shall see in the second part of this section, $h_0^2{\mathrm{\Omega }}_{gw}(\nu ,{\tau }_0)$ may exceed the signal of the concordance scenario and may even exceed the contribution of the single-graviton line at lower frequencies.¹¹⁰ As already suggested in the past^45,46,47 high-frequency detectors might resolve single-gravitons.^110,111 In the perspective of Ref. 111 this could happen by conversion to photons in a strong magnetic field¹¹² with experimental techniques very similar to the ones employed for the scrutiny of axion-like particles¹¹³ (see also Refs. 114, 115, 116, 117, 118 for some other papers with similar inspiration). In the case of Eq. (3.51) the relic gravitons cannot exceed the THz domain where coupled microwave cavities with superconducting walls,^{119,120,121,122,123,124} waveguides^{125,126,127,128,129} or even small interferometers^130,131,132 could be used for direct detection even if the current sensitivities should not be overestimated as often done in recent times.¹¹⁰ The observation of Ref. 110 raised a debate on the possibility of assessing the quantumness of the relic gravitons by looking at the analysis of the Bose–Einstein correlations.^100,101 In this second perspective the quantumness of the relic gravitons does not rest on the possibility of literally detecting a single gravitons (as sometimes misunderstood) but rather on the correlation properties of the underlying macroscopic quantum state.^t101

3.3. The expansion history and the spectral energy density

3.3.1. The maximal frequencies

While the bound on ${\nu }_{max}$ deduced in Eq. (3.51) follows from quantum mechanical considerations, in a classical perspective the maximal frequency is computed from the smallest wavelength that crosses the Hubble radius of 4 and immediately reenters; this is why Eq. (3.52) depends upon $H_1\simeq H$ and also upon the timeline of the post-inflationary expansion rate discussed in Sec. 2. Let us therefore start from the simplest situation where the post-inflationary evolution is dominated by radiation. In terms of the cartoons of Figs. 4 and 6 this means that all the ${\delta }_i\to 1$. Since in this case we already denoted the number of e-folds with an overline (e.g. ${\overline{N}}_{max}$, ${\overline{N}}_k$ and do on) we are now going to indicate the maximal frequency deduced in this case by ${\overline{\nu }}_{max}$:

If the transition to radiation dominance is almost sudden, ${\tau }_r$ coincides approximately with ${\tau }_1$; for notational convenience we also use the convention $H\simeq H_1$ where H indicates the inflationary expansion rate estimated from the power spectrum of the curvature inhomogeneities. Equation (3.53) does not assume a specific relation between ${ϵ}_k$ and $r_T$ however, if the consistency relations are enforced, we can always trade ${ϵ}_k$ for $r_T$ and the value of ${\overline{\nu }}_{max}$ becomes The impact of $\mathcal{𝒞}(g_s,g_{\rho },{\tau }_r,{\tau }_{\mathrm{\text{eq}}})$ on ${\overline{\nu }}_{max}$ is minor; for typical values of the late-time parameters (i.e. $g_{\rho ,r}=g_{s,r}=106.75$ and $g_{\rho ,\mathrm{\text{eq}}}=g_{s,\mathrm{\text{eq}}}=3.94$) $\mathcal{𝒞}(g_s,g_{\rho },{\tau }_r,{\tau }_{\mathrm{\text{eq}}})=0.7596$ and the determination of ${\overline{\nu }}_{max}$ of Eq. (3.54) moves from ${\overline{\nu }}_{max}=271.93$ to $206.53$MHz. When the timeline of the post-inflationary evolution is not dominated by radiation but by a generic sequence of stages expanding either faster or slower than radiation (see Figs. 4 and 6 and discussions therein) the maximal frequency can be related to ${\overline{\nu }}_{max}$ and is given by When all the ${\delta }_i\to 1$ the value of ${\nu }_{max}$ coincides with the ${\overline{\nu }}_{max}$ of Eq. (3.52). In case all the ${\delta }_i$ are equal (i.e. ${\delta }_i=\delta$) the post-inflationary evolution consists of a single stage. The product of all the ${\xi }_i$ then coincides with ${\xi }_r=H_r∕H$, as explained in Eqs. (2.55) and (2.56). Provided $\delta <1$ (i.e. when the expansion rate is slower than radiation) ${\nu }_{max}>{\overline{\nu }}_{max}$; conversely, when $\delta >1$ (and the expansion rate is faster than radiation) ${\nu }_{max}<{\overline{\nu }}_{max}$. According to Eq. (3.54), the value of ${\overline{\nu }}_{max}$ is $\mathcal{𝒪}(300)$MHz. This means that if the post-inflationary evolution is dominated by an expanding stage with $\delta \to 1∕2$ with $H_r=\mathcal{𝒪}(10^{-30})M_P$ the value of ${\nu }_{max}$ is going to be $\mathcal{𝒪}(10)$GHz. Similarly if $\delta \to 2$ (and with the same choice of $H_r$) ${\nu }_{max}=\mathcal{𝒪}(100)$kHz. In summary we can say that

•	in a model-dependent perspective the maximal frequency of the relic gravitons obeying the bound (3.51) is sensitive to the timeline of the post-inflationary expansion rate;
•	in the case of a radiation-dominated evolution extending throughout the post-inflationary stage the maximal frequency is of the order of 300MHz;
•	if the post-inflationary expansion rate is smaller than radiation for some time ${\nu }_{max}>\mathcal{𝒪}(300)$MHz;
•	if the expansion rate is instead faster than radiation ${\nu }_{max}<\mathcal{𝒪}(300)$MHz;
•	in general terms, recalling the considerations of Sec. 2, we have that $\mathcal{𝒪}(100)\mathrm{\text{kHz}}<{\nu }_{max}<\mathrm{\text{THz}}$.

Although the maximal frequency alone cannot be used to determine observationally the timeline of the expansion rate, Eqs. (3.53)–(3.55) suggest nonetheless that ${\nu }_{max}$ of the spectrum is sensitive to all the aspects of the post-inflationary evolution.^u

3.3.2. The intermediate frequencies

From Figs. 4 and 6, we have that the bunch of frequencies $\nu =\mathcal{𝒪}({\nu }_{max})$ corresponds to the wavelengths that left the horizon at the end of inflation and reentered immediately after. Depending on the timeline of the post-inflationary evolution there are other typical frequencies that can be explicitly computed. Moreover, since ${\overline{\nu }}_{max}$ depends on $r_T$, also all the other frequencies are sensitive to the specific value of the tensor-to-scalar-ratio. Rather than starting from the general considerations it is better to consider a specific example. Let us then suppose that, before the onset of radiation dominance, the post-inflationary epoch consists of thee separate phases. This means, according to Figs. 4 and 6, that the final spectrum is going to be characterized by the three typical frequencies ${\nu }_1={\nu }_{max}$, ${\nu }_2$ and ${\nu }_3={\nu }_r$. As already stressed after Eq. (2.55) we actually recall that we can always assume that $a_n=a_r$ and ${\nu }_r={\nu }_n$ so that the nth stage of expansion corresponds (by construction) to the radiation phase. In the case $n=3$ the expression of ${\nu }_{max}$ follows from Eq. (3.55) and it is

where ${\xi }_1=H_2∕H_1$ and ${\xi }_2=H_r∕H_2$. The intermediate frequencies ${\nu }_2$ and ${\nu }_r$ are related to ${\overline{\nu }}_{max}$ and they are where, by definition, ${\xi }_r={\xi }_1{\xi }_2=H_r∕H_1$. Equation (3.57) can be easily generalized so that when n intermediate phases are present prior to $a_r$ the generic intermediate frequencies ${\nu }_m$ and ${\nu }_r$ are

$${\nu }_m=\sqrt{{\xi }_1}\dots \sqrt{{\xi }_{m-1}}\prod _{m=1}^{n-1}{\xi }_i^{\frac{{\delta }_i-1}{2({\delta }_i+1})}{\overline{\nu }}_{max},$$(3.58)

$${\nu }_r={\nu }_n=\sqrt{{\xi }_1}\sqrt{{\xi }_2}\dots \sqrt{{\xi }_{n-2}}\sqrt{{\xi }_{n-1}}{\overline{\nu }}_{max}.$$(3.59)

Recalling the remarks presented before, since the different phases must not last below $H_r$, the product of all the ${\xi }_i$ equals $H_r∕H_1$, i.e. by definition ${\xi }_1{\xi }_2\dots {\xi }_{n-1}{\xi }_n={\xi }_r=H_r∕H_1$. Therefore, in case the consistency relations are enforced, Eqs. (3.56), (3.57) and (3.58), (3.59) show that both the maximal and the intermediate frequencies of the spectrum depend on $r_T$ through ${\xi }_r$. Since $m=1,\dots ,n-2$, if there are n different stages there are $(n-2)$ intermediate frequencies between ${\nu }_1$ and ${\nu }_r$. If $n=3$, as exemplified above, the only intermediate frequency is ${\nu }_2$ and it is given in Eq. (3.57).

3.3.3. The slopes of the spectra

In the previous subsection, we derived the typical frequencies of the spectrum in the case of a generic sequence of post-inflationary stages with expansion rates that can be either faster or slower than radiation. Within the same framework we could now discuss the slopes of ${\mathrm{\Omega }}_{gw}(k,\tau )$ within the various frequency domains. The calculation of the spectral energy density can be sometimes carried on in analytic terms but more often with appropriate numerical techniques. Here, we shall not review all these aspects but just remark that, for a sound estimate of the spectral slopes, it is sufficient to employ an approximate description that is based on the Wentzel–Kramers–Brillouin (WKB) solution of the mode functions (see, for instance, Ref. 37 and discussion therein). If the power spectra $P_T(k,\tau )$ and $Q_T(k,\tau )$ of Eq. (3.35) are inserted into Eq. (3.9) ${\mathrm{\Omega }}_{gw}(k,\tau )$ can be directly expressed in terms of the mode functions

We note that $F_k(\tau )$ and $G_k(\tau )$ can also be rescaled, i.e. $a(\tau )F_k(\tau )=f_k(\tau )$ and $a(\tau )G_k(\tau )=g_k(\tau )$; in this way Eq. (3.60) becomes Before a given wavelength exits the Hubble radius (see Fig. 6) the mode functions are simple plane waves normalized as in Eq. (3.26) to preserve the canonical commutation relations between field operators. After the wavelengths reenter the Hubble radius $f_k(\tau )$ and $g_k(\tau )$ are

$$f_k(\tau )=\frac{e^{-ik{\tau }_{\mathrm{\text{ex}}}}}{\sqrt{2k}}[A_k({\tau }_{\mathrm{\text{ex}}},{\tau }_{\mathrm{\text{re}}})cosk\mathrm{\Delta }\tau +B_k({\tau }_{\mathrm{\text{ex}}},{\tau }_{\mathrm{\text{re}}})sink\mathrm{\Delta }\tau ],$$(3.62)

$$g_k(\tau )=e^{-ik{\tau }_{\mathrm{\text{ex}}}}\sqrt{\frac{k}{2}}[-A_k({\tau }_{\mathrm{\text{ex}}},{\tau }_{\mathrm{\text{re}}})sink\mathrm{\Delta }\tau +B_k({\tau }_{\mathrm{\text{ex}}},{\tau }_{\mathrm{\text{re}}})cosk\mathrm{\Delta }\tau ],$$(3.63)

where $\mathrm{\Delta }\tau =(\tau -{\tau }_{\mathrm{\text{re}}})$. In Eqs. (3.62) and (3.63) ${\tau }_{\mathrm{\text{ex}}}$ and ${\tau }_{\mathrm{\text{re}}}$ denote, respectively, the moments where a given scale exits and reenters the Hubble radius; the two coefficients $A_k({\tau }_{\mathrm{\text{ex}}},{\tau }_{\mathrm{\text{re}}})$ and $B_k({\tau }_{\mathrm{\text{ex}}},{\tau }_{\mathrm{\text{re}}})$ are given by

$$A_k({\tau }_{\mathrm{\text{ex}}},{\tau }_{\mathrm{\text{re}}})=\left(\frac{a_{\mathrm{\text{re}}}}{a_{\mathrm{\text{ex}}}}\right){\mathcal{𝒥}}_k^{(t)}({\tau }_{\mathrm{\text{ex}}},{\tau }_{\mathrm{\text{re}}}),$$(3.64)

$$B_k({\tau }_{\mathrm{\text{ex}}},{\tau }_{\mathrm{\text{re}}})=\left(\frac{{\mathcal{\mathscr{H}}}_{\mathrm{\text{re}}}}{k}\right)\left(\frac{a_{\mathrm{\text{re}}}}{a_{\mathrm{\text{ex}}}}\right){\mathcal{𝒥}}_k^{(t)}({\tau }_{\mathrm{\text{ex}}},{\tau }_{\mathrm{\text{re}}})-\left(\frac{a_{\mathrm{\text{ex}}}}{a_{\mathrm{\text{re}}}}\right)\left(\frac{{\mathcal{\mathscr{H}}}_{\mathrm{\text{ex}}}+ik}{k}\right),$$(3.65)

$${\mathcal{𝒥}}_k^{(t)}({\tau }_{\mathrm{\text{ex}}},{\tau }_{\mathrm{\text{re}}})=1-({\mathcal{\mathscr{H}}}_{\mathrm{\text{ex}}}+ik){\int }_{{\tau }_{\mathrm{\text{ex}}}}^{{\tau }_{\mathrm{\text{re}}}}\frac{a_{\mathrm{\text{ex}}}^2}{a^2({\tau }^{\prime })}d{\tau }^{\prime }.$$(3.66)

In Eqs. (3.64)–(3.66), with obvious notations, ${\mathcal{\mathscr{H}}}_{\mathrm{\text{ex}}}=\mathcal{\mathscr{H}}({\tau }_{\mathrm{\text{ex}}})$, and ${\mathcal{\mathscr{H}}}_{\mathrm{\text{re}}}=\mathcal{\mathscr{H}}({\tau }_{\mathrm{\text{re}}})$. It can be immediately checked that Eqs. (3.64)–(3.66) together with Eqs. (3.62) and (3.63) imply the Wronskian normalization condition $f_k(\tau )g_k^{\ast }(\tau )-f_k^{\ast }(\tau )g_k(\tau )=i$. If the background expands between $a_{\mathrm{\text{ex}}}$ and $a_{\mathrm{\text{re}}}$ we have that all the terms containing the ratio $(a_{\mathrm{\text{re}}}∕a_{\mathrm{\text{ex}}})\gg 1$ superficially dominate against those proportional to $(a_{\mathrm{\text{ex}}}∕a_{\mathrm{\text{re}}})\ll 1$. If we use this logic in a simplified manner we would keep the dominant terms and completely neglect the subdominant ones; by following this logic the approximate expression of the mode functions becomes

$$f_k(\tau )\simeq \frac{e^{-ik{\tau }_{\mathrm{\text{ex}}}}}{\sqrt{2k}}{\mathcal{𝒥}}_k^{(t)}({\tau }_{\mathrm{\text{ex}}},{\tau }_{\mathrm{\text{re}}})\left(\frac{a_{\mathrm{\text{re}}}}{a_{\mathrm{\text{ex}}}}\right)\left\{\frac{{\mathcal{\mathscr{H}}}_{\mathrm{\text{re}}}}{k}sink\mathrm{\Delta }\tau +cosk\mathrm{\Delta }\tau \right\},$$(3.67)

$$g_k(\tau )\simeq e^{-ik{\tau }_{\mathrm{\text{ex}}}}\sqrt{\frac{k}{2}}{\mathcal{𝒥}}_k^{(t)}({\tau }_{\mathrm{\text{ex}}},{\tau }_{\mathrm{\text{re}}})\left(\frac{a_{\mathrm{\text{re}}}}{a_{\mathrm{\text{ex}}}}\right)\left\{\frac{{\mathcal{\mathscr{H}}}_{\mathrm{\text{re}}}}{k}cosk\mathrm{\Delta }\tau -sink\mathrm{\Delta }\tau \right\}.$$(3.68)

We stress that Eqs. (3.67) and (3.68) are quantitatively correct but they do not faithfully account for the unitary evolution of the field operators since the all the subleading terms have been neglected. These terms are essential if the unitarity is to be restored order by order in the perturbative expansion (see in this respect that discussion at the end of this section). Inserting finally Eqs. (3.67) and (3.68) in Eq. (3.61) we can deduce the explicit expression of the spectral energy density in critical units : Equations (3.67) and (3.68) are valid for $k\gg aH$ (when all the corresponding wavelengths are shorter than the Hubble radius) but they do not satisfy the Wronskian normalization condition owing to their approximate form. Equation (3.69) holds within the same approximations which are adequate for the estimates we are presenting here and in the following sections.

3.3.4. Spectral energy density, exit and reentry

According to Eq. (3.69) the slopes of ${\mathrm{\Omega }}_{gw}(k,\tau )$ in a given range of wavenumbers chiefly depend on the dynamics of the expansion rate at ${\tau }_{\mathrm{\text{ex}}}$ and ${\tau }_{\mathrm{\text{re}}}$. For illustrative purposes we can consider that all the wavelengths of spectrum exited the Hubble radius during a conventional inflationary stage; this is the viewpoint of Figs. 4 and 6. The exit may also occur as in Fig. 5 but this possibility is going to be separately examined in Sec. 5. Since the exit of all wavelengths of the spectrum occurs during inflation,

the scale factor in this regime can be deduced from Eq. (3.70) and it is approximately given by $a(\tau )={(-\tau ∕{\tau }_1)}^{-\mu }$. We shall assume that the reentry takes place during a decelerated stage of expansion not necessarily coinciding with a radiation epoch. If the reentry takes place in a generic $\delta$ phase for $\tau \ge -{\tau }_1$ we have that where the continuity of the scale factor and of the Hubble rate has been enforced. The condition $k{\tau }_{\mathrm{\text{re}}}=\mathcal{𝒪}(1)$ holds provided $\delta \ne 1$. This can be understood by appreciating that the condition of the crossing of a given wavelength is not simply given by $k\simeq aH$ but, more precisely, by This condition ultimately follows from the observation that the evolution of the mode functions can be decoupled in terms of $f_k(\tau )$ as $f_k^{\prime \prime }+[k^2-a^{\prime \prime }∕a]f_k=0$; it is easy to show that the crossing condition obtained from these considerations is $k^2\simeq |a^{\prime \prime }∕a|$ which can also be rewritten as $k^2\simeq a^2{H}^2(2-ϵ)$. When $ϵ\ne 2$ both turning points are regular and this implies that the two solutions of Eq. (3.72) are given, respectively, by $k{\tau }_{\mathrm{\text{ex}}}=\mathcal{𝒪}(1)$ and by $k{\tau }_{\mathrm{\text{re}}}=\mathcal{𝒪}(1)$. For instance when a given wavelength crosses the Hubble radius during inflation we have that $ϵ\ll 1$ and $k\simeq a_{\mathrm{\text{ex}}}{H}_{\mathrm{\text{ex}}}$ that also means, by definition, $k{\tau }_{\mathrm{\text{ex}}}\simeq 1$. Similarly if the given wavelength reenters in a decelerated stage of expansion different from radiation $k\simeq a_{\mathrm{\text{re}}}{H}_{\mathrm{\text{re}}}$. However, if the reentry occurs in the radiation stage (or close to it) we have that ${ϵ}_{\mathrm{\text{re}}}\to 2$ (i.e. $\delta ={\delta }_{\mathrm{\text{re}}}\to 1$) and the condition (3.72) implies that $k{\tau }_{\mathrm{\text{re}}}\ll 1$. Equation (3.69) can be further simplified in the following form : which is valid for $k\tau >1$, $\tau >{\tau }_{\mathrm{\text{re}}}$. Recalling then Eqs. (3.70) and (3.71) the WKB estimate of the spectral energy density becomes^v where the spectral index ${\overline{n}}_T$ determines the slope of $h_0^2{\mathrm{\Omega }}_{gw}(k,\tau )$. If the consistency relations are enforced the slow-roll parameter can be traded for the tensor-to-scalar ratio $r_T$ so that the slope is ultimately given by^w Equation (3.75) implies that the high-frequency spectral slope is, respectively, increasing or decreasing depending if the expansion rate is either slower or faster than radiation : The analysis leading to Eq. (3.74) determines the conventional low-frequency slope which is applicable for the frequencies $\nu <{\nu }_r\simeq a_r{H}_r$ and which is given by Eq. (3.75) evaluated in the limit $\delta \to 1$ : Equation (3.77) corresponds to the slope of the spectral energy density obtained for the transition between a conventional inflationary stage of expansion and a radiation phase. Thanks to the consistency relations, $r_T\simeq 16ϵ$ so that the result of Eq. (3.77) also implies that ${\overline{n}}_T=-2ϵ$. We can similarly have that the high-frequency slope is given, in practice, by ${\overline{n}}_T(\delta ,r_T)$ in the limit $r_T\to 0$ : The spectral index ${\overline{n}}_T(\delta ,r_T)$ can be expressed also in different ways; for instance if we consider that the post-inflationary stage is governed by a perfect barotropic fluid we can also write the spectral index as a function of w (the barotropic index) and $ϵ$ (the slow-roll parameter) : where the slow-roll parameter denoted by $ϵ$ coincides with ${ϵ}_{\mathrm{\text{ex}}}=ϵ({\tau }_{\mathrm{\text{ex}}})$ and it should be a slowly varying function of the conformal time coordinate; in the limit $ϵ\to 0$ and $w\to 1$ the spectral index in the high-frequency region is blue, i.e. ${\overline{n}}_T\to 1$. Different values of w slightly reduce the slope (e.g. for $w\to 2∕3$ we have ${\overline{n}}_T\to 2∕3$) so that $0<{\overline{n}}_T\le 1$ as long as $1∕3<w\le 1$; as expected, the results expressed by Eqs. (3.74)–(3.79) cannot be applied for $ϵ\to 1$ since, in this limit, the slow-roll approximation breaks down.^x

3.3.5. Approximate forms of the averaged multiplicities and unitarity

In the past there have been various attempts to justify the loss of quantum coherence of the relic gravitons by claiming that when particles are copiously produced the averaged multiplicities are very large (see e.g. Ref. 133 and references therein). The averaged multiplicity $\overline{n}(k,\tau )$ accounting for the pairs of gravitons with opposite three-momenta for each tensor polarization follows then from Eqs. (3.41) and (3.42)

The largeness of the averaged multiplicity is (incorrectly) used to argue that the final state of the evolution of the relic gravitons is classical and the argument is, in short, the following. Since $|{\beta }_k(\tau ){|}^2\gg 1$ we also have that $|{\alpha }_k(\tau ){|}^2\gg 1$ and this means that the field operators approximately commute^y : The heuristic argument of Eq. (3.81) is self-contradictory since it suggests that unitarity is approximately lost every time a large number of pairs is produced. This is markedly false. On the contrary when the approximation scheme is accurate the violations of unitarity are neither explicit nor implicit even if the averaged multiplicity of the produced pairs is very large. To clarify this point we first note that, in the WKB approach of Eqs. (3.64)–(3.66) the values of ${\alpha }_k(\tau )$ and ${\beta }_k(\tau )$ are given by

$${\alpha }_k(\tau )=\frac{1}{2}[A_k({\tau }_{\mathrm{\text{ex}}},{\tau }_{\mathrm{\text{re}}})+i{B}_k({\tau }_{\mathrm{\text{ex}}},{\tau }_{\mathrm{\text{re}}})]e^{-ik(\mathrm{\Delta }\tau +{\tau }_{\mathrm{\text{ex}}})},$$(3.82)

$${\beta }_k^{\ast }(\tau )=-\frac{1}{2}[A_k({\tau }_{\mathrm{\text{ex}}},{\tau }_{\mathrm{\text{re}}})-i{B}_k({\tau }_{\mathrm{\text{ex}}},{\tau }_{\mathrm{\text{re}}})]e^{ik(\mathrm{\Delta }\tau +{\tau }_{\mathrm{\text{ex}}})}.$$(3.83)

Following the logic of the argument (3.81) we could now naively argue that all the terms containing the ratio $(a_{\mathrm{\text{re}}}∕a_{\mathrm{\text{ex}}})\gg 1$ superficially dominate against those proportional to $(a_{\mathrm{\text{ex}}}∕a_{\mathrm{\text{re}}})\ll 1$, always in the approximation that the background expands between $a_{\mathrm{\text{ex}}}$ and $a_{\mathrm{\text{re}}}$. If we would use this logic we would simply keep the dominant terms and discard the subdominant ones. This approach violates unitarity and a the correct strategy is instead to expand systematically in a Laurent series ${\alpha }_k(\tau )$ and ${\beta }_k(\tau )$ with the constraint that $|{\alpha }_k(\tau ){|}^2-|{\beta }_k(\tau ){|}^2=1$. The result of this strategy is where $q_{\mathrm{\text{ex}}}=a_{\mathrm{\text{ex}}}{H}_{\mathrm{\text{ex}}}∕k$ and $q_{\mathrm{\text{re}}}=a_{\mathrm{\text{re}}}{H}_{\mathrm{\text{re}}}∕k$. We can immediately verify that, to the given order in the perturbative expansion (i.e. ${(a_{\mathrm{\text{re}}}∕a_{\mathrm{\text{ex}}})}^5$), Eq. (3.84) implies that $|{\alpha }_k(\tau ){|}^2-|{\beta }_k(\tau ){|}^2=1$ so that unitarity is not lost while keeping the leading terms of the expansion. The result of Eq. (3.84) can be further simplified by neglecting various factors so that, ultimately, ${\alpha }_k({\tau }_{\mathrm{\text{ex}}},{\tau }_{\mathrm{\text{re}}})$ and ${\beta }_k({\tau }_{\mathrm{\text{ex}}},{\tau }_{\mathrm{\text{re}}})$ appearing in Eqs. (3.41) and (3.42) can be evaluated as where, again, $|{\alpha }_k({\tau }_{\mathrm{\text{ex}}},{\tau }_{\mathrm{\text{re}}}){|}^2-|{\alpha }_k({\tau }_{\mathrm{\text{ex}}},{\tau }_{\mathrm{\text{re}}}){|}^2=1$. Equation (3.85) should be still used with some attention since, in the derivation, we artificially neglected some phases just for the benefit of a simple expression. However the approximation of Eq. (3.85), unlike other approaches, is at least consistent with unitarity. If the limit $a_{\mathrm{\text{re}}}\gg a_{\mathrm{\text{ex}}}$ is then taken at the very end of the calculation, the previous results for the spectral energy density are explicitly obtained without any violation of the unitary evolution.

4. The Expansion History and the Low-Frequency Gravitons

4.1. General considerations

4.1.1. Enhancements and suppressions of the inflationary observables

The low-frequency range of the relic gravitons falls in the aHz domain and it corresponds to the CMB wavelengths that left the Hubble radius $N_k$e-folds before the end of inflation. As already mentioned in section 2, these wavelengths are $\mathcal{𝒪}({\lambda }_p)$ where ${\lambda }_p=2\pi ∕k_p$ and $k_p=0.002{\mathrm{\text{Mpc}}}^{-1}$ is the pivot scale at which the spectra of the scalar and tensor modes of the geometry are normalized within the present conventions; note in fact that ${\nu }_p=k_p∕(2\pi )=\mathcal{𝒪}(3)$ aHz. The timeline of the post-inflationary evolution directly affects the values of the tensor to scalar ratio and of the other inflationary observables^z through their dependence upon $N_k$ which can be substantially different from $\mathcal{𝒪}(60)$. For instance a stage expanding faster than radiation has been suggested in the past with the purpose of enhancing the values of $r_T$ (see, for instance, Refs. 134, 135, 136). Indeed, if the expansion rate is faster than radiation $N_k$ gets eventually smaller than the value it would have when the post-inflationary evolution is dominated by radiation (see Eq. (2.57) and discussion therein). But since the inflationary observables and the tensor to scalar ratio are all suppressed by different powers of $N_k$, they might all experience a certain level of enhancement as long as the post-inflationary expansion rate is faster than radiation^aa and this is why this possibility has been employed to account for the BICEP2 excess.¹³⁷ Different mechanisms have been suggested for the same purpose like the violation of the consistency relations caused, for instance, by the quantum initial conditions in the case of a short inflationary stage.^138,139 A post-inflationary stage expanding faster than radiation efficiently enhances the value of $r_T$ especially in the case of single-field scenarios with monomial potentials. We now know that the BICEP2 measurements were seriously affected by foreground contaminations so that, at the moment, the current bounds suggests $r_T\le 0.06$^42,43,44; this also means that the observational evidence would suggest that $r_T(N_k)$ is comparatively more suppressed than in the case $N_k={\overline{N}}_k=\mathcal{𝒪}(60)$. In this respect an even earlier suggestion^45,46,47 (discussed well before the controversial BICEP2 observations¹³⁷) implies that the values of the inflationary observables can be further reduced (rather than enhanced) thanks to a stage expanding more slowly than radiation^107,108 (see also Refs. 64, 65); this is ultimately the punchline of the considerations of Sec. 2 (see in particular Eq. (2.57) and discussions therein). As pointed out in Ref. 70 a reduction of $r_T$ (such as the one suggested by current determinations^42,43,44) implies that the spectral energy density of relic gravitons is enhanced for frequencies larger than the kHz. This conclusion is particularly interesting since two widely separated frequency domains (i.e. the aHz and the MHz regions) may eventually cooperate in the actual determination of the post-inflationary expansion history, as originally pointed out in Refs. 107 and 108.

4.1.2. The number of e-folds and the potential

When the pivot scales cross the comoving Hubble radius the values of the inflationary observables can be directly expressed as a function of $N_k$ for $k=\mathcal{𝒪}(k_p)$. For this purpose the values of the slow-roll parameters (and their dependence on $N_k$) must be evaluated not simply for the conventional post-inflationary evolution dominated by radiation but in the case of different expansion rates. The total number of e-folds elapsed since the crossing of the CMB wavelengths follows from Eq. (2.14) and it is given by

where $t_i(k)$ coincides with the crossing of the CMB wavelengths and $t_f$ marks the end of the inflationary stage. In the case of plateau-like potentials it makes sense to normalize $V(\phi )$ in terms of its scale of variation conventionally denoted hereunder by M : Once Eq. (4.2) is inserted into Eq. (4.1) the number of e-folds elapsed since the crossing of the CMB wavelengths becomes where ${\mathrm{\Phi }}_k$ is now the value of the field when the scale k crosses the comoving Hubble radius while ${\mathrm{\Phi }}_f$ coincides with the end of inflation. Even if different approaches can be envisaged we are here suggesting that the end of inflation effectively occurs when From Eq. (4.4) it also follows that where, by definition, ${\rho }_{\phi }=V+{\dot{\phi }}^2∕2$. If ${ϵ}_f\to 1$ we have that ${\dot{\phi }}_f^2=V_f$ and ${\rho }_{\phi }^{(f)}=3V_f∕2$. Equations (4.3) and (4.4) imply a direct connection between ${\mathrm{\Phi }}_k$ and $N_k$ even if the scaling of the various inflationary observables and of the slow-roll parameters may be different. For instance the slow-roll parameter $ϵ(\tau )$ evaluated at $\tau \simeq 1∕k$ (i.e. $ϵ(\tau )=ϵ(1∕k)={ϵ}_k$) scales differently as a function of $N_k$: we could have ${ϵ}_k\propto 1∕N_k$ (as it happens in the case of monomial potentials^134,135,136) or ${ϵ}_k\propto 1∕N_k^2$ (a typical scaling of plateau-like potentials) or even other more complicated scalings like the ones of hill-top potentials^{140,141,142,143} (see also Refs. 144, 145, 146). As suggested in Eq. (2.57) the value of $N_k$ ultimately depends on the post-inflationary evolution and it cannot be reliably fixed without a specific knowledge of the early expansion rate right after the end of inflation. In particular, if the evolution of the background prior to nucleosynthesis is faster than radiation, some of the ${\delta }_i$ in Eq. (2.57) will be larger than 1 and $N_k$ gets comparatively smaller than ${\overline{N}}_k$. The opposite is true if the expansion rate is slower than radiation: in this case some of the ${\delta }_i$ in Eq. (2.57) will be smaller than 1 so that ${ϵ}_k$ will be more suppressed than in the case of radiation dominance where $N_k\to {\overline{N}}_k$. This perspective is further scrutinized and more concretely illustrated hereunder.

4.1.3. Illustrative examples and physical considerations

In the case of plateau-like potentials $v(\mathrm{\Phi })$ may be written as the ratio of two functions approximately scaling with the same power for $\mathrm{\Phi }\gg 1$; for instance we can have

In Eq. (4.6), $\beta$, p and q are the parameters of the potential and, for technical reasons related with the limit $\mathrm{\Phi }\gg 1$, it is practical to require $4q>p$. Under the conditions of spelled out in Eq. (4.6), an oscillating stage may arise in the limit $\mathrm{\Phi }\ll 1$ where $v(\mathrm{\Phi })={\beta }^p{\mathrm{\Phi }}^{2q}$. With the same strategy different concrete examples can be explicitly constructed : The potentials of Eqs. (4.6) and (4.7) depend upon three parameters (i.e. p, q and $\beta$). The examples of Eqs. (4.6) and (4.7) could be further simplified by choosing, for instance, With similar logic we may also consider The examples of plateau-like potentials can be (obviously) multiplied but instead of focusing on the peculiar features of various potentials it is more interesting to consider the effects of different post-inflationary evolutions on the interplay between the large-scale and small-scale constraints, as originally pointed out in Refs. 45, 46,47; along a more technical perspective the same suggestion has been reinstated also in Refs. 107 and 108. The form of the potential is then useful for illustration but from the physical viewpoint the production of the relic gravitons is determined by the evolution of the space–time curvature: while the scalar inhomogeneities (see Appendix A) may depend on the features of the potential, the production of the relic gravitons is directly sensitive to the expansion rate. It is nonetheless useful to consider which potential might lead to an invisible $r_T$ in the aHz range together with a very large signal in the MHz and GHz regions.⁷⁰ The potentials suggested above go along this perspective even though one can find other classes of potentials that may suppress $r_T$ for the CMB wavelengths without leading to a large signal at higher frequencies (see, for instance, Refs. 142, 143 and 147).

4.2. The tensor to scalar ratio

The amplitudes of the tensor and scalar power spectrum are related via $r_T$ which is, in general terms, both scale-dependent and time-dependent :

where $F_k^{(t)}(\tau )$ and $F_k^{(s)}(\tau )$ denote, respectively, the tensor and the scalar mode functions, i.e.

$$F_k^{(t)\prime \prime }+2\mathcal{\mathscr{H}}{F}_k^{(t)\prime }+k^2{F}_k^{(t)}=0,F_k^{(t)\prime }=G_k^{(t)},\mathcal{\mathscr{H}}=a^{\prime }∕a,$$(4.11)

$$F_k^{(s)\prime \prime }+2\mathcal{ℱ}{F}_k^{(s)\prime }+k^2{F}_k^{(s)}=0,F_k^{(s)\prime }=G_k^{(t)},\mathcal{ℱ}=z^{\prime }∕z,$$(4.12)

where, notational convenience, in the following discussion we posit $z_{\phi }(\tau )=z(\tau )$. The definition of $P_T(k,\tau )$ has been already introduced in Eq. (3.35) whereas the scalar power spectrum is also discussed in Appendix A. To avoid confusions the tensor mode functions (denoted by $F_k(\tau )$ in Eq. (3.35)) are distinguished from their scalar counterpart by a superscript (i.e. $F_k^{(t)}(\tau )$). With these notations, the tensor to scalar ratio (4.10) becomes We are now going to evaluate $r_T(k,\tau )$ before and after reentry. In Appendix A, the tensor to scalar ratio is discussed by using the exact solutions for the evolution of the mode functions during the inflationary stage; by construction the analysis of Appendix A applies for wavelengths larger than the Hubble radius during inflation.

4.2.1. The tensor to scalar ratio before reentry

The initial conditions of the temperature and polarization anisotropies of the CMB are set prior to matter–radiation equality (i.e. $\tau <{\tau }_{\mathrm{\text{eq}}}$) when the relevant wavelengths are larger than the Hubble radius. This means that Eq. (4.13) should be evaluated for ${\tau }_{\mathrm{\text{ex}}}\le \tau <{\tau }_{\mathrm{\text{re}}}$ and $k\ll aH$; as before ${\tau }_{\mathrm{\text{ex}}}(k)$ and ${\tau }_{\mathrm{\text{re}}}(k)$ denote, respectively, the moments at which a given wavelength either exits or reenters the Hubble radius (see Fig. 6 and discussions therein). In this approximation Eqs. (4.11) and (4.12) can be independently solved :

where the shorthand notations $z_{\mathrm{\text{ex}}}=z_{\phi }({\tau }_{\mathrm{\text{ex}}})$ and $a_{\mathrm{\text{ex}}}=a({\tau }_{\mathrm{\text{ex}}})$ have been used; note that $z_{\phi }(\tau )$ already appears right after Eq. (A.1) and its definition will not be repeated. In the cartoon of Fig. 6 the solutions of Eq. (4.14) hold for ${\tau }_{\mathrm{\text{ex}}}<\tau <{\tau }_{\mathrm{\text{re}}}$, i.e. for frequencies that are larger than the Hubble radius; both solutions of Eq. (4.14) have been correctly normalized to their respective quantum mechanical initial data. The functions ${\mathcal{𝒥}}_k^{(s)}({\tau }_{\mathrm{\text{ex}}},\tau )$ and ${\mathcal{𝒥}}_k^{(t)}({\tau }_{\mathrm{\text{ex}}},\tau )$ are defined as where, we remind, $\mathcal{ℱ}=z_{\phi }^{\prime }∕z_{\phi }=z^{\prime }∕z$ and $\mathcal{\mathscr{H}}=a^{\prime }∕a$; the integral ${\mathcal{𝒥}}_k^{(t)}({\tau }_{\mathrm{\text{ex}}},\tau )$ already appears in Eq. (3.66). When Eqs. (4.14) and (4.15) are inserted into Eq. (4.10), the explicit form of $r_T(k,\tau )$ is obtained and it is valid in the regime $k<aH$ and ${\tau }_{\mathrm{\text{ex}}}\le \tau <{\tau }_{\mathrm{\text{re}}}$. In the case of single-field inflationary models and for the timeline of the comoving horizon illustrated in Fig. 6 we deduce since ${\mathcal{𝒥}}_k^{(t)}({\tau }_{\mathrm{\text{ex}}},\tau )\simeq {\mathcal{𝒥}}_k^{(s)}({\tau }_{\mathrm{\text{ex}}},\tau )\to 1$; in Eq. (4.17) the second equality follows from $2Ḣ=-{\ell }_P^2{\dot{\phi }}^2$ (see Eqs. (2.18), (2.19) and discussion thereafter).

4.2.2. The tensor to scalar ratio after reentry

The expression of the scalar and tensor mode functions after reentry can be directly obtained from the previous results of Eqs. (3.62), (3.63) and from the subsequent discussion. In particular, within the same approximation leading to Eqs. (3.67), (3.68), the evolution of the tensor mode functions is approximately given by

where, as in Eqs. (3.62) and (3.63) $\mathrm{\Delta }\tau =(\tau -{\tau }_{\mathrm{\text{re}}})$; ${\mathcal{𝒥}}_k^{(t)}({\tau }_{\mathrm{\text{ex}}},{\tau }_{\mathrm{\text{re}}})$ has been already defined in Eq. (4.15) and it is now evaluated for $\tau \to {\tau }_{\mathrm{\text{re}}}$. Equation (4.18) holds when all the comoving frequencies are larger than the expansion rate (i.e. for $k\gg aH$) and in the same approximation $G_k^{(t)}(\tau )$ becomes Equations (4.18) and (4.19) assume an expanding background (i.e. $a_{\mathrm{\text{re}}}\gg a_{\mathrm{\text{ex}}}$) but they are otherwise general since the rates at ${\tau }_{\mathrm{\text{ex}}}(k)$ and ${\tau }_{\mathrm{\text{re}}}(k)$ have not been specified. Always in the limit of short wavelengths, the mode function for the curvature inhomogeneities becomes Equations (4.18) and (4.20) can be finally inserted into the definition of Eq. (4.13) to obtain the wanted form of $r_T(k,\tau )$ valid for $\tau \ge {\tau }_{\mathrm{\text{re}}}$ in the short-wavelength limit (i.e. for $k\tau >1$) : Recalling the discussion of Eq. (3.72), in the limit ${ϵ}_{\mathrm{\text{re}}}\to 2$ we also have ${\mathcal{\mathscr{H}}}_{\mathrm{\text{re}}}∕k\simeq {\mathcal{ℱ}}_{\mathrm{\text{re}}}∕k\gg 1$ and, in this case, $\mathcal{𝒢}(k\mathrm{\Delta }\tau )\to 1$. Conversely, when ${ϵ}_{\mathrm{\text{re}}}\ne 2$ we have instead that ${\mathcal{\mathscr{H}}}_{\mathrm{\text{re}}}∕k\simeq {\mathcal{ℱ}}_{\mathrm{\text{re}}}∕k=\mathcal{𝒪}(1)$; also in this situation $\mathcal{𝒢}(k\mathrm{\Delta }\tau )=\mathcal{𝒪}(1)$. From Eq. (4.21), we can therefore obtain After inflation $ϵ(\tau )$ is in practice piecewise constant and it is of the order of ${ϵ}_{\mathrm{\text{re}}}$ so that, ultimately, $r_T(k,\tau )\to 16{ϵ}_k$ even for short wavelengths.

4.2.3. Oscillating potentials

If the background expands as simple power law $ϵ(\tau )$ is constant; similarly, if the reentry of the given wavelength takes place when the inflaton potential is still dominant (and oscillating) $ϵ(\tau )$ remains approximately constant. To analyze this situation we can first write $ϵ(\tau )$ in terms of the inflaton potential $V(\phi )$, i.e.

As suggested long ago the coherent oscillations of the inflaton imply the approximate constancy of the corresponding energy density.^66,67,68,69 Indeed we have that the evolution of the inflaton energy density ${\rho }_{\phi }$ can be rephrased as During a stage driven by the inflaton oscillations we have that $3H{\dot{\phi }}^2\ll {\dot{\rho }}_{\phi }$ so that, approximately, the energy density is conserved, i.e. ${\dot{\rho }}_{\phi }\simeq 0$, i.e. $V_{max}=V({\phi }_{max})$. Equation (4.25) can be integrated further since the inflaton potential around its minimum can be parametrized as where $x=\phi ∕{\phi }_{max}$. We can now go back to Eq. (4.23); when the numerator and the denominator are averaged over one period of oscillations (say between $\phi =0$ and $\phi ={\phi }_{max}$) $ϵ(\tau )$ becomes Thus, from Eqs. (4.22)–(4.26) and Eq. (4.27), $ϵ(\tau )∕{ϵ}_{\mathrm{\text{re}}}\to 1$ also when the reentry occurs in a stage driven by the coherent inflaton oscillations. With the same approach the average expansion rate can be computed; in particular we can obtain where $\delta$ denotes the expansion rate in the conformal time coordinate (i.e. $a(\tau )={(\tau ∕{\tau }_1)}^{\delta }$). The evolution of the comoving horizon in Fig. 6 assumes a sequence of different expanding stages characterized by the constancy of the expansion rate. A fully equivalent strategy is to consider the continuous variation of $\delta$ implying where ${\rho }_t(a)$ denotes the total energy density governing the post-inflationary evolution prior to radiation. In the case of inflaton-dominated oscillations ${\rho }_t(a)={\rho }_{\phi }$ and By going back to Fig. 6 we therefore have that when the given wavelength crosses the Hubble radius prior to radiation dominance the value of $\delta$ is scale-dependent ${\delta }_k=\delta ({\tau }_{\mathrm{\text{re}}})=\delta (1∕k)$. This conclusion follows by recalling that, during the post-inflationary stage illustrated in the cartoon of Fig. 6, $\delta (a)\ne 1$ which also implies $ϵ(a)\ne 2$.

4.3. Consistency relations and inflationary observables

In this final subsection, we are going to analyze the dependence of the observables upon the post-inflationary timeline encoded in the value of $N_k$. We first introduce the standard form of the slow-roll parameters

and recall that in terms of the dimensionless variables of Eq. (4.2) $ϵ(\tau )$ and $\overline{\eta }(\tau )$ become During inflation all the slow-roll parameters are much smaller than 1 and the corresponding observables at the crossing time become where ${ϵ}_k=ϵ(1∕k)$ and ${\overline{\eta }}_k=\overline{\eta }(1∕k)$ denote the slow-roll parameters evaluated when the bunch of wavelengths corresponding to the CMB scales exited the comoving horizon approximately $N_k$ e-folds prior to the end of inflation. According to the current limits, the tensor-to-scalar ratio and the scalar spectral index are determined as^bb where ${\overline{r}}_T$ ranges between $\mathcal{𝒪}(0.06)$ and $\mathcal{𝒪}(0.03)$ while $0.96448<{\overline{n}}_s<0.96532$ with a central value corresponding to $0.9649$.^42,43,44 For the monomial potentials ${ϵ}_k$ and ${\overline{\eta }}_k$ are of the same order and these scenarios are practically excluded by current data. Let us then consider, as an example, the potential given in Eq. (4.8). In this case from the expression of the number of e-folds we obtain which also implies that

4.3.1. Scaling of the spectral indices with the number of e-folds

When the consistency relations are enforced the tensor to scalar ratio cannot be equally small for all the classes of inflationary potentials and while the monomials are clearly excluded, the plateau-like and the hill-top potentials may lead to $r_T$ that are comparatively smaller. In the case of Eq. (4.6) the explicit expressions of the slow-roll parameters follow from $ϵ(\mathrm{\Phi })$ and $\overline{\eta }(\mathrm{\Phi })$ are given by

In this case, according to Eq. (4.37), the tensor-to-scalar ratio and the scalar spectral index are given by The number of e-folds is ultimately given, in this case, by where we simply assumed ${\mathrm{\Phi }}_f\to 1$. Since the field value at ${\mathrm{\Phi }}_k$ is defined at the time of the crossing during inflation we can take the limit ${\mathrm{\Phi }}_k\gg 1$ in Eq. (4.39) and eventually determine the connection between ${\mathrm{\Phi }}_k$ and $N_k$ : Thanks to Eq. (4.40), Eqs. (4.37) and (4.38) can be directly expressed in terms of ${\mathrm{\Phi }}_k>1$ Finally using Eq. (4.40) into Eq. (4.41) we have In what follows the scaling of the inflationary observables with the number of e-folds and with the other parameters will be swiftly discussed in the case of the example of Eqs. (4.43)–(4.45).

4.3.2. An illustrative example

While the examples along these lines can be multiplied, for the present purposes, different functional forms of the potential do not radically modify the scaling of $r_T(k)$ and of ${\overline{\eta }}_k$. From Eq. (4.42), $n_s(k)$ and $r_T(k)$ becomes

where $n^{(\mathrm{\text{low}})}(k)=-2{ϵ}_k$ denotes the low-frequency spectral index deduced from the enforcement of the consistency relations. Unlike $n_T^{(\mathrm{\text{low}})}(k)$, the high-frequency spectral index may not be unique and it can also increase, as already discussed in Sec. 3. In Fig. 7, we discuss the case of Eqs. (4.43)–(4.45) for $p=1$. In each plot there are two shaded regions: the first one corresponds to the bounds on the scalar spectral index (i.e. $0.96448<{\overline{n}}_s<0.96532$) and it resembles to a vertical stripe that gets wider as q increases; the second shaded area illustrates the bound on $r_T$. When the two shaded regions overlap the constraints on $r_T$ and on $n_s$ are concurrently satisfied. In Fig. 7, we artificially lowered the bounds on $r_T$ (typically $r_T<0.03$) and also considered $N_k$ as a free parameter. A reduction of $r_T$ always entails large values of q; in this case the high-frequency slope of the spectral energy density is increasing, as already discussed in Sec. 3. In Fig. 8, we always illustrate Eqs. (4.43)–(4.45) but for $p=4$. According to Figs. 7 and 8, the region where the constraints are simultaneously satisfied moves toward large q-values where the inflaton oscillations effectively lead to a phase expanding at a rate that is slower than radiation. In this case the high-frequency bound are therefore essential⁷⁰ (see also Refs. 107 and 108).

In summary we have that the low-frequency region is sensitive to the post-inflationary expansion rate through the number of e-folds which can be either larger or smaller than $\mathcal{𝒪}(60)$. If the timeline of the expansion rate is faster than radiation $N_k$ gets smaller and therefore all the inflationary observables are comparatively less suppressed than in the radiation-dominated case. Thanks to the current measurements^42,43,44 we are however in the opposite situation and $N_k$ must comparatively larger than $\mathcal{𝒪}(60)$. In this case the inflationary observables are more suppressed than in the standard radiation-dominated case. Probably the most economical way of achieving this goal is to consider inflationary scenarios where the post-inflationary expansion rate is slower than radiation. In this case, following the considerations of Sec. 4, a high-frequency background of relic gravitons must be expected between the MHz and the THz.

5. The Expansion History and the Intermediate Frequencies

In the intermediate region of the spectrum (extending, approximately, between few pHz and the Hz) two important scales are related, respectively, to the BBN epoch (i.e. ${\nu }_{\mathrm{\text{bbn}}}$) and to the electroweak time (i.e. ${\nu }_{\mathrm{\text{ew}}}$). While ${\nu }_{\mathrm{\text{bbn}}}$ is three orders of magnitude smaller than the observational region of the pulsar timing arrays (PTA in what follows), ${\nu }_{\mathrm{\text{ew}}}$ is comparable with the window where space-borne interferometers might eventually operate a score year from now. During the last four years the PTA reported a series of evidences of gravitational radiation in the nHz range; it is then interesting to understand if these claimed signals are truly primordial or are just coming from diffuse backgrounds of gravitational radiation formed after matter–radiation equality. In any case the PTA set already an indirect constraint on the expansion history of the Universe. Within a similar perspective, the lack of detection between few $\mu$Hz and the Hz (i.e. $\nu \ge {\nu }_{\mathrm{\text{ew}}}$) sets an essential limit on the post-inflationary expansion rate.

5.1. The theoretical frequencies

5.1.1. Neutrino free-streaming

Given the expansion rate at the BBN time (when the temperature of the plasma was approximately $\mathcal{𝒪}(1)$MeV), the general expression of ${\nu }_{\mathrm{\text{bbn}}}$ is

where $g_{\rho ,\mathrm{\text{bbn}}}$ denotes the effective number of relativistic species at the nucleosynthesis epoch. Since $H_0=1.742\times h_{0}10^{-61}{M}_P$, Eq. (5.1) becomes Between the nHz domain and the audio band (with $\mathrm{\text{Hz}}<{\nu }_{\mathrm{\text{audio}}}<10\mathrm{\text{kHz}}$) the spectral energy density of inflationary origin is, at most, $\mathcal{𝒪}(10^{-17})$ in critical units and the deviations from scale-invariance in the direction of blue spectral indices are excluded at least in the conventional situation where the corrections to $h_0^2{\mathrm{\Omega }}_{gw}(\nu ,{\tau }_0)$ always lead to decreasing spectral slope.^cc In the case of the concordance paradigm the spectral energy density is further reduced by various sources of damping and, most notably, by the free-streaming of neutrinos^{160,161,162,163,164} exactly for frequencies below the nHz. The same phenomenon also affects the spectral energy density when the corresponding slopes are increasing^107,108; in both situations, however, the suppression due to the neutrinos operates for $\nu <{\nu }_{\mathrm{\text{bbn}}}$ and when the expansion rate is dominated by radiation.

5.1.2. Big bang nucleosynthesis bound

The frequency range associated with ${\nu }_{\mathrm{\text{bbn}}}$ is related to a set of direct limits on the expansion rate of the plasma at the BBN epoch when the expansion rate was $H_{\mathrm{\text{bbn}}}=\mathcal{𝒪}(10^{-44})M_P$. Any excess in the energy density of the massless species at the BBN time increases the value of $H_{\mathrm{\text{bbn}}}$. The additional massless species may be either bosonic or fermionic; however they are theoretical traditionally parametrized in terms of the effective number of neutrino species as $N_{\nu }=3+\mathrm{\Delta }{N}_{\nu }$. The standard BBN results are in agreement with the observed abundances for $\mathrm{\Delta }{N}_{\nu }\le 1$.^{165,166,167,168} The most constraining bound for the intermediate and high-frequency branches of the relic graviton spectrum is represented by BBN as argued long ago by Schwartzman.¹⁰² The increase in the expansion rate affects, in particular, the synthesis of ${}^4\mathrm{\text{He}}$ and to avoid its overproduction the expansion and rate the number of relativistic species must be bounded from above. All in all, if the additional species are relic gravitons^{102,103,104,105,106} the integral of the spectral energy density over the whole spectrum must satisfy the following bound :

where ${\nu }_{\mathrm{\text{bbn}}}$ is given by Eq. (5.2) and ${\nu }_{max}$ corresponds instead to the maximal frequency of the spectrum. For the relic gravitons produced within the concordance scenario ${\nu }_{max}$ is given by Eq. (5.2). As discussed in Sec. 3 the maximal frequency is model dependent but it is possible to deduce an absolute bound on ${\nu }_{max}$ and, according to this bound, ${\nu }_{max}<\mathrm{\text{THz}}$. Depending on the combined data sets (i.e. various light elements abundances and different combinations of CMB observations), the standard BBN scenario implies that the bounds on $\mathrm{\Delta }{N}_{\nu }$ range from $\mathrm{\Delta }{N}_{\nu }\le 0.2$ to $\mathrm{\Delta }{N}_{\nu }\le 1$. All the relativistic species present inside the Hubble radius at the BBN contribute to the potential increase in the expansion rate and this explains why the integral in Eq. (5.3) must be performed from ${\nu }_{\mathrm{\text{bbn}}}$ to ${\nu }_{max}$. The constraint of Eq. (5.3) can be relaxed in some nonstandard nucleosynthesis scenarios, but, in what follows, the validity of Eq. (5.3) will be enforced by adopting $\mathrm{\Delta }{N}_{\nu }\simeq 1$. The considerations discussed so far can be complemented by other bounds which are, however, less stringent. In particular the same logic employed for the derivation of Eq. (5.3) can be applied at the decoupling of matter and radiation^109,105 when the typical lower extremum of integration becomes ${\nu }_{\mathrm{\text{dec}}}=\mathcal{𝒪}(100)$aHz : The BBN limits examined so far can be relaxed in nonstandard BBN scenarios¹⁰³ (see also Ref. 49). In particular this may happen in the presence of matter–anti-matter domains; instead of being $\mathcal{𝒪}(10^{-5})$ the integral of Eq. (5.3) may get $\mathcal{𝒪}(10^{-4})$.

5.1.3. The electroweak frequency

The Standard Model of particle interactions (based on the $S{U}_L(2)\otimes U_Y(1)\otimes S{U}_c(3)$ gauge group) appears to be successful at least up to energy scales $\mathcal{𝒪}(\mathrm{\text{TeV}})$ and its basic correctness ultimately suggests that the electroweak phase transition cannot produce a detectable background of gravitational radiation for typical frequencies smaller than the Hz. To explain this viewpoint we start by remarking that the dynamics of the electroweak phase transition has been studied since the early 1970s and while it is plausible that spontaneously broken symmetries are restored at high-temperatures, the order of the electroweak phase transition determines the physical features of the purported gravitational signal. The symmetry breaking phase transitions may cause departures from local thermal equilibrium (and from homogeneity) but, according to the current experimental evidence, the electroweak phase transition does not lead to large anisotropic stresses that could eventually produce a diffuse background of gravitational radiation. A large anisotropic stress can only be produced if the electroweak phase transition is of first order and proceeds through the formation of bubbles of the new phase. It was clear already from the first (perturbative) estimates that the electroweak phase transition cannot be strongly first order^169,170,171; however a definite conclusion on this issue was delayed because of the hope that, by using nonperturbative techniques,¹⁷² the essence of the perturbative result could be somehow disproved. The phase diagram of the electroweak theory at high-temperature has been first analyzed by reducing the theory from four to three dimensions and by subsequently simulating on the lattice the lower-dimensional theory with compactified time coordinate.^173,174,175 These analyses have been later corroborated by genuine four-dimensional lattice simulations discussing the $SU(2)$-Higgs system.^176,177 The main results relevant for the present discussion can be summarized, in short, as follows. For approximate values of the Higgs mass $m_H$ smaller than the W-boson mass the phase diagram of the electroweak theory contains a line of first-order phase transitions but for $m_H\ge \mathcal{𝒪}(75)$ GeV the phase transition if of higher order and when $m_H\gg m_W$ (as it is the case from an experimental viewpoint) the phase transition disappears since we can pass from the symmetric to the broken phase in continuous manner. In this cross-over regime there large deviations from homogeneity do not arise and diffuse backgrounds of gravitational radiation are absent.

Although the electroweak phase transition is of higher order, strongly first-order phase transitions may anyway lead to bursts of gravitational radiation and, for this reason, the production of gravitational waves has been investigated in a number of hypothetical first-order phase transitions. Provided the phase transition proceeds thanks to the collision of bubbles of the new phase, the lower frequency scale of the burst is (at most) comparable with the Hubble radius at the corresponding epoch. Denoting by ${\nu }_b$ the frequency of the purported burst, we should always require that ${\nu }_b\ge \mathcal{𝒪}({\nu }_{\mathrm{\text{ew}}})$ where ${\nu }_{\mathrm{\text{ew}}}$ is the typical frequency corresponding to the electroweak horizon. This condition follows directly from the observation that gravitational waves should be formed inside the Hubble radius when the expansion rate of the Universe was approximately $\mathcal{𝒪}(H_{\mathrm{\text{ew}}})$. Assuming the electroweak plasma is dominated by radiation between $H_{\mathrm{\text{ew}}}$ and $H_{\mathrm{\text{bbn}}}$ the electroweak frequency is given by

where the second equality follows from the approximate expansion history and, as usual, $H_0$ denotes the current value of the Hubble rate. Since during the radiation stage $a^{2}H$ is just constant, Eq. (5.5) implies The result of Eq. (5.6) can be further simplified by recalling that If Eq. (5.7) is now inserted into Eq. (5.6) we get the following estimate : where the value of $T_{\mathrm{\text{ew}}}$ has been chosen to be slightly above the value of the top quark mass just to make sure that all the species of the Standard Model are in local thermal equilibrium. Strictly speaking the adiabatic evolution only implies the constancy of $a^3{T}^3{g}_s(T)$ so that result of Eq. (5.8) should be slightly corrected : where $\mathcal{𝒞}(g_s,g_{\rho },{\tau }_r,{\tau }_{\mathrm{\text{eq}}})$ has been already introduced in Eq. (2.43). Equation (5.9) also implies that where the explicit estimate follows by recalling that, for $T_{\mathrm{\text{ew}}}>m_t$, $g_{s,ew}=g_{\rho ,ew}=106.75$. Moreover, since for $T=T_{\mathrm{\text{eq}}}$ the values of $g_{s,\mathrm{\text{eq}}}$ and $g_{\rho ,\mathrm{\text{eq}}}$ are slightly different (i.e. $3.94$ and $3.36$, respectively) the typical value of ${\nu }_{\mathrm{\text{ew}}}$ given in Eq. (5.8) passes from $7.98$ to $6.06\mu \mathrm{\text{Hz}}$.

5.2. Pulsar timing arrays and the expansion history

In the last few years a set of direct observations potentially related with the diffuse backgrounds of gravitational radiation have been reported for a typical benchmark frequency $\mathcal{𝒪}(30)$nHz. This range of frequencies is between 3 and 4 orders of magnitude larger than ${\nu }_{\mathrm{\text{bbn}}}$ and it is currently probed by the Pulsar Timing Arrays (PTAs in what follows). As recently pointed out¹⁷⁸ the signals possibly observed by the PTA may be the result of the pristine variation of the space–time curvature. The specific features of the current observations seem to suggest, however, that $h_0^2{\mathrm{\Omega }}_{gw}(\nu ,{\tau }_0)$ in the nHz domain may only depend on the evolution of the comoving horizon at late, intermediate and early times. This is also, in a nutshell, the systematic perspective swiftly outlined hereunder.

5.2.1. Basic terminology and current evidences

A PTA is just a series of millisecond pulsars that are monitored with a specific cadence that ultimately depends on the choices of the given experiment. We refer here, in particular, (i) to the NANOgrav collaboration,^38,39 (ii) to the Parkes PTA (PPTA)^40,41 and (iii) to the European PTA (EPTA).^179,180 The PTA data have been also combined in the consortium named International PTA (IPTA).¹⁸¹ The data of the PTA collaborations have been released^39,41,180 together with the first determinations of the Chinese PTA (CPTA).¹⁸² As suggested long ago 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.^183,184,185 The signal coming from diffuse backgrounds of gravitational radiation, unlike other noises, should be correlated across the baselines. The 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.¹⁸⁵ If the gravitational waves are not characterized by stochastically distributed Fourier amplitudes the corresponding signal does not necessarily follow the Hellings–Downs correlation. In the past various upper limits on the spectral energy density of the relic gravitons in the nHz range have been obtained^{186,187,188,189} and during the last five years the PTA reported an evidence that could be attributed to isotropic backgrounds of gravitational radiation. The observational collaborations customarily assign the chirp amplitude at a reference frequency ${\nu }_{\mathrm{\text{ref}}}=31.68\mathrm{\text{nHz}}$ that corresponds to ${\mathrm{\text{yr}}}^{-1}$ :

To avoid confusions we stress that the $\beta$ appearing in Eq. (5.11) has nothing to do with the quantity characterizing the inflaton potential (see Eq. (4.6) and discussion thereafter). Recalling now the relation between the spectral energy density and the chirp amplitude, we have where the second equality follows from Eq. (5.11). If we now multiply Eq. (5.12) by $h_0^2$ (where $h_0$ denotes the indetermination in the present value of the Hubble rate) and take into account the explicit value of ${\nu }_{\mathrm{\text{ref}}}$, the expression for the spectral energy density becomes³⁷ In Eq. (5.13), Q is parametrized as $Q=q_0\times 10^{-15}$ (where $q_0$ is a number of order 1) since this is basically the observational evidence. Clearly, for $\nu \to {\nu }_{\mathrm{\text{ref}}}$ implying $h_0^2{\mathrm{\Omega }}_{gw}({\nu }_{\mathrm{\text{ref}}},{\tau }_0)=\mathcal{𝒪}(2.57)\times 10^{-8}$ in the case of Ref. 39 (for $q_0=6.4$) and $h_0^2{\mathrm{\Omega }}_{gw}({\nu }_{\mathrm{\text{ref}}},{\tau }_0)=\mathcal{𝒪}(6.04)\times 10^{-9}$ for Ref. 41 (for $q_0=3.1$). With the same logic we can also deduce the explicit relation between the spectral and the chirp amplitudes : It is also customary to employ $\sqrt{S_h(\nu ,{\tau }_0)}=5.61\times 10^{-12}{q}_0{(\nu ∕{\nu }_{\mathrm{\text{ref}}})}^{\beta -1∕2}{\mathrm{\text{Hz}}}^{-1∕2}$ for a direct comparison it with the spectral amplitude of the signal. Recalling however the considerations of Sec. 3 the use of a spectral amplitude implicitly assumes that the signal can be mimicked by a stationary and homogeneous stochastic process; this is not the case for the relic gravitons.⁷⁶ Bearing in mind the results and the notations of Eqs. (5.12)–(5.15) the main statements of the observers can be summarized, in short, as follows:

•	The pivotal class of models analyzed in Refs. 38, 39, 40, 41, 179–182 always assume $\beta =-2∕3$ (i.e. $\overline{\gamma }=13∕3$); recall, in this respect, that the relation between $\overline{\gamma }$ and $\beta$ is simply given by $\beta =(3-\overline{\gamma })∕2$.
•	In the former data releases the $q_0$ ranged between $1.92$ and $5.13$ depending on the values of $\beta$.38,40,179,181
•	The latest data releases of the Parkes and the NANOgrav collaborations41,39,180 seem to suggest different origins of the diffuse background of gravitational radiation.
•	In particular, after considering 30 ms pulsars spanning 18 years of observations, the PPTA collaboration estimates $q_0=3.1_{-0.9}^{1.3}$ with a spectral index $\beta =-0.45\pm 0.20$41; for a spectral the pivotal model $\beta =-2∕3$ the collaboration suggests instead $q_0=2.04_{-0.22}^{0.25}$ which is compatible with the determinations of the previous data releases40; the PPTA collaboration does not clearly claim the detection of the Hellings–Downs correlation41 and carefully considers possible issues related to time-dependence of the common noise.
•	The conclusions of the PPTA seem significantly more conservative than the one of the NANOgrav collaboration examining 67ms pulsars in the last 15 years.
•	The NANOgrav experiment claims the detection of the Hellings–Downs correlation39 but the inferred values of the spectral parameters are slightly different from the ones of PPTA since $q_0=6.4_{-2.7}^{+4.2}$ and $\beta =-0.10\pm 0.30$.39

5.2.2. The comoving horizon after inflation

The measurements of the PTA set a number of relevant constraints on the spectrum of the relic gravitons and on the expansion rate of the Universe. If the observed excess in the nHz range is just a consequence of the primeval variation of the space–time curvature the spectral energy density of the relic gravitons in the nHz domain only depends on the evolution of the comoving horizon at late, intermediate and early times.¹⁷⁸ Two complementary aspects of the problem will now be addressed. In the first part of the discussion, we are going to see if a post-inflationary modification of the expansion rate can account for the nHz excess. In the second part of the analysis, we consider instead the possibility of explaining the observed PTA excess through the evolution of the effective horizon at early times.

A first general observation is that, in the concordance paradigm, the PTA results do not set any further constraint besides the ones of the aHz region already discussed in Sec. 4. This happens because the spectral energy density of Eqs. (5.15) and (5.14) always exceeds the the one of the concordance paradigm in the nHz region. Indeed, if the expansion rate is dominated by radiation after inflation, $h_0^2{\mathrm{\Omega }}_{gw}(\nu ,{\tau }_0)<\mathcal{𝒪}(10^{-17})$ for typical frequencies larger than ${\nu }_{\mathrm{\text{bbn}}}$. Furthermore, in the concordance paradigm, $h_0^2{\mathrm{\Omega }}_{gw}(\nu ,{\tau }_0)$ is a monotonically decreasing function of the comoving frequency between the aHz and the MHz domain. This means that in the nHz range the signal of the relic gravitons produced within the conventional lore is always 10 orders of magnitude smaller than the one suggested by Eqs. (5.15) and (5.14). If the expansion history is modified in comparison with the concordance paradigm the relevant time-scale of the problem must coincide with ${\tau }_k$, i.e. the moment at which the wavelength associated with ${\nu }_{\mathrm{\text{PTA}}}\simeq {\nu }_{\mathrm{\text{ref}}}=\mathcal{𝒪}(30)\mathrm{\text{nHz}}$ crossed the comoving Hubble radius after the end of inflation (see Fig. 6). The actual value of ${\tau }_k$ represents in fact a fraction of the time-scale associated with BBN :

In Eq. (5.16), we are actually assuming, for simplicity, that ${\nu }_{\mathrm{\text{PTA}}}=\mathcal{𝒪}({\nu }_{\mathrm{\text{ref}}})$ and if ${\nu }_{\mathrm{\text{PTA}}}>{\nu }_{\mathrm{\text{ref}}}$ the corresponding wavelength crossed the comoving horizon even earlier. Besides Eq. (5.16), the second relevant scale of the problem follows from the ratio between ${\nu }_{\mathrm{\text{PTA}}}$ and the expansion rate at the end of inflation : where ${\mathcal{𝒜}}_{\mathcal{ℛ}}$ denotes, as usual, the amplitude of the curvature inhomogeneities at the pivot scale $k_p$ (see Eq. (2.24) and discussion thereafter). Already from Eqs. (5.16) and (5.17) it follows that any modification of the post-inflationary evolution is unlikely to produce a hump for frequencies $\mathcal{𝒪}({\nu }_{\mathrm{\text{PTA}}})$: the value of ${\nu }_{\mathrm{\text{PTA}}}$ in units of the expansion rate is too small. It is on the contrary more likely that a hump will be produced over larger frequencies $\nu >\mu \mathrm{\text{Hz}}$, as we are going to see later on in this section.

To substantiate the previous statement we now consider a generic post-inflationary expanding stage (i.e. a single $\delta$-phase in the language of Sec. 2). When the wavelengths $\lambda =\mathcal{𝒪}({\lambda }_{\mathrm{\text{PTA}}})$ cross the comoving Hubble radius during the $\delta$-phase we have

where, once more, $H_r$ and $H_1$ denote, respectively, the Hubble rates at the onset of the radiation stage and at the end of inflation. As $H_r<H_1$ the comoving horizon at its minimum is comparatively larger for $\delta >1$ than for $\delta \to 1$; for the same reason the opposite is true when $\delta <1$. Since the post-inflationary evolution is modified the spectral energy density of the relic gravitons gets larger, as it follows from Eq. (3.75). When the PTA wavelength crosses the Hubble radius during the $\delta$-phase $h_0^2{\mathrm{\Omega }}_{gw}(\nu ,{\tau }_0)$ exhibits a twofold slope:

•	in the low-frequency regime the slope is simply given by $n_T^{(\mathrm{\text{low}})}=-r_T∕8$; this is true when the consistency relations are enforced as we are assuming throughout;
•	if the wavelength corresponding to ${\nu }_{\mathrm{\text{PTA}}}$ reenters the Hubble radius when $\delta \ne 1$ the high-frequency slope follows from Eqs. (3.75)–(3.77) and it is $n_T^{(\mathrm{\text{high}})}=2(1-\delta )+\mathcal{𝒪}(r_T)$.

To compare $n_T^{(\mathrm{\text{high}})}$ with the potential excesses suggested by the PTA we may recall Eq. (5.13) and then consider the theoretical estimate of the spectral energy density in critical units⁶⁵

where $\mathcal{𝒩}(r_T,\nu )$ includes the effects of the low-frequency suppressions associated with the transfer function, with the neutrino free-streaming^161,162,163 and with the other late-time sources of damping (like the one related with the dark-energy dominance³⁷). In connection with $\mathcal{𝒩}(r_T,\nu )$ see also Eqs. (6.3) and (6.4) and the discussion therein. For $r_T=0.03$ we can numerically estimate that $\mathcal{𝒩}(r_T,\nu )=10^{-16.8}$ and since above ${\nu }_{\mathrm{\text{bbn}}}$ the value of $\mathcal{𝒩}(r_T,\nu )$ has a mild frequency dependence controlled by the value of the low-frequency slope (i.e. $n_T^{(\mathrm{\text{low}})}\simeq -r_T∕8$) and by the low-frequency transfer function, for the present ends we can assume $\mathcal{𝒩}(r_T,\nu )\simeq \mathcal{𝒩}=\mathcal{𝒪}(10^{-17})$ which is the value already quoted before. According to Eq. (5.19), the theoretical amplitude at ${\nu }_{\mathrm{\text{PTA}}}\simeq {\nu }_{\mathrm{\text{ref}}}$ ultimately depends upon ${\nu }_r={\overline{\nu }}_{max}\sqrt{\xi }$ where ${\overline{\nu }}_{max}=\mathcal{𝒪}(300)\mathrm{\text{MHz}}$ has been already computed in Eq. (3.54). In spite of the specific values of ${\overline{\nu }}_{max}$ we have that ${\nu }_r$ cannot be smaller than ${\nu }_{\mathrm{\text{bbn}}}$. In the general case^dd (i.e. when the special value $\beta =-2∕3$ is not preliminarily selected) the PPTA collaboration⁴¹ suggests that $\beta =-0.45\pm 0.20$. This determination is marginally compatible with the value of Eq. (5.20) in the limit $\delta \ge 1∕2$ and the discrepancy between the observational determination of $\beta$ and the values predicted by Eq. (5.20) becomes even more significant if we look at that NANOgrav data suggesting³⁹ $\beta =-0.10\pm 0.30$. All in all, when $\delta \ge 1∕2$, the relation between $\delta$, $\beta$ and $r_T$ is so that, from Eq. (5.20), $\delta =-\beta +\mathcal{𝒪}(r_T)$. Both in the previous^38,40 and in the most recent^39,41 data releases the value of $2(1+\beta )$ is always positive definite (i.e. $1+\beta >0$). In the special case $\beta \to -2∕3$, Eq. (5.20) implies $\delta =2∕3+\mathcal{𝒪}(r_T)$. If we would now assume that the post-inflationary evolution is driven by a relativistic and irrotational fluid we would have $\delta =2∕(3w+1)$ implying that $\beta \to -2∕3$ for $w\to 2∕3$. Another possibility would be that the effective expansion rate is dictated by an oscillating scalar field (like the inflaton) with potential $V(\phi )=V_0{(\phi ∕{\overline{M}}_P)}^{2q}$; in this case the expansion rate during the oscillating phase would be given by $\delta =(q+1)∕(2q-1)$ suggesting that $q=\mathcal{𝒪}(5)$ for $\beta =\mathcal{𝒪}(-2∕3)$. Although, for specific values of $\delta$, the theoretical and the observed slopes can be compatible the corresponding amplitudes involve orders of magnitude that are grossly different; to analyze this aspect we then impose that Eqs. (5.13) and (5.19) should coincide at ${\nu }_{\mathrm{\text{ref}}}$ If the two sides of Eq. (5.21) would be mutually consistent the post-inflationary modification of the comoving horizon might indeed explain the observed PTA excess. But unfortunately the left-hand side of Eq. (5.21) is systematically smaller than the right-hand side; the two contributions are of the same order only when ${\nu }_r\ll {\nu }_{\mathrm{\text{ref}}}$ while, at the same time, $n_T^{(\mathrm{\text{high}})}=2+2\beta$ is sufficiently large and positive. A large (and positive) value of $n^{(\mathrm{\text{high}})}$ guarantees a sharp increase of the spectral energy density while the condition ${\nu }_r\ll {\nu }_{\mathrm{\text{ref}}}$ widens the frequency range for a potential growth of $h_0^2{\mathrm{\Omega }}_{gw}(\nu ,{\tau }_0)$. Since the minimal value of ${\nu }_r$ is provided by ${\nu }_{\mathrm{\text{bbn}}}$ we can select the most favorable situation and posit ${\nu }_r=\mathcal{𝒪}({\nu }_{\mathrm{\text{bbn}}})$. For different values of $\mathcal{𝒩}(r_T,{\nu }_{\mathrm{\text{ref}}})$ (see Eq. (5.19) and discussion thereafter) Eq. (5.21) leads therefore to a specific relation between $\beta$ and $log{q}_0$ :

The result of Eq. (5.22) must then be compared in the plane $(log{q}_0,\beta )$ with the ranges of $\beta$ and $q_0$ determined by the PTA collaborations. The two filled rectangles in Fig. 9 correspond to the observational ranges of $q_0$ and $\beta$; in the same plot the relation between $\beta$ and $log{q}_0$ has been illustrated as it follows from Eq. (5.22) for two neighboring values of $\mathcal{𝒩}(r_T)$. The three diagonal lines of Fig. 9 imply that the values of $\beta$ required to obtain $h_0^2{\mathrm{\Omega }}_{gw}({\nu }_{\mathrm{\text{ref}}},{\tau }_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. 9 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. 9. For the specific case $\beta =-2∕3$ Eq. (5.21) becomes

Again to maximize the potential growth of the spectral energy density we set ${\nu }_r=\mathcal{𝒪}({\nu }_{\mathrm{\text{bbn}}})$ and obtain that the left-hand side of Eq. (5.23) is $2.015\times 10^{-15}$ whereas the right-hand side is always larger than $\mathcal{𝒪}(10^{-9})$. As in the previous case, larger values of ${\nu }_r$ only reduce the left-hand side of Eq. (5.23) and ultimately increase the mismatch between Eqs. (5.13) and (5.19). The argument based on a single $\delta$-phase can be generalized to include different stages expanding either faster or slower than radiation. In this case there is the possibility of developing a hump in the spectrum which is however always much smaller than the excess observed by the PTA.¹⁷⁸

5.2.3. The comoving horizon during inflation

The previous analysis demonstrated that the PTA excess cannot be explained by a post-inflationary modification of the expansion rate. However, if the effective comoving horizon is modified during inflation (as suggested in 5) it is possible to explain the PTA excess in terms of a relic signal.¹⁷⁸ To implement an effective modification of the comoving horizon without obliterating the inflationary expansion we actually consider an effect suggested almost 10 years ago: a dynamical refractive index associated with the propagation of the tensor modes of the geometry in curved backgrounds naturally leads to an increasing spectral energy density at intermediate frequencies.¹⁹⁰ The tensor modes of the geometry may indeed acquire an effective index of refraction when they travel in curved space–times^191,192 and the blue spectral slopes (compatible with the PTA excesses) arise from the variation of the refractive index even if the background geometry evolves according to a conventional stage of expansion possibly supplemented by a standard decelerated epoch¹⁹⁰ (see also Refs. 193, 194, 195). When the refractive index of the relic gravitons is dynamical ($n(a)$ in what follows) the conditions associated with the crossing of a given wavelength are different; the action of the tensor modes of the geometry in the case of a dynamical refractive index^194,195 is given by

see also Eq. (B.26) and discussion therein. The analysis of Eq. (5.24) simplifies if the conformal time coordinate is redefined from $\tau$ to $\eta$ where the relation between the new and the old time parametrizations follows from $n(\eta )d\eta =d\tau$. Equation (5.24) becomes then canonical in terms of a redefined scale factor conventionally denoted hereunder by $b(\eta )$ : The result of Eq. (5.25) explains how and why the evolution of the tensor modes is modified even during a conventional stage of inflationary expansion. The evolution of the tensor amplitude can be directly deduced from Eq. (5.25) and it is where the overdot now denotes a derivation with respect to the $\eta$-time. Equation (5.26) also implies that the standard crossing condition $k^2=a^{\prime \prime }∕a$ is now replaced by $k^2=\ddot{b}∕b$. Between these two conditions the former seems superficially equivalent to the latter but this is not the case.¹⁹⁰ The evolution dictated by Eq. (5.26) ultimately leads to a spectral energy density that increases over intermediate frequencies provided the effective phase velocity^ee of the relic gravitons remains sub-luminal. Although the phase velocity of the relic gravitons is not required to be sub-luminal we impose, for consistency, that $n(a)\ge 1$; in particular we 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 : where $a_i$ and $a_1$ denote, respectively, the beginning and the end of the inflationary epoch; $a_{\ast }$ indicates the boundary of the refractive stage. Equation (5.27) also implies the refractive index is not dynamical in the post-inflationary stage. Some other possibilities have been considered in Refs. 190, 193 and 194 but they will not be examined here. In what follows the parametrization of Eq. (5.27) is regarded as the minimal example that successfully produces a nHz excess. Equation (5.27) can be analyzed in three relevant physical limits: (i) for $a\gg a_1$ we have that $n(a)\to 1$ and the sharpness of the transition depends on the parameter $\overline{s}\ge 1$; (ii) in the range $a_{\ast }<a<a_1$$n(a)$ is constant but still larger than 1 (i.e. $n(a)\simeq n_{\ast }>1$) and, finally, when (iii) $a<a_{\ast }$ the refractive index is truly dynamical since $n(a)\simeq n_{\ast }{(a∕a_{\ast })}^{\alpha }$. When $a<a_{\ast }$ we can directly compute $b(\eta )$ in case the dynamics of the geometry is given by a conventional inflationary stage where $aH=-1∕[(1-ϵ)\tau ]$. By direct integration from $n(\eta )d\eta =d\tau$ we obtain the relation between $\eta$ and $\tau$ and then compute $b(\eta )=a(\eta )∕\sqrt{n(\eta )}$; the result is In the $\eta$-time coordinate the evolution of the tensor modes can be quantized in the standard manner as where, as usual, the sum over $\lambda$ runs over the two tensor polarizations. From Eq. (5.26), the mode functions obey ${\ddot{f}}_{k,\lambda }+(k^2-\ddot{b}∕b)f_{k\lambda }=0$ and recalling the expression of Eq. (5.28), for each tensor polarization the solution of the mode function reads where $|\mathcal{𝒩}|=\sqrt{\pi ∕2}$ and $H_{\nu }^{(1)}(-k\eta )$ is the Hankel function of first kind.^78,79 Once more, the solution of Eq. (5.30) is relatively simple in terms of the $\eta$-time but gets more cumbersome in the conformal time coordinate. While in the case of the concordance paradigm the low-frequency spectral index is red (i.e. $n_T^{(\mathrm{\text{low}})}<0$) when the refractive index is dynamical $n_T^{(\mathrm{\text{low}})}>0$. Indeed if we compute the tensor power spectrum in the long wavelength limit (i.e. $|k\eta |<1$) from Eq. (5.30) we obtain where ${\overline{\mathcal{\mathscr{H}}}}_{\ast }^2=1∕{\eta }_{\ast }^2$ and $n_T^{(\mathrm{\text{low}})}=3-2\nu =2(1-\zeta )$ is, by definition, the low-frequency spectral index which is generically blue as long as $\alpha >0$ in Eq. (5.27) : The spectral energy density for typical wavenumbers $k<a_{\ast }{H}_{\ast }$ can also be computed once the expression of $b(\eta )$ is known; within the WKB approach already outlined in Sec. 3 is given by Equation (5.33) has been computed by assuming that the reentry of the relevant wavelengths occurs when the refractive index is not dynamical and this implies that when the relevant wavelength reenters the $\eta$-time and the conformal time coordinates coincide, i.e. ${\eta }_{\mathrm{\text{re}}}={\tau }_{\mathrm{\text{re}}}$. Furthermore in the simplest situation ${\tau }_{\mathrm{\text{re}}}$ falls within the radiation phase (i.e. $a^{\prime \prime }\to 0$) so that $k{\tau }_{\mathrm{\text{re}}}\ll 1$ in Eq. (5.33). Since any wavelength exiting for $\eta <-{\eta }_{\ast }$ does its first crossing during the inflationary phase, the corresponding refractive index is $n=n_{\ast }{(a∕a_{\ast })}^{\alpha }$; the explicit expression of ${\mathcal{𝒬}}_k({\eta }_{\mathrm{\text{ex}}},{\eta }_{\mathrm{\text{re}}})$ is which is the analog of the expression already obtained in Sec. 3 when the refractive index is not dynamical. Finally, using Eq. (5.28) an even more explicit expression of the spectral energy density can be deduced:

$$h_0^2{\mathrm{\Omega }}_{gw}(\nu ,{\tau }_0)={\left(\frac{H_1}{M_P}\right)}^2{\mathcal{𝒟}}_{\ast }(\alpha ,n_T^{(\mathrm{\text{low}})}){\left(\frac{\nu }{{\nu }_{\ast }}\right)}^{n_T^{(\mathrm{\text{low}})}},{\nu }_{\mathrm{\text{eq}}}<\nu <{\nu }_{\ast },$$(5.35)

$${\mathcal{𝒟}}_{\ast }(\alpha ,n_T)=\frac{4n_{\ast }^3{h}_0^2{\mathrm{\Omega }}_{R0}}{3\pi }{\left(1+\frac{\alpha }{1-{ϵ}_k}\right)}^2\left(\frac{g_{\rho ,r}}{g_{\rho ,\mathrm{\text{eq}}}}\right){\left(\frac{g_{s,\mathrm{\text{eq}}}}{g_{s,r}}\right)}^{4∕3}{\left(\frac{{\mathrm{\Omega }}_{M0}}{{\mathrm{\Omega }}_{\mathrm{\Lambda }}}\right)}^2.$$(5.36)

As usual ${\mathrm{\Omega }}_{M0}$ and ${\mathrm{\Omega }}_{\mathrm{\Lambda }}$ denote the present critical fractions of matter and dark energy; it is actually well known that the dominance of dark energy suppresses the spectrum by a factor ${({\mathrm{\Omega }}_{M0}∕{\mathrm{\Omega }}_{\mathrm{\Lambda }})}^2=\mathcal{𝒪}(0.1)$ (see, for instance, Ref. 37). In Eq. (5.35), ${\nu }_{\ast }$ denotes the frequency of the spectrum associated with ${\eta }_{\ast }$ and since $k_{\ast }=1∕{\eta }_{\ast }$ the corresponding comoving frequency is In Eq. (5.37), $N_{\ast }=ln(a_{\ast }∕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, as before, ${\overline{\nu }}_{max}$ indicates the maximal frequency of the spectrum and it coincides with Eq. (3.54) since, so far, the radiation dominance starts right after the end of inflation. The tensor spectral index of Eq. (5.32) applies in the low and intermediate frequency ranges when the corresponding wavelengths exit during inflation and reenter in the radiation phase; in Eq. (5.32) $\alpha$ is always much larger than ${ϵ}_k\simeq r_T∕16\le 0.03∕16\ll 1$ so that the exact result can be accurately evaluated in the limit ${ϵ}_k\ll 1$. While Eqs. (5.32) and (5.35) hold for $\nu <{\nu }_{\ast }$, the spectral energy density can also be evaluated in the range ${\nu }_{\ast }<\nu <{\overline{\nu }}_{max}$ (i.e. $a_{\ast }{H}_{\ast }<k<a_1{H}_1$) corresponding to wavelengths that exited the comoving horizon when the refractive index was already constant (i.e. $n\to n_{\ast }$ and ${\eta }_{\mathrm{\text{ex}}}={\tau }_{\ast }∕n_{\ast }$); in this case the spectral energy density becomes where the spectral index is given by $n_T^{(\mathrm{\text{high}})}=-2{ϵ}_k=-r_T∕8$ and

Equation (5.38) evaluated for $\nu ={\nu }_{\ast }$ reproduces Eq. (5.35) computed at the same reference frequency and the equivalence of the two expressions ultimately follows from Eq. (5.37). Furthermore, in Eqs. (5.35) and (5.32) ${(H_1∕M_P)}^2$ can be traded for $\pi {ϵ}_k{\mathcal{𝒜}}_{\mathcal{ℛ}}$ where ${\mathcal{𝒜}}_{\mathcal{ℛ}}$ is the amplitude of curvature inhomogeneities at the pivot scale $k_p$. It is finally worth recalling that, for a standard thermal history, $g_{s,\mathrm{\text{eq}}}=3.94$ while $g_{\rho ,r}=g_{s,r}=106.75$ in Eqs. (5.36)–(5.39). In Eqs. (5.37) and (5.39), $N_{\ast }$ measures the range of variation of the refractive index during inflation and, for this reason, $N_{\ast }<N_t$. As we shall see in a moment, the relatively short inflationary stages (where $N_t\le \mathcal{𝒪}(61)$) seem to be preferred for a potential explanations of the PTA excesses. Equations (5.33), (5.35) and (5.38) are now compared with the parametrizations of the PTA signal given in Eqs. (5.13) and (5.14). Since, by definition, the intermediate spectral index is given as $2+2\beta =n_T^{(\mathrm{\text{low}})}$ Eq. (5.32) implies a relation that determines $\alpha$ as a function of ${ϵ}_k$ (or $r_T$) and $\beta$ :

Moreover, given that $q_0$ depends on all the other parameters determining the amplitude of ${\mathrm{\Omega }}_{gw}(\nu ,{\tau }_0)$ (see Eqs. (5.35) and (5.38)), we can demand that $\beta$ and $q_0$ fall within the phenomenologically allowed ranges and check if the results of Eqs. (5.35)–(5.38) are compatible with the empirical determinations of the PTA. According to the PPTA the values of $\beta$ and $q_0$ fall, respectively, in the following intervals : Equation (5.41) constrains the spectral energy density and the corresponding region of the theoretical parameters is illustrated in the left plot of figure where we report $q_0(\beta ,N_{\ast })$ for different values of $N_t$; the shape of each shaded region directly follows by requiring $2.2<q_0(\beta ,N_{\ast })<4.4$ for the various $N_t$ mentioned in the plot. On a technical side we note that Eq. (5.40) has been used with the purpose of trading directly $\alpha$ for $\beta$ at a fixed value of ${ϵ}_k$. The same analysis illustrated in the case of the PPTA can be repeated for the NANOgrav determinations with slightly different results; the analog of Eq. (5.41) is now given by³⁹ While the range of $\beta$ given in Eq. (5.42) is narrower than in Eq. (5.41), in the case of $q_0$ we observe the opposite: the allowed values of $q_0$ of Eq. (5.42) are comparatively larger than the ones of Eq. (5.41). A second class of constraints determining the shaded allowed regions is related to the direct bounds from the operating wide-band detectors; in particular we remind that the LIGO, Virgo and Kagra (LVK) collaborations reported a constraint implying^{33,34,35,36,37} : in the case of a flat spectral energy density; in the present notations ${\nu }_L$ indicates the LVK frequency. The limit of Eq. (5.43) improves on a series of bounds previously deduced by the wide-band interferometers (see Ref. 37 for a review of the older results). In the present notations the parametrization of ${\mathrm{\Omega }}_{gw}(\nu ,{\tau }_0)$ adopted by Ref. 36 reads

and the three specific cases constrained in Refs. 35 and 36 are reminded in Table 1. As the value of $\sigma$ increases the bound becomes more restrictive for a fixed reference frequency and the three previous results are summarized by the following interpolating formula :

Since in the present case the bound (5.45) should be applied at high frequencies we will have $\sigma =-2{ϵ}_k∕(1-{ϵ}_k)$ with $ϵ\ll 0.1$; to leading order in ${ϵ}_k$, Eq. (5.45) implies that $log\overline{\mathrm{\Omega }}({ϵ}_k)<-8.236-0.335{ϵ}_k-0.393{ϵ}_k^2$. As an example in the two plots of Fig. 10 we considered two different values of $\beta$ (i.e. $\beta =-0.65$ and $\beta =-0.55$). If $N_{\ast }$ 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 BBN constraint. Conversely, when $N_{\ast }<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 quasi-de Sitter stage of expansion. In both plots of Fig. 10, to ease the comparison, we selected $N_t=61$ while different values of $N_{\ast }$ are illustrated. In both plots, for the same choice of the parameters, we also illustrated (with an arrow) the PTA excess and the LVK bound. The PTA signal occurs for typical frequencies $\mathcal{𝒪}({\nu }_{\mathrm{\text{ref}}})$ while the LVK bound applies approximately between 25 and 100Hz.

The previous discussion does not exclude the possibility of two concurrent modifications of the comoving horizon operating before and after the end of inflation. This viewpoint is explored in Fig. 11 where we consider the possibility that the refractive index stops its evolution well before the end of inflation (i.e. $N_{\ast }\ll N_t$); The spectral energy density in critical units will therefore have three different slopes for $\nu >{\nu }_{\mathrm{\text{eq}}}$. In both plots of Fig. 11, at intermediate frequencies $h_0^2{\mathrm{\Omega }}_{gw}(\nu ,{\tau }_0)$ has the same intermediate slopes appearing in Fig. 10 (see also Eqs. (5.32) and (5.40)). 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. 11 the high-frequency spectral slope is $\mathcal{𝒪}(1)$ since during the post-inflationary stage the evolution is described by a stiff fluid with $\delta \simeq 1∕2$ implying that ${(aH)}^{-1}\propto a^2$. The main difference between the plots of Figs. 10 and 11 comes from the high-frequency shape where the bounds coming from BBN must be taken into account (see Eq. (5.3) and discussion therein). The theoretical perspective explored in this discussion 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 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 post-inflationary expansion rate prior to BBN. This possibility can be safely ruled out. 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. This explanation seems viable in the light of the current observations. We may 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. This option may also be viable with some caveats and has been explored in Ref. 178.

5.3. Space-borne interferometers and the expansion history

The direct measurements in the range ${\nu }_{\mathrm{\text{ew}}}\le \nu <\mathrm{\text{Hz}}$ may primarily clarify the nature of the post-inflationary expansion rate. Indeed, after inflation, the expanding stage could include a sequence of stages expanding either faster or slower than radiation; in this situation a hump in $h_0^2{\mathrm{\Omega }}_{gw}(\nu ,{\tau }_0)$ is generically expected below the a fraction of the Hz where the relic gravitons may exceed (even by eight orders of magnitude) the signals obtained under the hypothesis of radiation dominance throughout the whole expansion history prior to the formation of light nuclei.

5.3.1. The conventional wisdom

An old and conventional viewpoint stipulates that between a fraction of the mHz and few Hz the spectral energy density of the inflationary gravitons can be disregarded even assuming the most optimistic sensitivities of the space-borne detectors. On the contrary, always within the standard lore, in the region between the $\mu$Hz and few Hz the signals coming from the electroweak physics (or from some other phase transition) should represent the dominant contribution of cosmological origin. This perspective is not completely consistent for (at least) two independent reasons.

•

The first one (already mentioned earlier on in this section) is that the electroweak phase transition does not proceed through the formation of bubbles of the new phase and does not imply large deviations from homogeneity as required for the formation of a diffuse secondary background of gravitational radiation. This statement hods given the measured values of the Higgs and W masses.

•

The usual counterargument is that we might expect strongly first-order phase transitions from new physics which did not show up so far from collider searches. This assumption is however ad hoc since there are no tangible signals of a new electroweak physics from colliders; it is therefore not clear why the purported new physics should always lead to a burst of gravitational radiation in a range compatible with ${\nu }_{\mathrm{\text{ew}}}$.

We are now going to discuss how a modified expansion history may lead to a hump in the frequency domain compatible with ${\nu }_{\mathrm{\text{ew}}}$. This is why any limit on the spectral energy density of the relic gravitons between few $\mu$Hz and the Hz indirectly constrains the timeline of the post-inflationary expansion rate.

5.3.2. Chirp amplitudes and frequency dependence

The direct bounds on the relic gravitons from the audio band ultimately depend upon the spectrum of the signal. For a nearly scale-invariant spectrum, $h_0^2{\mathrm{\Omega }}_{gw}(\nu ,{\tau }_0)<\mathcal{𝒪}(10^{-9})$ between 10Hz and 80Hz^33,34,35,36 (see also Ref. 37 for a recent review including earlier bounds). To compare the ground-based detectors and the space-borne interferometers it is useful to express the spectral energy density in terms of the chirp amplitude $h_c(\nu ,{\tau }_0)$³⁷ when the typical frequencies fall in the audio band :

From left to right Eq. (5.46) implies that to probe $h_0^2{\mathrm{\Omega }}_{gw}(\nu ,{\tau }_0)=\mathcal{𝒪}(10^{-9})$ we should have a sensitivity in the chirp amplitude $\mathcal{𝒪}(10^{-24})$ for a typical frequency $\nu =\mathcal{𝒪}(100)$Hz. From right to left the same relation suggests instead that, for comparable sensitivities in $h_c(\nu ,{\tau }_0)$, the minimal detectable $h_0^2{\mathrm{\Omega }}_{gw}(\nu ,{\tau }_0)$ gets comparatively smaller with the frequency; besides the absence of seismic noise this is probably one of strongest arguments in favor of space-borne detectors for typical frequencies ranging between a fraction of the mHz and the Hz. This is why the minimal detectable spectral energy density could be $h_0^2{\mathrm{\Omega }}_{gw}(\nu ,{\tau }_0)=\mathcal{𝒪}(10^{-11})$ or even $h_0^2{\mathrm{\Omega }}_{gw}(\nu ,{\tau }_0)=\mathcal{𝒪}(10^{-15})$ under the hypothesis that the same sensitivity reached in the audio band for the chirp amplitude can also be achieved in the mHz range. With this great hope, various space-borne detectors have been proposed so far: the Laser Interferometric Space Antenna (LISA),^196,197 the Deci-Hertz Interferometer Gravitational Wave Observatory (DECIGO),^198,199 the Ultimate-DECIGO²⁰⁰ (conventionally referred to as U-DECIGO), the Big Bang Observer (BBO).²⁰¹ This list has been recently enriched by the Taiji^202,203 and by the TianQin^204,205 experiments. Since these instruments are not yet operational (but might come into operation within the next 20 years) their actual sensitivities are difficult to assess, at the moment. However, without dwelling on the specific nature of the noise power spectra, Eq. (5.46) shows that, as long as $h_c=\mathcal{𝒪}(10^{-23})$ the space-borne detectors might probe $h_0^2{\mathrm{\Omega }}_{gw}(\nu ,{\tau }_0)=\mathcal{𝒪}(10^{-14})$ for ${\nu }_S=\mathcal{𝒪}(0.01)$Hz and this is, roughly speaking, the daring expectation of DECIGO^198,199 and of U-DECIGO.²⁰⁰

The fiducial frequency interval of space-borne interferometers ranges from a fraction of the mHz to the Hz and, within this interval, the minimal detectable spectral energy density (denoted hereunder by $h_0^2{\mathrm{\Omega }}_{gw}^{(\mathrm{\text{min}})}(\nu ,{\tau }_0)$) defines the potential sensitivity of the hypothetical instrument. The LISA interferometers might hopefully probe the following region of the parameter space :

In the case of the DECIGO the minimal detectable spectral energy density could be smaller The values of Eq. (5.48) are still quite hypothetical so that it is prudent to choose $h_0^2{\mathrm{\Omega }}_{gw}^{(\mathrm{\text{min}})}(\nu ,{\tau }_0)$ between the standard values of the hoped sensitivity of the DECIGO project (suggesting $h_0^2{\mathrm{\Omega }}_{gw}^{(\mathrm{\text{min}})}(\nu ,{\tau }_0)=\mathcal{𝒪}(10^{-13.1})$) and the optimistic figure reachable by the Ultimate-DECIGO²⁰⁰ (conventionally referred to as U-DECIGO) where $h_0^2{\mathrm{\Omega }}_{gw}^{(\mathrm{\text{min}})}(\nu ,{\tau }_0)=\mathcal{𝒪}(10^{-17.5})$. For the record, the BBO²⁰¹ might reach sensitivities There finally exist also recent proposals such as Taiji^202,203 and TianQin^204,205 leading to figures that are roughly comparable with the LISA values. In summary for the typical frequency of the space-borne detectors we consider the broad range $0.1\mathrm{\text{mHz}}<{\nu }_S<0.1\mathrm{\text{Hz}}$ and suppose that in this range $h_0^2{\mathrm{\Omega }}_{gw}^{(\mathrm{\text{min}})}(\nu ,{\tau }_0)$ may take the following two extreme values : While the two values of Eq. (5.50) are both quite optimistic, they are customarily assumed by the observational proposals and, for this reason, they are used here only for illustration.

5.3.3. Humps in the spectra from the modified expansion rate

The expansion rates can be bounded by requiring that for frequencies of the order of ${\nu }_S$ the corresponding spectral energy density exceeds $h_0^2{\mathrm{\Omega }}_{gw}^{(\mathrm{\text{min}})}({\nu }_S,{\tau }_0)$; all the other constraints on the diffuse backgrounds of gravitational radiation must also be satisfied. With this strategy it is possible to constrain the unconventional post-inflationary expansion histories by simultaneously obtaining a large signal for frequencies $\mathcal{𝒪}({\nu }_S)$. The spectral energy density of the relic gravitons might exhibit various successive local maxima but the simplest case consists in a single hump for frequencies comparable with ${\nu }_S$. The $h_0^2{\mathrm{\Omega }}_{gw}(\nu ,{\tau }_0)$ can be expressed in this case as

where, as usual, $n_T^{(\mathrm{\text{low}})}$ is the spectral index associated with the wavelengths leaving the Hubble radius during the inflationary phase and reentering during the radiation stage. In Eq. (5.51), ${\nu }_p$ and ${\nu }_{\mathrm{\text{eq}}}$ define the lowest frequency range of the spectral energy density : The transfer function of Eq. (5.51) also includes the dependence on the spectral slopes $n_T^{(1)}$ and $n_T^{(2)}$; up to corrections $\mathcal{𝒪}(r_T)$ the depend directly on the expansion rate expressed in the conformal time parametrization : where $n_T^{(1)}<0$ and $n_T^{(2)}>0$ since we consider the situation where during the first stage the Universe expands faster than radiation (i.e. ${\delta }_1>1$) while in the second stage it is slower than radiation (i.e. ${\delta }_2<1$). In the simplest case where the consistency relations are enforced we have In Eq. (5.52), ${\mathcal{𝒯}}_{\mathrm{\text{low}}}(\nu ∕{\nu }_{\mathrm{\text{eq}}})$ is the low-frequency transfer function of the spectral energy density³⁷ : The high-frequency transfer function ${\mathcal{𝒯}}_{\mathrm{\text{high}}}(\nu ,{\nu }_2,{\nu }_r,{\delta }_1,{\delta }_2)$ appearing in Eq. (5.51) is specifically discussed in Ref. 206. In Fig. 12, the spectral energy density has been explicitly illustrated for a selection of the parameters. In the left plot of Fig. 12, we selected $\xi =10^{-36}$ and ${\xi }_2=10^{-10}$ for different values of ${\delta }_1>1$ and ${\delta }_2<1$. We recall that, by definition,

As expected the value of ${\nu }_r$ is always larger than $10^{-10}$. The parameters of the dot-dashed and of the dashed curves of the left plot in Fig. 12 have been selected in order to get an artificially large signal that is in fact excluded both by the BBN constraint and by the limit of ground-based detectors. The results of the right plot in Fig. 12 correspond instead to a slightly different choice of the parameters, namely, $\xi =10^{-34}$ and ${\xi }_2=10^{-8}$. For illustration we have chosen ${\delta }_1\to 1$ implying that between ${\nu }_{max}$ and ${\nu }_2$ the spectral energy density is quasi-flat. This is the most constraining case from the viewpoint of the limits coming from wide-band detectors.^{33,34,35,36,37}

5.3.4. Complementary considerations

So far we saw that different frameworks motivate the presence of post-inflationary stages expanding at rate either faster or slower than radiation and this is why the model independent perspective of Ref. 45 (see also Refs. 107 and 108) is, in our opinion, the most useful. We remind here that stiff post-inflationary phase is dynamically realized in different situations and the first speculations along this direction probably date back to the ideas Zeldovich,¹⁴⁹ Sakharov⁷¹ and Grishchuk.¹⁵ After the formulation of conventional inflationary models Ford¹⁵⁰ noted that gravitational particle production at the end of inflation could account for the entropy of the present Universe and observed that the backreaction effects of the created quanta constrain the length of a stiff post-inflationary phase by making the expansion dominated by radiation. These effects typically lead, in our notations, to a pivotal frequency ${\nu }_r$ of the order of the mHz. It has been later argued by Spokoiny¹⁵¹ that various classes of scalar field potentials exhibit a transition from inflation to a stiff phase dominated by the kinetic energy of the inflaton. In more recent times it became increasingly plausible to have a single scalar field acting as inflaton in the early Universe and as quintessence field in the late Universe.^152,153 A generic signature of a post-inflationary phase stiffer than radiation is the production of relic gravitons with increasing spectral energy density.⁴⁵ In quintessential inflationary models the inflaton and the quintessence field are identified in a single scalar degree of freedom⁴⁸ and various concrete forms of the inflaton-quintessence potential $V(\phi )$ have been proposed and scrutinized through the years. The transition between an inflationary phase and a kinetic phase can be realized both with power-law potentials and with exponential potentials. See also Refs. 154, 155, 156 for further applications. We pointed out so far that the expected signal coming from the phase transitions is probably rather small; however, as suggested in the past, a strong hypermagnetic background may be present in the symmetric phase of the electroweak theory^157,158,159 because of the symmetries of the plasma at finite density and finite conductivity. The overall magnitude of the spectra of gravitational radiation induced by a hypermagnetic background have been estimated, for the first time, in Refs. 157, 158 and turn out to be generally different from the ones associated with a modified post-inflationary evolution.¹⁵⁹

6. The Expansion History and the High-Frequency Gravitons

The high-frequency region of the spectrum ranges between few Hz and the THz since, as already discussed in section 3, their maximal frequency cannot exceed the $\mathrm{\text{THz}}$ domain. Only for practical reasons, in this broad region we distinguish the ultra-high frequency domain (between the MHz and the THz) and the high-frequency band ranging from the Hz to the MHz. To analyze the bounds on the post-inflationary expansion rate it is simpler to address first the THz domain and then focus on the MHz region that also includes the operating window of ground-based interferometers.

6.1. Spikes in the GHz domain

If the post-inflationary evolution consists of a single stage, the results of Eqs. (3.51) and (3.77), (3.78) suggest that the maximal signal should always be concentrated between the GHz and the THz. This happens when the expansion rate is slower than radiation (i.e. $\delta <1$, see Sec. 2 and notations therein). If the expansion rate is instead faster than radiation (i.e. $\delta >1$) the high-frequency slope is negative (i.e. $n_T^{(\mathrm{\text{high}})}<0$) so that the spectral energy density is ultimately decreasing and potentially even smaller than the signal of the concordance paradigm (i.e. $h_0^2{\mathrm{\Omega }}_{gw}(\nu ,{\tau }_0)\le \mathcal{𝒪}(10^{-17})$) in the same range of frequencies.

6.1.1. General considerations

The high-frequency branch of the spectrum bears the mark of the post-inflationary expansion rate and from the frequency profile of the spectral energy density we can directly infer the post-inflationary expansion rate, the maximal frequency and the other pivotal frequencies of the spectrum (including the approximate curvature scale of radiation dominance). In the left plot of Fig. 13, we report the spectral energy density in critical units as a function of the frequency for a selection of examples (common logarithms are employed on both axes); note that, for the reported spectra, the post-inflationary expansion rate is slower than radiation. In the right plot of the same figure the parameter space is illustrated in the plane $(log\xi ,\delta )$ where $\xi =H_r∕H_1$ estimates the overall duration of the post-inflationary stage of expansion. The shaded region in the right plot of Fig. 13 denotes instead the allowed portion of the parameter space; in particular, the darker sector above the dashed curve accounts for the bounds coming from BBN which are ultimately the most constraining. The region with lighter shading below the dashed curve corresponds to the current limits coming from wide-band interferometers. Finally the dashed curve itself is deduced through a semianalytic approximation discussed hereunder. The specific features of ${\mathrm{\Omega }}_{gw}(\nu ,{\tau }_0)$ illustrated Fig. 13 allow for a quantitative reconstruction of the expansion rate if and when the sensitivities of the dedicated detectors (both in the audio band and in the GHz region) will be able to resolve the class of signals suggested here. Along this perspective the main features of Fig. 13 motivate, in short, the following observations.

•	In the aHz region the spectral energy density decreases as ${\nu }^{-2}$ while we can appreciate the suppression due to the neutrino free streaming close to ${\nu }_{\mathrm{\text{bbn}}}$.160,161,162,163,164 Other sources of suppression taken into account in Fig. 13 and in the remaining plots include the late-time dominance of dark energy and the evolution of relativistic species. The spectra of Fig. 13 have been deduced by using for the fiducial parameters the last Planck data release in the case of three massless neutrinos where $R_{\nu }={\rho }_{\nu }∕({\rho }_{\gamma }+{\rho }_{\nu })=0.405$, as indicated on top of each plots; this is the choice of the minimal $\mathrm{\Lambda }$CDM scenario. If the radiation would dominate the whole post-inflationary evolution the quasi-flat plateau (decreasing because of the slow-roll corrections) would last up to frequencies $\mathcal{𝒪}(300)$MHz.
•	When the expansion rate is faster than radiation (i.e. $\delta >1$ in the notations of Fig. 13) the spectral energy density further decreases between ${\nu }_r$ and ${\nu }_{max}$: this timeline implies that $h_0^2{\mathrm{\Omega }}_{gw}(\nu ,{\tau }_0)\ll \mathcal{𝒪}(10^{-17})$ (in particular in the audio band). No further constraints (besides the low-frequency limits that translate into the upper bound on $r_T$42,43,44) appear when $\delta >1$.
•	When the post-inflationary expansion rate is slower than radiation (i.e. $\delta <1$ in Fig. 13) the spectral energy density grows for $\nu >{\nu }_r$ and eventually reaches a maximum that roughly corresponds to the onset of the exponential suppression taking place for $\nu >{\nu }_{max}$.

To trace the origin of the high-frequency spike we remark that $h_0^2{\mathrm{\Omega }}_{gw}(\nu ,{\tau }_0)$ can be written, with compact notations, as

where ${\nu }_p$ and ${\nu }_{\mathrm{\text{eq}}}$ are, respectively, the pivot and the equality frequencies already introduced in Eq. (5.52). As usual $r_T=r_T({\nu }_p)$ is the tensor to scalar ratio evaluated at the pivot scale whereas ${\mathcal{𝒯}}_{\mathrm{\text{low}}}^2(\nu ∕{\nu }_{\mathrm{\text{eq}}})$ and ${\mathcal{𝒯}}_{\mathrm{\text{high}}}^2(\nu ∕{\nu }_r,\delta )$ denote the transfer functions directly computed for the spectral energy density^107,108; the value of ${\mathcal{𝒩}}_{\rho }=\mathcal{𝒪}(4)\times 10^{-15}$ can be accurately fixed both analytically and numerically.^ff Since for $\nu >{\nu }_r$ the high-energy transfer function has the slope $n_T^{(\mathrm{\text{high}})}$ (i.e. ${\mathcal{𝒯}}_{\mathrm{\text{high}}}^2\to {(\nu ∕{\nu }_r)}^{n_T^{(\mathrm{\text{high}})}}$) for the analytic estimates of the limits imposed on the spectral energy density we can express $h_0^2{\mathrm{\Omega }}_{gw}(\nu ,{\tau }_0)$ in the following approximate form : Equation (6.2) rests on the observation that ${\mathcal{𝒯}}_{\mathrm{\text{low}}}({\nu }_r∕{\nu }_{\mathrm{\text{eq}}})\to 1$ for $\nu \ge {\nu }_r$ while, in the same limit, it is also true that ${\overline{n}}_T\ll 1$. Since the overall normalization mildly depends on $\nu$ and $r_T$ we can express the spectral energy density as where ${\mathcal{𝒩}}_{\rho }(r_T,\nu )$ is

Although ${\mathcal{𝒩}}_{\rho }(r_T,\nu )$ exhibits a mild frequency dependence (mainly coming from neutrino free-streaming), for simplified analytic estimates this dependence can be approximately ignored. Along this perspective we may estimate ${\mathcal{𝒩}}_{\rho }=\mathcal{𝒪}(10^{-16.5})$ for $r_T\le \mathcal{𝒪}(0.06)$. In case a spectral energy density compatible with the one of Fig. 13 we may deduce various pieces of information on the early expansion rate and on the various transitions that occurred throughout the evolution of the plasma.

6.1.2. Invisible gravitons in the aHz region

The results of Fig. 13 do not rely on the specific value of $r_T$ and when $r_T\ll \mathcal{𝒪}(0.06)$ the high-frequency spike gets modified but does not disappear while the large-scale limits applicable to ${\mathrm{\Omega }}_{gw}(\nu ,{\tau }_0)$ in the aHz region also affect the small-scale constraint as suggested for the first time in Refs. 107 and 108 (see also the discussion of Sec. 4). In Fig. 14, we illustrate a sharp reduction of $r_T$ and a consequent suppression in the aHz region. When $r_T$ is reduced also the high-frequency signal gets suppressed although this effect is easily counterbalanced by a smaller value of $H_r$. To clarify this point we first observe that the values of $\xi =H_r∕H_1$ (illustrated both in Figs. 13 and 14) ultimately depend upon the assumed values of $r_T$: this happens since $H_1$ is sensitive to the inflationary expansion rate so that eventually $\xi$ scales as $\xi \propto r_T^{-1∕2}$ and it increases when $r_T$ gets progressively reduced. But ${\nu }_{max}={\xi }^{(\delta -1)∕[2(\delta +1)]}{\overline{\nu }}_{max}$ depends on $\xi$ also because ${\overline{\nu }}_{max}$ itself scales as $r_T^{1∕4}$ (see Eqs. (3.54) and (3.55) and discussion therein). These different effects can be combined with the purpose of deducing the scaling of $h_0^2{\mathrm{\Omega }}_{gw}({\nu }_{max},{\tau }_0)$ with $r_T$; up to a numerical factor that depends on $\delta$ the result is

where $\mathcal{ℬ}(\delta )={\mathcal{𝒞}}^4(g_{\rho },g_s,{\tau }_r,{\tau }_{\mathrm{\text{eq}}}){(16∕\pi )}^{(\delta -1)∕(\delta +1)}∕3$ is just a numerical factor that is not strictly essential in the forthcoming considerations. According to Eq. (6.5), for the same $H_r$ a reduction of $r_T$ entails an overall suppression of $h_0^2{\mathrm{\Omega }}_{gw}({\nu }_{max},{\tau }_0)$. Conversely, when $r_T$ is kept fixed, a reduction of $H_r$ increases $h_0^2{\mathrm{\Omega }}_{gw}({\nu }_{max},{\tau }_0)$ when $\delta <1$; when $\delta >1$ a reduction of $H_r$ (for fixed $r_T$) further suppresses the spectral energy density. This means, as anticipated, that a reduction of $r_T$ may be compensated by an appropriate reduction of $H_r$ in the case when the post-inflationary expansion rate is slower than radiation. In the left plot of Fig. 14 the high-frequency spectral indices have been chosen exactly with the purpose of demonstrating that lower values of $r_T$ do not necessarily suppress the high-frequency signal that remains exceedingly large in comparison with $\mathcal{𝒪}(10^{-17})$. In the right plot of Fig. 14, we illustrate the parameter space in the plane defined by $H_r$ and $\delta$. As in Fig. 13 the shaded area corresponds to the region allowed by the constraints stemming from BBN while the darker region comes from the bounds of the wide-band detectors operating in the audio band. In the right plot Fig. 14 $r_T=3\times 10^{-2}$ while in the two plots of Fig. 15 we selected instead $r_T=3\times 10^{-4}$ and $r_T=3\times 10^{-6}$, respectively. The shaded regions in Fig. 14 illustrate the intersection between the BBN bounds and the limits following from wide-band interferometers. Indeed, a general requirement determining the lowest value of $r_T$ is obtained from the current limits (summarized in Table 1) on the presence of relic graviton backgrounds in the audio band.^{33,34,35,36,37} By following here this approach we adopted the condition

where ${\nu }_{\mathrm{\text{LVK}}}$ denotes the Ligo–Virgo–Kagra frequency which can be estimated in terms of ${\nu }_{\mathrm{\text{ref}}}$. The most sensitive region for the detection of relic gravitons in the audio band is, grossly speaking, below $0.1$ kHz since, in this band, the overlap reduction function has its first zero.³⁷ Equation (73) requires, in practice, that the bounds coming from wide-band interferometers are satisfied while, in the same frequency range, $h_0^2{\mathrm{\Omega }}_{gw}(\nu ,{\tau }_0)$ is larger than $\mathcal{𝒪}(10^{-13})$. We cannot foresee when the corresponding sensitivity will be reached by wide-band detectors but the requirement of Eq. (73) follows from some of the optimistic claims suggested by the observational collaborations.^gg36

6.1.3. Bounds on the expansion rate

In terms of Eqs. (6.3) and (6.4) the BBN constraint assumes a particularly simple analytical form and since the largest contribution to the integral comes from the bunch of frequencies $\mathcal{𝒪}({\nu }_{max})$, Eqs. (6.3) and (6.4) can be used to set a limit on the integral of Eq. (3.47); if we require, for instance, $h_0^2{\mathrm{\Omega }}_{gw}({\nu }_{max},{\tau }_0)<10^{-6}$ we obtain the following constraint in the $(\xi ,\delta )$ plane :

where we used the scaling of $({\nu }_{max}∕{\nu }_r)$ with $\xi$, i.e. $({\nu }_{max}∕{\nu }_r)\propto {\xi }^{-1∕(\delta +1)}$. As in Eq. (6.7) it is always true that $r_T\ll \delta$, Eq. (6.7) translates into $log\xi >-5.25(1+\delta )∕(1-\delta )$ (where we took $r_T=0.06$ and consequently estimated $log{\mathcal{𝒩}}_{\rho }=-16.5$). Equation (6.7) implies then a lower bound on $\xi$. Indeed it can be argued that $\delta$ cannot get smaller than $1∕2$ (see below Eq. (6.19) and discussion thereafter) and, in this case, we would have $log\xi >-15.75$. If $\delta$ would decrease below $1∕2$ the lower bound on $\xi$ would get larger: when $\delta =1∕3$ the lower bound is given by $\xi >10^{-10.5}$, and so on. We finally remind, as already pointed out in Sec. 5, that a further lower bound on $\xi$ is obtained by requiring that ${\nu }_r>{\nu }_{\mathrm{\text{bbn}}}$; but then the bound is much less restrictive and it only demands $\xi >10^{-38}$.

The limits obtained from Eqs. (6.3), (6.4) and (6.7) can be checked by direct numerical evaluation of the integral appearing in Eq. (3.47). In the right plot of Fig. 13 the shaded region illustrates the BBN constraint directly computed from Eq. (3.47) and, in the same plot, the dashed curve describes the analytic bound coming from Eq. (6.7) for $r_T\to 0.03$. The two determinations compare quite well and corroborate the approximation schemes of Eqs. (6.3) and (6.4). We point out that in the right plot of Fig. 13 the darker region corresponds to the BBN whereas the area defined by the lighter shading accounts for the LVK bounds of Table 1. The limits illustrated in Figs. 13, 14 and 15 are two-dimensional slices of a three-dimensional parameter space where the values of $r_T$ are consistently reduced. The allowed region in three dimensions is represented by volume in the space ($\delta$, $r_T$, $H_r$). To deduce the three-dimensional bounds we first observe, once more, that in spite of the complicated expansion timeline, the frequency ${\nu }_r$ is always related to ${\overline{\nu }}_{max}$ as ${\nu }_r=\sqrt{\xi }{\overline{\nu }}_{max}$. If we now require

we obtain, in practice, that $H_r\ge H_{\mathrm{\text{bbn}}}$. Recalling the considerations of Eqs. (3.45) and (3.46) the simplest way of obtaining a bound on the expansion rate is to appreciate since, by definition, $\overline{n}({\nu }_{max},{\tau }_0)=\mathcal{𝒪}(1)$. From Eq. (6.9), we can write Because we are always requiring that ${\nu }_r\ge {\nu }_{\mathrm{\text{bbn}}}$ the integral of Eq. (3.47) can be approximated as follows : so that the BBN bound is now compactly expressed as A single post-inflationary stage expanding at a rate slower than radiation has fewer parameters in comparison with multiple stages of expansion (see e.g. Figs. 4 and 6). In the simplest situation of a single post-inflationary stage, the previous discussion clarifies that the three relevant parameters are: the tensor to scalar ratio $r_T$, the expansion rate during the post-inflationary evolution (related to $\delta$) and the Hubble rate at the onset of the radiation stage (i.e. $H_r∕M_P$). We can eventually trade $(H_r∕M_P)$ for the reheating temperature; the relation between the two quantities is obtained by assuming complete thermal equilibrium at $T_r$ and it is The full three-dimensional parameter space is illustrated in Fig. 16: if the parameters fall within the shaded volume of Fig. 16 all the relevant constraints are satisfied. The illustrative examples reported in Figs. 15 and 17 can be viewed as two-dimensional projections of the three-dimensional parameter space of Fig. 16. From the shape of the spectral energy density it is then possible to infer the post-inflationary expansion rate and for a single post-inflationary phase the maximum of $h_0^2{\mathrm{\Omega }}_{gw}(\nu ,{\tau }_0)$ falls in the GHz region. If the expansion rate is more complicated the maximum can be present from the GHz region to the audio band and this is the possibility examined in the following subsection.

6.2. Spikes in the kHz domain

When a single post-inflationary stage precedes the radiation epoch, $h_0^2{\mathrm{\Omega }}_{gw}(\nu ,{\tau }_0)$ consists of three separate branches. If the timeline of the expansion rate contains different stages of expansion the spectral energy density may include multiple frequency domains and a maximum also develops below the MHz. In the simplest situation there are two intermediate stages preceding the radiation-dominated phase. Besides the standard aHz region and part of the intermediate branch (for ${\nu }_{\mathrm{\text{eq}}}<\nu <{\nu }_r$), the slopes in the two supplementary ranges (i.e. ${\nu }_r<\nu <{\nu }_2$ and ${\nu }_2<\nu <{\nu }_{max}$) depend on the values of the expansion rates (i.e. ${\delta }_1$ and ${\delta }_2$) well before the electroweak epoch.

6.2.1. Maxima in the audio band

In Fig. 17, we illustrated few examples and the selected parameters also account for possible reductions of $r_T$. With a unified notation the spectral slopes (denoted in Fig. 17 by $n_1^{(\mathrm{\text{high}})}$ and $n_2^{(\mathrm{\text{high}})}$) are

The profiles of $h_0^2{\mathrm{\Omega }}_{gw}(\nu ,{\tau }_0)$ given in Fig. 17 follow from the shape of the comoving horizon where, prior to radiation dominance, the post-inflationary evolution consists of two successive stages where the background first expands faster than radiation (i.e. ${\delta }_1>1$) and then slows down (i.e. ${\delta }_2<1$). We have from Eq. (6.14) that the spectral energy density decreases for $\nu >{\nu }_2$ (i.e. $n_1^{(\mathrm{\text{high}})}<0$) while it increases at lower frequencies (i.e. $n_2^{(\mathrm{\text{high}})}>0$ for $\nu <{\nu }_2$). If $r_T\ll 0.03$^42,43,44 Eq. (6.14) reduces to

and $n_1^{(\mathrm{\text{high}})}=2(1-{\delta }_1)<0$ for the wavelengths reentering before $a_2$ while for the wavelengths reentering between $a_2$ and $a_r$ we would have $n_2^{(\mathrm{\text{high}})}=2(1-{\delta }_2)>0$. In case the timeline is reversed (and ${\delta }_1<1$ while ${\delta }_2>2$) instead of a spike $h_0^2{\mathrm{\Omega }}_{gw}(\nu ,{\tau }_0)$ exhibits a trough but this timeline would be comparatively less constrained than the one of Fig. 17. All in all, recalling the parametrization of Eqs. (6.3) and (6.4), the two high-frequency branches of the spectral energy density can be parametrized as

$$h_0^2\mathrm{\Omega }(\nu ,{\tau }_0)={\mathcal{𝒩}}_{\rho }(r_T,\nu ){\left(\frac{\nu }{{\nu }_r}\right)}^{n_2^{(\mathrm{\text{high}})}},{\nu }_r<\nu <{\nu }_2,$$(6.16)

$$h_0^2\mathrm{\Omega }(\nu ,{\tau }_0)={\mathcal{𝒩}}_{\rho }(r_T,\nu ){\left(\frac{{\nu }_2}{{\nu }_r}\right)}^{n_2^{(\mathrm{\text{high}})}}{\left(\frac{\nu }{{\nu }_2}\right)}^{-|n_1^{(\mathrm{\text{high}})}|},{\nu }_2<\nu <{\nu }_{max},$$(6.17)

where we are implicitly assuming that $n_1^{(\mathrm{\text{high}})}<0$ and $n_2^{(\mathrm{\text{high}})}>0$. The spectral energy density given of Eqs. (6.16) and (6.17) exhibits a maximum for $\nu =\mathcal{𝒪}({\nu }_2$) but when ${\delta }_1\to 1$ the maximum is replaced by a plateau since $h_0^2{\mathrm{\Omega }}_{gw}(\nu ,{\tau }_0)$ flattens out (i.e. $n_1^{(\mathrm{\text{high}})}\to 0$ for $\nu >{\nu }_2$).⁶⁵ We then illustrated the situations that are phenomenologically more constraining; on this basis it is now possible to derive further limits on $r_T$ under the hypothesis of an expansion history including at least two different post-inflationary stages different from radiation (i.e. ${\delta }_i\ne 1$).

6.2.2. Again on the maximal frequency

The maximal frequency of the relic gravitons depends on $H_1$ but a modified post-inflationary evolution may artificially increase the value of ${\nu }_{max}$ by few orders of magnitude and potentially contradict the quantum bound of Eq. (3.51). The expansion histories leading to ${\nu }_{max}\gg \mathrm{\text{THz}}$ must then be rejected since the violations of the quantum bound also entail a violation of the limits set by BBN in the vicinity of ${\nu }_{max}$. To be more specific we now assume that between the end of inflation and the dominance of radiation there are n different stages of expansion that are arbitrarily different from radiation; this is, again, the general case illustrated in Figs. 4 and 6. We know from Eq. (3.52) that the value of the maximal frequency becomes, in this case :

where ${\xi }_i$ and ${\beta }_i=({\delta }_i-1)∕[2({\delta }_i+1)]$ measure, respectively, the duration of each of the post-inflationary stages and the corresponding expansion rate. When all the ${\beta }_i\to 0$ (i.e. ${\delta }_i\to 1$), the evolution is dominated by radiation from $H_1$ down to $H_{\mathrm{\text{eq}}}$ and ${\nu }_{max}\to {\overline{\nu }}_{max}=\mathcal{𝒪}(300)\mathrm{\text{MHz}}$. Conversely the value of ${\nu }_{max}$ given in Eq. (6.18) may exceed 300MHz provided at least one of the various ${\delta }_i$ gets smaller than 1. In the case of n intermediate stages preceding the dominance of radiation at $a_r$, between ${\nu }_{max}$ and ${\nu }_r$ there will be $n-2$ intermediate frequencies corresponding to specific breaks in the spectral energy density. The post-inflationary contribution to ${\nu }_{max}$ is then maximized when the ${\delta }_i$ and ${\xi }_i$ take their minimal values : The common value of the various ${\delta }_i$ corresponds to the slowest expansion rate of the primeval plasma. For instance in a perfect fluid the maximal value of the barotropic index (be it $w_{max}$) corresponds to the expansion rate, i.e. ${\delta }_{min}=2∕(3w_{max}+1)$ and since, at most, $w_{max}\to 1$ we obtain, as suggested in Eq. (6.19) ${\delta }_{min}\to 1∕2$. The expansion rate can also be slower than radiation when the energy–momentum tensor is dominated by the oscillations of the inflaton and if assume that the minimum of the potential is located in $\phi =0$$V(\phi )$, can be parametrized as $V(\phi )\simeq V_1{\mathrm{\Phi }}^{2q}$ (where, as usual, $\mathrm{\Phi }=\phi ∕M_P$). The averaged evolution of the comoving horizon can then mimic the timeline of a stiff epoch and the graviton spectra. Recalling Eqs. (4.24) and (4.25), during the coherent oscillations of $\phi$ the energy density of the scalar field is roughly constant^66,67,68,69 and, in average, the expansion rate is $\delta =(q+1)∕(2q-1)$. Thus, ${\delta }_{min}$ is still $\mathcal{𝒪}(1∕2)$ and this happens when $q\gg 1$. When all the ${\delta }_i$ are equal the product of all the ${\xi }_i$ (denoted by ${\xi }_r$ in Eq. (6.19)) is ultimately raised to the same common power implying that the contribution of the whole decelerated stage of expansion of Eq. (6.18) is maximized by a single expanding stage characterized by $\overline{\delta }={\delta }_{min}<1$. Thanks to Eqs. (6.19) we therefore obtain the following bound on ${\nu }_{max}$ : where it has been assumed that $H_r\to H_{\mathrm{\text{bbn}}}=10^{-42}{M}_P$. If we now consider together Eqs. (6.20) and (3.51) we must conclude that the quantum bound of Eq. (3.51) is always more constraining.¹⁰⁰

6.3. Interplay between low-frequency and high-frequency constraints

The previous considerations suggest an interplay between the low-frequency constraints and the high-frequency bound. We are going to examine first the bounds on the inflationary potential coming from the high-frequency region and their connection with the low-frequency limits of Sec. 4. In the second part of this discussion, we swiftly describe some notable quantum mechanical aspects of the relic gravitons at high frequencies.

6.3.1. General bounds on the inflationary potential

Let us suppose, as suggested in Sec. 4, that the inflationary potential interpolates between two complementary regimes: it is inflationary for $\mathrm{\Phi }=\phi ∕{\overline{M}}_P\gg 1$ while it oscillates as $V(\mathrm{\Phi })=V_0{\mathrm{\Phi }}^{2q}$ in the limit $\mathrm{\Phi }\ll 1$. Few examples of this class of potentials have been illustrated in section 4 (see, in particular, Eqs. (4.6), (4.7), (4.8) and (4.9)). In this situation there is no absolute bound on the value of q but the parameter space of the model is effectively three-dimensional: $r_T$ controls the low-frequency normalization, $H_r∕M_P$ determines the reheating scale and q fixes the high-frequency spectral index of $h_0^2{\mathrm{\Omega }}_{gw}(k,{\tau }_0)$ according to

From the specific form of the spectral energy density at high frequencies we may require that the BBN constraint is satisfied while, in the audio band we may require This condition roughly guarantees the enforcement of the constraints of Table 1 together with a potentially detectable signal (in the far future); Eq. (6.22) has the same content of Eq. (73) with the difference that we consider here a slightly larger interval in the spectral energy density. For different values of $r_T$ the bounds are modified and this aspect illustrates once more the interplay between the constraints coming from different frequency regimes, as originally suggested in Refs. 107 and 108. The bounds on $r_T$ appear in the aHz region whereas the bounds related to $n_T^{(\mathrm{\text{high}})}$ come from the high-frequency range. In Fig. 18, we consider two different values of $r_T$ and draw the allowed region in the plane $(q,H_r∕M_P)$. The value of $H_r$ must always be larger than $10^{-44}{M}_P$ roughly corresponding to the BBN scale. In Fig. 18, we considered $r_T$ as a free parameter even though its potential suppression occurs via the total number of e-folds which is always larger than 60 as long as $q>1$. This result should compared also with the low-frequency limits on q, $r_T(k)$ and $n_s(k)$ discussed in Sec. 4; see also, in this respect, the analysis of Ref. 70.

6.3.2. Quantum sensing and the relic gravitons

We already established that the quantum bound is more constraining than the classical limit of Eq. (6.20) and this is true in general terms since Eq. (3.51) does not depend on the specific timeline of the post-inflationary evolution but just on the observation that at the maximal frequency only one graviton pair is produced. It makes then sense to normalize the chirp amplitude directly in the THz domain^hh; with this logic the bound on ${\nu }_{max}$ of Eq. (3.51) can be converted into a limit on $h_c$. If the spectral energy is normalized in the THz domain with a putative high-frequency slope ${\nu }^{n_T^{(\mathrm{\text{high}})}}$, the minimal chirp amplitude required for the direct detection of cosmic gravitons must comply with the following limit :

This means that a sensitivity $\mathcal{𝒪}(10^{-20})$ or even $\mathcal{𝒪}(10^{-24})$ in the chirp amplitude for frequencies in the MHz or GHz regions is irrelevant for a direct or indirect detection of high-frequency gravitons. It has been suggested long ago that microwave cavities^{119,120,121,122,123,124} operating in the MHz and GHz regions could be employed for the detection of relic gravitons.^45,46,47 The same class of instruments has been also invoked in Refs. 114, 115, 116, 117 with the difference that, unlike previous studies (more aware of the potential sources and of the instrumental noises), the required chirp amplitudes are now optimistically set in the range $h_c^{(\mathrm{\text{min}})}=\mathcal{𝒪}(10^{-20})$ for arbitrarily high frequencies.ⁱⁱ Equation (6.23) also clarifies why $h_c^{(min)}$ must be at least $\mathcal{𝒪}(10^{-32})$ (or smaller) for a potential detection of cosmic gravitons in the THz domain. In a more optimistic perspective, for $n_T^{(\mathrm{\text{high}})}>2$ the largest signal occurs at the largest frequency, for $n_T^{(\mathrm{\text{high}})}\le 2$ the frequencies smaller than the THz are observationally convenient. If we consider, for instance, the case $n_T^{(\mathrm{\text{high}})}\to 1$ (which is, incidentally, typical of a post-inflationary stiff phase when we neglect here all the possible logarithmic enhancements) we would have that the chirp amplitude at in the MHz range could be $\mathcal{𝒪}(10^{-28})$ (as also proposed in Refs. 123 and 124 on the basis of more experimental considerations). Furthermore, when $n_T^{(\mathrm{\text{high}})}\to 2$ (typical of the ekpyrotic scenario) we would have instead that $h_c(\nu ,{\tau }_0)$ is the same at higher and smaller frequencies.^212,213 Finally for $n_T^{(\mathrm{\text{high}})}\to 3$ (as it happens in the case of the pre-big bang scenario^214,215) the chirp amplitude at lower frequencies gets even smaller. We have therefore a trade-off between the optimal frequency, the features of the signal and the noises (especially the thermal one) indicating that the highest possible frequency (close to ${\nu }_{max}$) is not always the most convenient. Also this aspect should be taken into account if the goal is really an accurate assessment of the required sensitivities of high-frequency instruments.

The limits following from Eq. (3.51) are also relevant for the analysis of the statistical properties of the relic gravitons and, in particular, of their degrees of first- and second-order coherence. These observables follow by generalizing the appropriate Glauber correlators^216,217 to the expectation values of tensor fields (see Ref. 218 and discussions therein); besides the physical aspects (discussed over a decade ago²¹⁹) the main technical difference between the gravitons and the photons involves the polarization structure of the correlation functions. Mutatis mutandis the physical idea is however similar: if cosmic gravitons are detected by independent interferometers the correlated outputs are employed to estimate the degrees of second-order coherence. The analysis of the interplay between the Hanbury Brown–Twiss (HBT) interferometry and the high-frequency gravitons has been recently discussed in literature (see also Refs. 218 and 219); for the present purposes we avoid the polarization dependence and introduce the single-particle (inclusive) density^220,221

where $Â(\overrightarrow{k},\tau )$ (and its Hermitian conjugate) are just a set of creation and annihilation operators that are nonzero inside the volume of the particle source associated with the three-dimensional integral (in real space) of an appropriate window function $\mathcal{𝒲}(\overrightarrow{x})$. By definition, $[Â(\overrightarrow{k}),{Â}^{†}(\overrightarrow{k})]=\int d^{3}x|\mathcal{𝒲}(\overrightarrow{x}){|}^2$. In the theory of Bose–Einstein interference^220,221 is the two-particle inclusive density and according to Eqs. (6.24) and (6.25) the normalized second-order correlation function estimates the degree of second-order coherence.^218,219 The value of $C_2({\overrightarrow{k}}_1,{\overrightarrow{k}}_2)$ is always enhanced in comparison with so-called Poissonian limit so that the statistics of the relic gravitons is always super-Poissonian and generally super-chaotic. Indeed in the limit of a large number of graviton pairs $C_2({\overrightarrow{k}}_1,{\overrightarrow{k}}_2)\to 3$ whereas $C_2({\overrightarrow{k}}_1,{\overrightarrow{k}}_2)\to 2$ in the case of a chaotic mixture. This result is slightly refined by taking into account the polarization structure of the correlators, as already discussed in the past^218,219; in this case $C_2({\overrightarrow{k}}_1,{\overrightarrow{k}}_2)\le 3$ but the statistics always remains super-Poissonian. While the statistical properties of the relic gravitons determine the degrees of first- and second-order coherence, their potential detection depends from the achievable $h_c^{(min)}$ which is not the same in different ranges of comoving frequency. It is therefore not surprising that the analyses of the Bose–Einstein correlations overlooking the physical properties of the cosmic gravitons are inconclusive and often superficial. All in all we demonstrated, both at the classical and quantum level, that the largest frequency of the relic gravitons never exceeds the THz band while the minimal detectable chirp amplitude should be at least $\mathcal{𝒪}(10^{-32})$ (or smaller) if the (hypothetical) detectors in the THz domain could claim (even in principle) the detection of a relic signal. However, if the pivotal frequencies of the instruments are reduced from the THz to the GHz (or even MHz) band the minimal required chirp amplitude may increase.

6.3.3. The quantumness of relic gravitons

The relic gravitons are characterized by autocorrelation functions that are not invariant under a shift of the time coordinate (see Sec. 3 and discussion therein); this is why Eqs. (3.27)–(3.28) do not only depend upon ${\tau }_1-{\tau }_2$ but also upon ${\tau }_1+{\tau }_2$. This property is rooted into the quantum mechanical origin or the corresponding particles: the initial travelling waves associated with the quantum fluctuations turn eventually into a collection of standing waves because of the evolution of the underlying background geometry. The formation of standing waves (also called Sakharov oscillations) simply means that relic gravitons are produced in entangled states of opposite (comoving) three-momenta according to the unitary process summarized by Eqs. (3.41) and (3.42). Although the field is initially in a pure state its entropy may increase if some information is lost and, for this reason, quantum measurements are somehow intrinsically associated with a loss of information. When observations are performed (for instance by means of HBT interferometry^218,219) the sign of the three-momentum cannot be determined; in other words only one of the members of pair is observed while the other one is in practice unobservable. The operators associated with the opposite momenta of a graviton pair effectively act on separated subspaces of the total Hilbert space of the problem. We can then focus on a single pair of gravitons so that the associated operators will be ${\widehat{b}}_{+}$ and ${\widehat{b}}_{-}$ (i.e. the signal and the idler mode in a quantum optical context^92,93). In this two-mode approximation the final state of the particle production process schematically corresponds to

where $[{\widehat{b}}_I,{\widehat{b}}_J^{†}]={\delta }_{I,J}$; here $I,J=+,-$ and the $\pm$ are related to the sign of a (single) comoving three-momentum. The operator $\mathrm{\Sigma }(z)$ can be factorized as the product of the exponentials of $L_0$ and $L_{\pm }$ where $L_0$ and $L_{\pm }$ are the generators of the $SU(1,1)$ Lie algebra : obeying the corresponding commutation relations $[L_{+},L_{-}]=-2L_0$ and $[L_0,L_{\pm }]=\pm L_{\pm }$. The operator ${\widehat{b}}_{+}^{†}$ creates a graviton of momentum $+\overrightarrow{k}$ while ${\widehat{b}}_{-}^{†}$ creates a graviton with momentum $-\overrightarrow{k}$; the Fock states are an appropriate basis for the irreducible representations^jj of the $SU(1,1)$ When the relic gravitons are produced in pairs of opposite three-momenta we have that $n_{+}=n_{-}$; furthermore the action of the group generators on the two-mode vacuum is given by $L_{-}|0_{+}{0}_{-}〉=0$ while for $L_0$ we have instead $L_0|0_{+}{0}_{-}〉=|0_{+}{0}_{-}〉∕2$. The density matrix associated with the state given in Eq. (6.28) is $\widehat{\rho }=|z〉〈z|$ and since the states are correctly normalized we can conclude that $\mathrm{\text{Tr}}{\widehat{\rho }}^2=\mathrm{\text{Tr}}\widehat{\rho }=1$. In the Fock basis of Eq. (6.30) the explicit form of $|z〉$ is Equation (6.31) seems unnecessarily complicated: on the one hand, we summed over $n_{+}$ and $n_{-}$ while, on the other hand, we included the ${\delta }_{n_{+}{n}_{-}}$ that effectively cancels one of the two sums and enforces the conservation of the three-momentum. The redundant form of Eq. (6.31) is however convenient in what follows since the action of the operators acting over the different subspaces of the total Hilbert space is immediately clear. Indeed, the signal and the idler modes (i.e. ${\widehat{b}}_{+}$ and ${\widehat{b}}_{-}$, respectively) act on two different Hilbert subspaces; they may actually arise as ingredients of a quantum measurement but they correspond here to gravitons with opposite three-momenta. We can always construct a set of Hermitian observables acting on one of the two Hilbert subspaces; for instance ${\widehat{N}}_{+}={\widehat{b}}_{+}^{†}{\widehat{b}}_{+}$ is the averaged multiplicity of the signal whereas ${Î}_{+}={\widehat{b}}_{+}^{†}{\widehat{b}}_{+}^{†}{\widehat{b}}_{+}{\widehat{b}}_{+}$ measures the intensity of the signal. Similar operators can be introduced for the idler mode by replacing $+\to -$. If we eventually average these operators and take their ratio we obtain, always in the case of the signal, $g_{+}^{(2)}=〈{Î}_{+}〉∕{〈{\widehat{N}}_{+}〉}^2$; $g_{+}^{(2)}$ which is the degree of second-order coherence appearing in the analysis of the HBT correlations^92,93 (see also Refs. 222 and 223). Let us now pretend to measure $〈z|{\widehat{N}}_{+}|z〉$; we have, from Eq. (6.31), that The result of Eq. (6.32) can also be expressed in a more transparent form way by introducing the averaged multiplicity $\overline{n}={sinh}^{2}r$ where $p_n$ now denotes the Bose–Einstein (geometric) distribution. The same discussion of Eq. (6.32) can be generalized to ${Î}_{+}$ implying that the analog of Eq. (6.32) becomes

6.3.4. The entanglement entropy

Equations (6.32) and (6.34) ultimately suggest that from the total density matrix $\widehat{\rho }$ a reduced density matrix can be obtained by tracing over the idler mode. To simplify the phases we can introduce

so that the total density matrix in the Fock basis reads Equation (6.36) still describes a pure quantum state but if we now trace over the idler oscillator we obtain a reduced density operator that only depends on the signal : where by definition ${\mathrm{\text{Tr}}}_{-}[\cdots ]$ denotes the trace over the idler mode. The final result for the reduced density operator becomes therefore All the Krönecker deltas appearing in Eq. (6.37) can be used by recalling that $〈k_{-}|m_{-}{m}_{+}〉={\delta }_{k_{-}{m}_{-}}|m_{+}〉$ and that, similarly, $〈n_{+}{n}_{-}|k_{-}〉={\delta }_{n_{-}{k}_{-}}〈n_{+}|$. Note also that in Eqs. (6.37) and (6.38) the summation index (i.e. $n_{+}\to n$) has been renamed. As in Eq. (6.32) the statistical weights of Eq. (6.38) are $p_n={\overline{n}}^n∕{(\overline{n}+1)}^{n+1}$ and they correspond to the Bose–Einstein probability distribution even if the averaged multiplicity $\overline{n}={sinh}^{2}r$ is nonthermal. The reduced density matrix can be used to compute all the correlation functions relevant for the description of the signal but it carries no information of the idler mode. The loss of information associated with the trace over the idler mode is measured by the von Neumann entropy computed from Eq. (6.31) : When a portion of the system is unobservable (or when observations are confined to a subset of the degrees of freedom) information is lost and the total density matrix can then be reduced. Equation (6.39) quantitatively describes the loss of information associated with the trace over the idler mode. Furthermore, from the result of Eq. (6.39) we see that in the limit $\overline{n}\gg 1$ the entropy gets proportional to $ln\overline{n}$, in other words We may now recall that in the process of particle production described by Eqs. (3.41) and (3.42) the averaged multiplicity for each k mode is proportional to $|{\beta }_k(\tau ){|}^2$ which we can estimate from Eq. (3.85); this means that the result of Eq. (6.40) can also be expressed as The result of Eqs. (6.40) and (6.41) hold, strictly speaking, for two oscillators with opposite three-momenta; the previous results are however valid for each pair of $\overrightarrow{k}$-modes of the field so that the density matrix of Eq. (6.31) can also be written as where, for completeness, we have also considered the contribution of a further phase ${\alpha }_k$ that is different for each k-mode. The density matrix can be written, in the Fock basis, as The multimode probability distribution appearing in Eq. (6.43) is given by where ${\overline{n}}_k$ is the average multiplicity of each Fourier mode. Furthermore, following the standard notation, $|\{n\}〉=|n_{{\overrightarrow{k}}_1}〉|n_{{\overrightarrow{k}}_2}〉|n_{{\overrightarrow{k}}_3}〉\dots$ where the ellipses stand for all the occupied modes of the field. We also note that the density matrix can be reduced by considering the phases of the final multiparticle state to be unobservable. In this case the right-hand side of Eq. (6.42) can be averaged averaging over ${\alpha }_k$. The reduced density matrix would be given, in this case, by This observation represents a further reduction scheme of the density matrix that ultimately leads to the same entropy of Eq. (6.39) in the limit of the large averaged multiplicities. To obtain the total entropy we must integrate over the whole spectrum for a fiducial Hubble volume at the present time; in this case the total entanglement entropy of the gravitons becomes The integral appearing in Eq. (6.46) can be estimated by using either the general form of the averaged multiplicity deduced in in Eq. (3.50) or the result of Eq. (3.85). For The power-law parametrization of Eq. (3.50) is actually compatible with Eqs. (3.85) where the averaged multiplicity for $k\ll k_{max}$ depends on $(a_{\mathrm{\text{re}}}∕a_{\mathrm{\text{ex}}})\gg 1$. If assume that between $a_{\mathrm{\text{ex}}}$ and $a_{\mathrm{\text{re}}}$ the background expands in this case the averaged multiplicity is given by^kk where the subleading contributions have been neglected since they do not affect the final value of the integral. Thanks to the explicit form of Eq. (6.31) the integration variable appearing in Eq. (6.46) can be rescaled and the total entropy of the gravitons is where $\mathcal{𝒦}(m_T)(4-m_T){\int }_{{\nu }_{min}∕{\nu }_{max}}^1{x}^{2}lnxdx=\mathcal{𝒪}(1)$ is a numerical factor that depends on the spectral slope $m_T=n_T^{(\mathrm{\text{high}})}$. Since in the integral $\mathcal{𝒦}(m_T)$ the value of $m_T$ is not essential and the lower limit of integration goes to zero (at least in practice since ${\nu }_p\simeq {\nu }_{min}=310^{-18}\mathrm{\text{Hz}}$) we have, from Eq. (6.48), that $S_{g0}=\mathcal{𝒪}(10){({\nu }_{max}∕H_0)}^3$. The result for $S_{g0}$ can now be compared with the well-known result of the thermal entropy of the CMB computed within the same fiducial Hubble volume Since $T_{\gamma 0}=T_{\gamma 0}=(2.72548\pm 0.00057)\mathrm{\text{K}}$ we also have that If we now require that $S_{g0}\le S_{\gamma 0}=\mathcal{𝒪}(10^{90})$ we have from Eqs. (6.48), (6.49) and (6.50) that The maximal frequency of the spectrum deduced in Eq. (3.51) implies that the entanglement entropy of the gravitons cannot exceed the total entropy of the CMB. The quantum theory of parametric amplification however provides a natural cosmological arrow associated with an entanglement entropy.^{225,226,227,228,229,133} In our context the growth of the averaged multiplicity of the produced gravitons is naturally associated with the increase of the entanglement entropy of the gravitational field. In an idealized experiment based, for instance, on the analysis of the gravitational analog of the HBT correlations only one of the two gravitons of the pair is typically detected. In this situation the resulting entanglement entropy is proportional to the logarithm of the averaged multiplicity. Since the final multiplicity is always large, the entropy is effectively proportional to the squeezing parameter r.^133,229 The reduction of the density matrix can be performed in different bases^{133,227,228,229}; the results depend on the basis but the asymptotic limit when the averaged multiplicity is large is generally proportional to $2r$ and hence to the logarithm of $\overline{n}$. The coarse grained entropy employed here is quite close to the notion of information theoretic entropy introduced many years ago^230,231 with the purpose of reformulating the indetermination relations. This idea actually inspired the reduction scheme discussed in Refs. 227 and 228. We finally stress that the value of the entropy of the gravitons for each single mode is equal to $2r$ and this result follows by considering the entropy as a measure of uncertainty in the result of a measurement or preparation of a given quantum mechanical observable. Let us consider the following canonical transformation of a two-mode harmonic oscillator : with $r>0$. Since the transformation is canonical the conjugate momenta will transform $p\to p_r=e^{r}p$ and as $\tilde{p}\to {\tilde{p}}_r=e^r\tilde{p}$. It is clear from Eq. (6.52) that operators x and $\tilde{x}$ fluctuate above the quantum noise since for a quantum state $|\psi 〉$ the corresponding normalized wave function reads From Eq. (6.53) it is simple to compute ${(\mathrm{\Delta }x)}^2=〈x^2〉-{〈x〉}^2$ (and similarly for $(\mathrm{\Delta }{\widetilde{x)}}^2$); thanks to Eq. (6.53) $\mathrm{\Delta }x=\sqrt{{(\mathrm{\Delta }x)}^2}=e^r∕\sqrt{2}$ and $\mathrm{\Delta }\tilde{x}=\sqrt{{(\mathrm{\Delta }\tilde{x})}^2}=e^r∕\sqrt{2}$. Both operators x and $\tilde{x}$ fluctuate above the quantum noise and a natural measure of uncertainty in the result of a measurement of the superfluctuant operators is given by Inserting Eq. (6.53) into Eq. (6.54) we obtain $s_{\mathrm{\text{super}}}=2r+ln(e\pi )$ which coincides with $2r$ in the limit of large produced particles. The coarse grained entropy analyzed here is then quite close to the notion of information theoretic entropy introduced many years ago^230,231 with the purpose of reformulating the indetermination relations. This idea actually inspired the reduction scheme discussed in Refs. 227 and 228.

7. Concluding Remarks

Since the Universe became transparent to the propagation of electromagnetic disturbances only after matter–radiation equality, the photons coming from the primeval stages of the evolution of the plasma cannot be detected so that earlier tests on the expansion history are actually related to the remarkable successes of BBN taking place when the expansion rate was of the order of $10^{-44}{M}_P$. This figure should be compared with the approximate value of the inflationary expansion rate (i.e. $\mathcal{𝒪}(10^{-6})M_P$) inferred from the amplitudes of the curvature inhomogeneities that affect the CMB temperature and polarization anisotropies. Between these two scales the expansion rate spanned 38 orders of magnitude where the evolution of the plasma could have been rather different from radiation.

During the last 50 years the interplay between high-energy physics and cosmology has been guided by the assumption that radiation should be the dominant component of the plasma well before the onset of BBN and immediately after the end of inflation. This conventional wisdom is consistent both with an early stage of inflationary expansion and with the concordance scenario at late times but it is not unique. The physical foundations of this paradigm are not corroborated by direct observations and they could be either partially or totally refuted in the years to come. Since during inflation the particle horizon diverges (while the event horizon is finite) all the wavelengths that are currently shorter than the Hubble radius were in causal contact during inflation provided the overall duration of inflation was sufficiently long. The length of the inflationary stage is customarily assessed in terms of the number of e-folds which should be $\mathcal{𝒪}(60)$ if the post-inflationary expansion rate is dominated by radiation. This estimate can be either reduced (down to $\mathcal{𝒪}(45)$) or increased (up to $\mathcal{𝒪}(75)$) depending on the post-inflationary expansion rate that may become either faster or slower than radiation, respectively.

Any presumption about the timeline of the expansion rate should necessarily acknowledge that every variation of the space–time curvature produces shots of gravitons with specific averaged multiplicities. After the actual detection of gravitational radiation there are no direct physical limitations forbidding the empirical scrutiny of the spectra of the relic gravitons (either in the audio band or in higher frequency domains) within the following score year. Since different timelines ultimately correspond to specific profiles of $h_0^2{\mathrm{\Omega }}_{gw}(\nu ,{\tau }_0)$ (for frequencies ranging between the aHz and the THz), the expansion rate can be systematically inferred from the slopes of the observed spectra and from their pivotal frequencies. The results outlined here specifically address the interplay between the expansion history of the plasma and the spectral energy density of the relic gravitons in the concrete situations inspired by the current phenomenological lore at low, intermediate and high frequencies.

•	The inflationary observables in the aHz region depend on the timeline of the post-inflationary evolution. In single-filed inflationary scenarios this means, in particular, that the tensor to scalar ratio and the scalar spectral index are more or less suppressed if the timeline of the expansion rate is either slower or faster than radiation respectively. At higher frequencies the PTAs (operating in the nHz range) are now setting interesting bounds on the post-inflationary expansion rate. The apparent excesses appearing in the last data releases of two PTAs could actually come from an increasing spectrum of relic gravitons at intermediate frequencies.
•	Between the $\mu \mathrm{\text{Hz}}$ and the Hz various space-borne detectors might be operational in the far future although the signals expected in the mHz region are dominated by astrophysical sources (e.g. galactic white dwarves, solar-mass black holes, supermassive black holes coming from galaxy mergers). The only cosmological sources customarily considered in this framework are associated with the phase transitions at the TeV scale although perturbative and nonperturbative estimates consistently suggest that the standard electroweak theory leads to a cross-over regime where drastic deviations from homogeneity (and the consequent bursts of gravitational radiation) should not be expected. The inflationary signal (often regarded as irrelevant between few $\mu$Hz and the Hz) could be in fact much larger that the purported signal coming from a realistic dynamics at the electroweak scale. Moreover, since the slopes of $h_0^2{\mathrm{\Omega }}_{gw}(\nu ,{\tau }_0)$ obtained in the case of a putative strongly first-order phase transition are much steeper than the ones associated with a modified expansion history, the most severe phenomenological bounds on the relic gravitons between the $\mu \mathrm{\text{Hz}}$ and the Hz arise (by continuity in frequency) from the audio band and from the operating ground-based detectors.
•	The window of wide-band detectors notoriously ranges between few Hz and 10kHz. The current limits imply that the sensitivity of correlated interferometers for the detection of a flat spectral energy density of relic gravitons is approximately $h_c^{(min)}=\mathcal{𝒪}(10^{-24})$ for typical frequencies in the audio band. Sharp deviations from scale-invariance lead to similar orders of magnitude and while these figures may improve in the years to come, the frequency domain of ground-based interferometers will remain the same. For this reason it is important to promote new instruments operating in higher frequency domains where the potential signals coming from the past history of the plasma are dominant. More than twenty years ago it was suggested that microwave cavities (operating between the MHz and the GHz regions) could be used for the detection of relic gravitons associated with post-inflationary phases stiffer than radiation. While forty years ago the typical sensitivities of these instruments were $h_c^{(min)}=\mathcal{𝒪}(10^{-17})$ they improved later on and reached $h_c^{(min)}=\mathcal{𝒪}(10^{-20})$ in the early 2000s. Similar prototypes aimed at the detection of dark matter could be used as high-frequency detectors of gravitational waves. The target sensitivities of these instruments are often set by requiring in the MHz (or even GHz regions) the same sensitivities reached (today) by the interferometers in the audio band. This means that the features of the instruments are not guided by the signals of the available sources in the corresponding frequency domain. To detect directly relic gravitons with high-frequency instruments operating between the MHz and the GHz the minimal detectable chirp amplitude should be $h_c^{(\mathrm{\text{min}})}=\mathcal{𝒪}(10^{-32})$ (or smaller). However, if the pivotal frequencies of the instruments are reduced from the THz to the GHz (or even MHz) band the minimal required chirp amplitude may increase. With these specifications, the detectors in the MHz and GHz domains may be able to probe directly the relic gravitons and their quantumness.

Both at the classical and quantum level, the largest frequency of the relic gravitons never exceeds the THz band and above the maximal frequency the averaged multiplicity is exponentially suppressed so that ${\nu }_{max}$ ultimately corresponds to the production of a single graviton pair. Since the relic gravitons are inherently quantum mechanical, their quantumness can be measured in terms of an entanglement entropy that is caused by the loss of the complete information on the underlying quantum field. The reduction of the density matrix in different bases leads to the same von Neumann entropy whose integral over all the modes of the spectrum is dominated again by the maximal frequency. Whenever the THz bound is applied, it turns out that the total integrated entropy of the relic gravitons is comparable with the entropy of the CMB but not larger. A potential detection of relic gravitons both at low and high frequencies may therefore represent a direct evidence of macroscopic quantum states associated with the gravitational field. For this reason the detectors operating in the MHz and GHz regions are quantum sensitive to the second-order interference effects. As in the case of optical photons, the interferometric techniques pioneered by Hanbury Brown and Twiss in the 1950s could be applied to high-frequency gravitons with the purpose of distinguishing the statistical properties of thermal and nonthermal gravitons.

Acknowledgments

I wish to acknowledge relevant discussions with the late Ph. Bernard, G. Cocconi and E. Picasso on high-frequency gravitons and microwave cavities. It is a pleasure to thank A. Gentil-Beccot, P. Birtwistle, A. Kohls, L. Pieper, S. Rohr and J. Vigen of the CERN Scientific Information Service for their kind help along the different stages of this investigation. Some of the discussions presented here have been developed on the occasion of few seminars and of a set of lectures; I thank the questions and the remarks of students and colleagues.

Notes

^a Indeed, if we reach the Planck time, the blueshifted value of the Hubble radius is of the order of $4\mu \mathrm{\text{m}}=4\times 10^{-4}\mathrm{\text{cm}}$. But since at the Planck time the particle horizon is $d_p(t_P)\simeq t_P\simeq 10^{-33}$cm, the ratio between $4\times 10^{-4}\mathrm{\text{cm}}$ and $10^{-33}\mathrm{\text{cm}}$ is approximately $\mathcal{𝒪}(10^{29})$ and the number of disconnected volumes is $\mathcal{𝒪}(10^{87})$.

^b Although this point is often ignored we like to point out that the limit $t_{min}\to -\infty$ is not well defined; strictly speaking an ever expanding inflationary evolution is not past geodesically complete.⁴⁹ The limit $t_{min}\to -\infty$ can be better defined by introducing a geodesically complete extension of the de Sitter space–time. This problem has been discussed in the past but will not be specifically addressed here.

^c Throughout this paper, $ln$ denotes the natural (Naperian) logarithm, $log$ indicates instead the common logarithm.

^d The inhomogeneity of the total pressure $\delta p_t$ can be decomposed as $\delta p_t=c_{st}^2\delta {\rho }_t+\delta p_{\mathrm{\text{nad}}}$ where $c_{st}$ is the sound speed of the plasma.

^e We privilege the approximate expressions for the evolution of the mode functions but the final results coincide with more accurate (and conventional) strategies such as the ones based on the exact evolution of the mode functions during the inflationary stage (see, in this respect, the discussion of Appendix A).

^f In general we have $r_T(k,\tau )$ but when $k=\mathcal{𝒪}(k_p)$ and $\tau \simeq 1∕k$ we obtain the standard value of $r_T$ which is customarily quoted in the literature.^42,43,44

^g Throughout this paper, the scale factor is normalized as $a_0=1$. This remark is quite relevant since by choosing $a_0=1$ we will have that comoving and physical frequencies of the relic gravitons coincide at the present time.

^h The difference due to $g_s$ and $g_{\rho }$ in the final results is actually negligible for the present purposes and it involves a factor $1.3$ (instead of 1) at the level of Eq. (2.42). However, from the conceptual viewpoint this difference is certainly relevant and this is why it will be taken into account.

ⁱ The result of Eq. (2.49) is in fact obtained from Eq. (2.48) by recalling that $H_k∕M_P=\sqrt{\pi {ϵ}_k{\mathcal{𝒜}}_{\mathcal{ℛ}}}$. From the consistency relations we also have that $r_T\simeq 16{ϵ}_k$ so that, for $r_T=0.06$ Eq. (2.49) demands that the value of ${\overline{N}}_k$ is given by ${\overline{N}}_k=59.7384$ (while all the other parameters are kept fixed at their typical values).

^j To avoid confusions $w_{min}$ and $w_{max}$ indicate, respectively, the minimal and the maximal values of w.

^k We stress, in this respect, that the indetermination on $N_{max}$ is not related to the considerations discussed in Eq. (2.44): in that context N denoted the total number ofe-folds which may be, for different reasons, larger than $N_{max}$.

^l For the moment it is sufficient to note that, for monomial inflationary potentials, the tensor to scalar ratio scales as $N_k^{-1}$ whereas for plateau-like potential the same quantity scales as $N_k^{-2}$. Both values may get eventually larger or smaller than in the radiation phase depending on the post-inflationary expansion rate. Moreover, as we shall see, the value of $N_k$ ultimately affects the value of the maximal frequency of the relic graviton spectrum.

^m It can also happen that the background evolution at early times is genuinely different from a stage of inflationary expansion. Both possibilities will be swiftly mentioned later on in Sec. 5.

ⁿ Since the tensor amplitude $h_{ij}(\overrightarrow{x},\tau )$ is real the corresponding Fourier amplitude must obey $h_{ij}^{\ast }(\overrightarrow{k},\tau )=h_{ij}(-\overrightarrow{k},\tau )$. Moreover $h_{ij}(\overrightarrow{x},\tau )$ is also solenoidal and traceless; thus $h_{ij}^{\ast }(\overrightarrow{k},\tau )$ must obey ${\widehat{k}}^i{h}_{ij}(\overrightarrow{k},\tau )={\widehat{k}}^j{h}_{ij}(\overrightarrow{k},\tau )=0$ and $h_i^i=0$. See also the considerations developed in Appendix B.

^o Sometimes in the literature Eq. (3.10) is taken as definition of ${\mathrm{\Omega }}_{gw}(k,\tau )$. This is also incorrect since Eq. (3.10) is only an approximation that holds for wavelengths that are sufficiently small in comparison with the effective horizon (or, in equivalent terms, wavenumbers much larger than the expansion rate).

^p In particular the intrinsic noises of the instruments are customarily assumed to be stationary, Gaussian, uncorrelated, much larger in amplitude than the gravitational strain, and statistically independent on the strain itself. The stationarity and the homogeneity are also conjectured for the signals associated with the diffuse background of gravitational radiation.⁸⁸ So far we demonstrated that the diffuse backgrounds of relic gravitons are homogeneous in space but to address the stationarity it is instead essential to take into account the quantum mechanical aspects of the problem.

^q These standing oscillations are in fact related to the tensor analog of the Sakharov oscillations^71,72,73 (see also Ref. 49). Both during the radiation phase and in the matter epoch the standing oscillations appearing in the power spectrum lead to nonstationary features in the autocorrelation function.⁷⁶

^r While it is debatable if the nonstationary features associated with the diffuse backgrounds of relic gravitons are (or will be) directly detectable, the spectral amplitude following from the Wiener–Khintchine theorem is generally inappropriate for a consistent description of the relic signal.

^s The rationale for the bound of Eq. (3.47) is discussed in Sec. 5.

^t The analyses of the Bose–Einstein correlations, however, cannot be pursued in spite of the properties of the sources; this is why to overlook the physical properties of the cosmic gravitons leads to conclusions that are ambiguous and ultimately superficial. It makes actually little sense to consider potentials signals coming from diffuse backgrounds for arbitrarily large frequencies (possibly much larger than the THz) without bothering about the underlying physical constraints. This approach is probably motivated by the need of claiming large sensitivities for potential instruments but has no physical basis unless the class of bounds related to Eq. (3.51) is understood and acknowledged. We shall get back to the quantumness of the relic gravitons at the end of Sec. 6.

^u This also means that the maximal frequency, the intermediate frequencies and the shape of ${\mathrm{\Omega }}_{gw}(\nu ,{\tau }_0)$ can all be employed, in different combinations, to infer timeline of the expansion rate as we are going to see in the following sections.

^v In Eq. (3.74), we restored $M_P$ by recalling its relation with ${\overline{M}}_P$, i.e. ${\overline{M}}_P=M_P∕\sqrt{8\pi }$.

^w The result of Eq. (3.75) holds, strictly speaking, in the case $\delta \ge 1∕2$. When $\delta <1∕2$ the contribution of ${\mathcal{𝒥}}^{(t)}({\tau }_{\mathrm{\text{ex}}},{\tau }_{\mathrm{\text{re}}})$ must be carefully evaluated and contributes to the slope which becomes ${\overline{n}}_T=2\delta +\mathcal{𝒪}(r_T)$.

^x The analysis leading to the results discussed above can be generalized to the situation where there are many post-inflationary stages characterized by different rates ${\delta }_i$. It is also possible to use different approximation schemes that will not be specifically discussed here (see however Ref. 49).

^y This conclusion would follow by appreciating that ${\alpha }_{k,\alpha }(\tau )=[k{f}_{k,\alpha }(\tau )+i{g}_{k,\alpha }(\tau )]∕\sqrt{2k}$ and also that ${\beta }_{k,\alpha }^{\ast }(\tau )=-[k{f}_{k,\alpha }(\tau )-i{g}_{k,\alpha }(\tau )]∕\sqrt{2k}$.

^z For a generic quantity that is both scale-dependent and time-dependent (be it for instance $W(k,\tau )$) we the have that its value is given by $W_k=W(k,\tau )=W(k,1∕k)$ when the CMB wavelengths cross the Hubble radius during inflation.

^aa For instance the BICEP2 observations¹³⁷ suggested $r_T=\mathcal{𝒪}(0.2)$ that looked rather large for single-field inflationary models with monomial potentials. If $N_k\ll 60$ the value of $r_T(N_k)$ is comparatively larger than in the case $N_k={\overline{N}}_k=\mathcal{𝒪}(60)$.

^bb We stress once more that $r_T(k,{\tau }_{\mathrm{\text{ex}}})=r_T(k,1∕k)=r_T(k)$ and similarly $n_s(k,{\tau }_{\mathrm{\text{ex}}})=n_s(k,1∕k)=n_s(k)$.

^cc This happens, for instance, in the single-field case where, thanks to the consistency relations, the tensor spectral index $n_T^{(\mathrm{\text{low}})}$ is related to the tensor to scalar ratio $r_T$ as $n_T^{(\mathrm{\text{low}})}\simeq -r_T∕8$. Since $r_T$ is currently assessed from the analysis of the temperature and polarization anisotropies of the CMB^42,43,44 $n_T^{(\mathrm{\text{low}})}$ cannot be positive.

^dd Unless the relic gravitons would lead exactly to the same slope of the astrophysical foregrounds associated with black-hole binary systems, the value $\beta =-2∕3$ is not particularly compelling in a cosmological setting.

^ee After Ref. 190 appeared in the form of a preprint, some authors made exactly the same speculation and talked about the sound speed (or sound velocity) of the relic gravitons. While this terminology makes little sense in the context of the propagation of massless particles, the idea is exactly the same (see Refs. 193, 194 and references therein). In the present context we prefer to discuss this class of phenomena in terms of an effective refractive index, as originally suggested in Refs. 190–192.

^ff The low-frequency transfer function ${\mathcal{𝒯}}_{\mathrm{\text{low}}}(\nu ∕{\nu }_{\mathrm{\text{eq}}})$ has a definite form^107,108 the high-frequency transfer function ${\mathcal{𝒯}}_{\mathrm{\text{high}}}(\nu ∕{\nu }_r,\delta )$ depends on the value of $\delta$ so that it does not have a general expression.⁶⁵ It should be stressed that we refer here to the transfer function of the spectral energy density^65,107,108 which is numerically more accurate (when estimating ${\mathrm{\Omega }}_{gw}(k,{\tau }_0)$) than the transfer function for the amplitude.^{207,208,209,210,211}

^gg Alternatively we may suppose that the relic gravitons backgrounds will not be accessible in the audio band; in what follows we shall entertain a less pessimistic attitude which is mainly motivated by the steady increase of the sensitivity to relic gravitons in the last 20 years. We must actually recall that in 2004 wide-band detectors gave limits implying $h_0^2{\mathrm{\Omega }}_{gw}(\nu ,{\tau }_0)<\mathcal{𝒪}(1)$³³ while today the same limits improved by roughly 10 orders of magnitude.^34,35

^hh The spectral energy density in critical units at the present time and the chirp amplitude $h_c(\nu ,{\tau }_0)$ are related as ${\mathrm{\Omega }}_{gw}(\nu ,{\tau }_0)=2{\pi }^2{\nu }^2{h}_c^2∕(3H_0^2)$.

ⁱⁱ To achieve $h_c^{(\mathrm{\text{min}})}=\mathcal{𝒪}(10^{-20})$ is technologically interesting; from the physical viewpoint this minimal sensitivity is more than 10 orders magnitude larger than the requirements associated with the direct detection of cosmic gravitons.

^jj An equivalent basis for the irreducible representations of $SU(1,1)$ is provided by the vectors $|Q{n}_t〉$ where $Q=n_{+}-n_{-}$ is the total charge and $n_t=n_{+}+n_{-}$ is the total number of charged species. The vectors $|Q{n}_t〉$ are the standard basis of the irreducible representations $T^{+k}$ of $SU(1,1)$ where k is the principal quantum number and m is the magnetic quantum number, i.e. the eigenvalue of $L_0$. The Casimir operator of the $SU(1,1)$ group can be notoriously written as $C=L_0(L_0-1)-L_{+}{L}_{-}$ so that, eventually, $C|km〉=k(k-1)|km〉$. The commuting set of observables is formed in this case by the Casimir operator and by $L_0$; k is usually referred to as the Bargmann parameter.²²⁴ The negative series $T^{-k}$ is symmetric under the exchange $n_{+}\to n_{-}$ while the principal (continuous) series will not play a specific role in the present considerations. In terms of k and m we have that the total charge and the total number of particles are given, respectively, by $Q=2k-1$ and by $n_t=2m-1$. Note finally that the Bargmann parameter²²⁴ should not be confused with the modulus of the comoving three-momentum; this is actually impossible since the basis of the irreducible representations employed here is the one given in Eq. (6.30) and not the Bargmann basis.

^kk The explicit form of Eq. (6.47) follows by assuming that the relevant wavelengths cross the Hubble radius for the first time during inflation (i.e. $k{\tau }_{\mathrm{\text{ex}}}=k∕{\mathcal{\mathscr{H}}}_{\mathrm{\text{ex}}}=\mathcal{𝒪}(1)$) when the scale factor is given, approximately, by $a_{\mathrm{\text{ex}}}={(-{\tau }_1∕{\tau }_{\mathrm{\text{ex}}})}^{{\beta }_{\mathrm{\text{ex}}}}\simeq |k{\tau }_1{|}^{{\beta }_{\mathrm{\text{ex}}}}$; the reentry takes place instead when $a_{\mathrm{\text{re}}}={({\tau }_{\mathrm{\text{re}}}∕{\tau }_1)}^{{\beta }_{\mathrm{\text{re}}}}\simeq |k{\tau }_1{|}^{-{\beta }_{\mathrm{\text{re}}}}$. In the conventional situation ${\beta }_{\mathrm{\text{ex}}}=1∕(1-ϵ)$ (where $ϵ$ is the slow-roll parameter) and ${\beta }_{\mathrm{\text{ex}}}=\mathcal{𝒪}(1)$: this means that ${\overline{n}}_k\simeq {(k∕k_{max})}^{-4}$. If the wavelengths reenters the Hubble radius during a maximally stiff phase we have instead that ${\overline{n}}_k\simeq {(k∕k_{max})}^{-3}$ since ${\beta }_{\mathrm{\text{re}}}=\mathcal{𝒪}(1∕2)$. The considerations based on Eq. (6.47) are then consistent with Eq. (3.50).

ORCID

Massimo Giovannini  https://orcid.org/0000-0001-6854-2306