Ability of LIGO and LISA to probe the equation of state of the early Universe
AAbility of LIGO and LISA to probe theequation of state of the early Universe
Daniel G. Figueroa and Erwin H. Tanin
Institute of Physics, Laboratory of Particle Physics and Cosmology (LPPC), ´Ecole Poly-technique F´ed´erale de Lausanne (EPFL), CH-1015 Lausanne, Switzerland.E-mail: daniel.figueroa@epfl.ch, [email protected]
Abstract.
The expansion history of the Universe between the end of inflation and the onsetof radiation-domination (RD) is currently unknown. If the equation of state during thisperiod is stiffer than that of radiation, w > /
3, the gravitational wave (GW) backgroundfrom inflation acquires a blue-tilt d log ρ GW d log f = w − / w +1 / > f (cid:29) f RD corre-sponding to modes re-entering the horizon during the stiff-domination (SD), where f RD isthe frequency today of the horizon scale at the SD-to-RD transition. We characterized indetail the transfer function of the GW energy density spectrum, considering both ’instant’and smooth modelings of the SD-to-RD transition. The shape of the spectrum is controlledby w , f RD , and H inf (the Hubble scale of inflation). We determined the parameter spacecompatible with a detection of this signal by LIGO and LISA, including possible changes inthe number of relativistic degrees of freedom, and the presence of a tensor tilt. Consistencywith upper bounds on stochastic GW backgrounds, however, rules out a significant fractionof the observable parameter space. We find that this renders the signal unobservable byAdvanced LIGO, in all cases. The GW background remains detectable by LISA, though onlyin a small island of parameter space, corresponding to scenarios with an equation of state inthe range 0 . (cid:46) w (cid:46) .
56 and a high inflationary scale H inf (cid:38) GeV, but low reheatingtemperature 1 MeV (cid:46) T RD (cid:46)
150 MeV (equivalently, 10 − Hz (cid:46) f RD (cid:46) . · − Hz).Implications for early Universe scenarios resting upon an SD epoch are briefly discussed. a r X i v : . [ a s t r o - ph . C O ] A ug ontents dof during RD 254.3.3 Free-streaming of relativistic particles 28 Observations of the cosmic microwave background (CMB) and the large scale structures of theUniverse strongly support the idea that the Universe underwent a period of cosmic inflationat its early stages. Apart from solving the horizon and flatness problems, and providing theappropriate initial conditions for primordial density perturbations, inflation is a very naturalphenomenon. The most recent measurement of the B-mode polarization anisotropies of theCMB [1, 2] sets a bound on the inflationary Hubble rate which corresponds to an energyscale (cid:46) GeV. At some point below this energy scale, the energy budget of the universemust be converted into a thermal bath of radiation in order to switch to the standard hotBig Bang cosmology, from which point the expansion history is well known. To explain theobserved abundance of light elements in the universe, Big Bang Nucleosynthesis (BBN) musttake place during the radiation domination (RD) epoch, starting at a temperature of around T BBN ∼ (cid:46) GeV and the onset of RD occurring at leastabove ∼ ∼
19 orders of magnitude in energy scale.During inflation, tensor metric perturbations generated from quantum fluctuations arespatially stretched to scales exponentially larger than the inflationary Hubble radius. Thisprocess results in a (quasi-)scale invariant tensor power spectrum at superhorizon scales [3–6], with a very small red tilt. After inflation, the horizon begins to grow faster than theredshifting of length scales, and the tensor modes re-enter the horizon successively as theHubble radius catches up with each wavelength. Once a tensor mode crosses inside thehorizon, it becomes part of the stochastic background of gravitational waves (GWs). Different– 1 –ensor modes re-enter the horizon at different times, and hence propagate through differentperiods of the evolution of the Universe after they have become sub-horizon. For tensormodes crossing during RD, the resulting present GW energy density spectrum is (quasi-)scale invariant, with the same tilt as the original super-horizon inflationary tensor powerspectrum. For periods of expansion when the energy budget of the Universe is not dominatedby relativistic species, the (quasi-)scale invariance is broken, and the resulting GW spectrumbecomes significantly tilted in the frequency range corresponding to the modes crossing thehorizon during such period(s). The expansion history of the Universe is imprinted in thisway in the power spectrum of the freely-lingering primordial GW background that we canobserve today. In other words, we might detect or constrain the post-inflation expansionhistory by attempting to detect the relic GWs of inflationary origin, see e.g. [7–14]. In turn,the expansion history will inform us about the matter fields driving the expansion.To characterize the post-inflation pre-BBN expansion history, we consider an interme-diate epoch taking place between the end of inflation and the onset of RD, with an EoSparameter w (cid:54) = 1 /
3. In scenarios where the inflaton oscillates around the minimum of its po-tential after inflation, an effective EoS (averaged over inflaton oscillations) emerges [15, 16],lying within the range 0 < w < / V ( φ ) ∝ φ p , p <
4. There is however a priori no rea-son, neither theoretical nor observational, to exclude a stiff case / < w ≤ . A period witha stiff EoS can actually be achieved naturally in a universe dominated by the kinetic energyof a scalar field after inflation, either through inflaton oscillations [15] under a steep potential(e.g. V ( φ ) ∝ φ p with p > w (cid:39)
1. We note that w ≤ c = ∂p/∂ρ ≡ w ≤ w can approach unity only frombelow, so we will consider only w <
1. Scenarios with a stiff dominated (SD) epoch betweeninflation and RD are particularly appealing from observational point of view: the stiff periodinduces a blue-tilt in the GW energy spectrum at large frequencies [7, 8, 11, 17–24], openingup the possibility of detection of this GW background by upcoming GW direct-detection ex-periments. Furthermore, the possibility of a stiff epoch is also well motivated on the theoryside, by various early Universe scenarios for which the implementation of a post-inflationarystiff period is crucial. For instance, in
Quintessential inflation [25–33] the inflationary epochis followed by a period where the universe is dominated by the kinetic energy of the inflaton,with the potential adjusted to describe the observed dark energy as a quintessence field. Inthe original gravitational reheating formulation [34, 35] the Universe is reheated by the de-cay products of inflationary spectator fields if the Universe undergoes a sufficiently long stiffepoch after inflation. Even though basic implementations of such gravitational reheatingscenarios have been shown to be inconsistent with BBN/CMB constraints [36], unnatural ad hoc constructions are still viable. Furthermore, variants that can naturally avoid theinconsistency have been also proposed, like the
Higgs-reheating scenario [37], where the Stan-dard Model (SM) Higgs is a spectator field with a non-minimal coupling to gravity and theUniverse is reheated into SM relativistic species (decay products of the Higgs), if inflationis followed by an SD period. The same mechanism can actually be realized with genericself-interacting scalar fields [38], other than the SM Higgs. See [39] for a recent re-analysisof the idea.Models with blue-tilted inflationary GW spectrum due to the presence of an SD epochhave been studied in various contexts [7–9, 11, 18–20, 40, 41]. The resulting GW energy– 2 –pectrum can be characterized by three quantities: the Hubble rate H inf during inflation,the equation state parameter w during the stiff epoch, and the Hubble/energy scale at theSD-to-RD transition, parametrized by the redshifted frequency today f RD correspondingto such scale. In general, a blue-tilted GW background can be probed with a variety ofexperiments, see e.g. [42–52]. The aim of our present work is to assess the ability of the Advanced Laser Interferometric Gravitational wave Observatory (aLIGO), as well as of the(will-be) first generation space-based GW detector, the
Laser Interferometric Space Antenna (LISA), to probe the parameter space spanned by H inf , w , and f RD . We assume inflation iswell described by a quasi- de Sitter stage, so that the inflationary tensor tilt is expected tobe only slightly-red tilted, according to the standard consistency condition [c.f. Eq. (2.7)].A study in this spirit was initiated in [9, 11], but with the SD-to-RD energy scale fixed toits smallest possible value, T BBN ∼ − GeV, so that f RD was fixed to its lowest possiblevalue f BBN ∼ − Hz. In our present work we present a systematic exploration of theobservability of the full parameter space { w, H inf , f RD } .As a consistent cosmological history must preserve the success of BBN, we need toprevent the presence of a stiff epoch from changing significantly the expansion rate duringBBN. There are two ways by which the latter can happen. First, if the stiff epoch does notend in time before the start of BBN the expansion rate would certainly deviate significantlyfrom that of the ΛCDM. Second, even if RD starts in time, GWs should not carry too muchenergy, or otherwise the expansion rate during BBN (or CMB decoupling for this matter)would still be sufficiently altered. Thus, in order to be consistent with BBN/CMB, boundsmust be put on the parameters { w, H inf , f RD } , as these control both the duration of the stiffperiod and the shape of the GW spectrum (and hence the energy carried by the GWs). Aswe will see, these constraints will restrict severely the ability of detectors to observe the bluetilted GW background due to an SD epoch. We find that the above constraints render thesignal completely inaccessible to the observational window of aLIGO, independently of theparameter space. Whilst the background remains detectable by LISA, it is only observable ifthe inflationary scale is as large as H inf (cid:38) GeV (corresponding to at least O (10)% of itscurrent upper bound), f RD lies in the range 10 − Hz (cid:46) f RD (cid:46) . · − Hz, or equivalentlythe Universe becomes RD at sufficiently low temperatures 1 MeV (cid:46) T RD (cid:46)
150 MeV (i.e. theSD spoch spans many decades in energy scale), and the stiff EoS is confined within the narrowrange 0 . (cid:46) w (cid:46) .
56. This corresponds to a small island in the parameter space .The paper is divided as follows. In Section 2, we briefly review the form of the spectrumof GWs from inflation. In Section 3, we derive the energy density spectrum of the inflationaryGW background in the presence of an SD epoch, considering two possible modelings for theSD-to-RD transition: an ’instantaneous’ transition with a sharp jump in the EoS, and a’smooth’ transition modeled by the evolution of two fluid components, one made of radiationand another by a scalar field dominated by its own kinetic energy. In Section 4, we firstquantify the full parameter space { w, H inf , f RD } that can be probed by LISA and aLIGO.We then restrict such parameter space to be consistent with upper bounds on stochastic GWbackgrounds from BBN and CMB considerations. We also extend our analysis to includepossible changes in the number of relativistic degrees of freedom ( dof ), and the presence ofa small red tensor tilt, as motivated in slow-roll inflation. In Section 5, we summarize ourfindings and discuss briefly the implications for some early Universe scenarios resting uponthe presence of an SD epoch.From now on, m p = 1 / √ πG (cid:39) . · GeV is the reduced Planck mass, a ( τ ) is the scalefactor, τ ≡ (cid:82) dta ( t ) is the conformal time, and we use the Friedman-Lemaˆıtre-Robertson-Walker– 3 –FLRW) metric ds = a ( τ ) η µν dx µ dx ν as the background metric. A tensor-perturbed FLRW metric can be written as ds = a ( τ ) (cid:2) dτ − ( δ ij + h ij ) dx i dx j (cid:3) , (2.1)where we assume the perturbations to be transverse and traceless ∂ j h ij = h ii = 0, so thatthey can be identified with gravitational waves (GWs). Expanding the Einstein equations tolinear order O ( h ) gives h (cid:48)(cid:48) ij + 2 a (cid:48) a h (cid:48) ij − ∇ h ij = 0 , (2.2)where primes denote derivatives with respect to the conformal time τ . To bring (2.2) to amore useful form, we perform a spatial Fourier- and polarization-mode decomposition h ij ( τ, x ) = (cid:88) λ (cid:90) d k (2 π ) e i k . x (cid:15) λij ( k ) h λ k ( τ ) , (2.3)where λ = + , × stands for the polarization states and (cid:15) λij ( k ) are a basis of polarization tensorssatisfying (cid:15) λij ( k ) = (cid:15) λji ( k ), (cid:15) λii ( k ) = 0, k i (cid:15) λij ( k ) = 0, (cid:15) λij ( k ) = (cid:15) λ ∗ ij ( − k ), and (cid:15) λij ( k ) (cid:15) σ ∗ ij ( k ) = 2 δ λσ .These bring us to the GW equation of motion h (cid:48)(cid:48) k + 2 a (cid:48) a h (cid:48) k + k h k = 0 , (2.4)where we have suppressed the polarization indices λ , as we assume that the GW spectrumis unpolarized (cid:10) | h + k | (cid:11) = (cid:10) | h × k | (cid:11) ≡ (cid:104)| h k |(cid:105) . Assuming that the background metric is isotropic,we also write h k = h k , with k ≡ | k | .A useful quantity to characterize a GW background is the tensor power spectrum∆ h ( τ, k ), defined through the following relation (cid:10) h ij ( τ, x ) h ij ( τ, x ) (cid:11) ≡ (cid:90) dkk ∆ h ( τ, k ) ⇐⇒ ∆ h ( τ, k ) = 2 k π (cid:10) | h k ( τ ) | (cid:11) , (2.5)where (cid:104) ... (cid:105) denotes a statistical ensemble average. In this paper we focus on the tensor modes generated during inflation as they were spatiallystretched past the inflationary Hubble radius. At the end of inflation, these modes representsuperhorizon tensor perturbations with an almost scale invariant power spectrum [23]∆ h, inf ( k ) (cid:39) π (cid:18) H inf m p (cid:19) (cid:18) kk p (cid:19) n t , (2.6)with n t a spectral tilt, k p a pivot scale, and H inf the Hubble rate when the mode k p exitedthe horizon during inflation. The presence of the tilt stems from the fact that the inflationary– 4 –hase cannot be perfectly de Sitter . Nevertheless, the spectrum is expected to be only slightlyred-tilted in slow-roll inflation, with the spectral index ’slow-roll suppressed’ as n t (cid:39) − (cid:15) (cid:39) − r p , (2.7)where r p is the tensor-to-scalar ratio evaluated at the scale k p , constrained by the most recentanalysis of the B-mode polarization anisotropies of the CMB at a scale k p /a = 0 .
002 Mpc − ,as r p ≤ .
064 [1, 2]. This bound actually implies an upper bound on the energy scale ofinflation, which must be constrained as H inf (cid:46) H max (cid:39) . · GeV . (2.8)Furthermore, the upper bound on r p also implies that the red-tilted spectral index canonly be very small − n t ≤ . (cid:28)
1. The tensor spectrum is therefore very close to beexactly scale-invariant, at least at around the CMB scales. Actually, in the absence ofrunning of the spectral index, the amplitude of the tensor spectrum would fall only by afactor ∼ (10 ) − . ∼ . e ) ∼
26 orders of magnitude separating theCMB scales and the Hubble radius at the end of inflation. Therefore, for simplicity, we willconsider from now on an exact scale-invariant inflationary spectrum, as this gives an excellentapproximation. We will comment on deviations from this assumption in Sect. 4.3.1.From a theoretical perspective, it is convenient to work with the power spectrum ∆ h ( k ),as it is precisely this quantity that is predicted by inflation to be approximately scale invari-ant. During the evolution of the Universe after inflation, when the tensor modes cross insidethe Hubble radius, they become a stochastic background of gravitational waves (GWs). Inorder to quantify the ability of GW direct detection experiments to measure the inflationaryGW background, it is customary to express the amount of GWs in terms of their energy den-sity spectrum (at sub-horizon scales) Ω GW , defined as the GW energy density ρ GW per unitlogarithmic comoving wavenumber interval, normalized to the critical density ρ crit = 3 m p H [23], Ω GW ( τ, k ) ≡ ρ crit dρ GW ( τ, k ) d ln k = k a ( τ ) H ( τ ) ∆ h ( τ, k ) , (2.9)It is customary to factorize the tensor power spectrum at arbitrary times as a function of theprimordial inflationary spectrum ∆ h, inf ( k ) [c.f. Eq. (2.6)] by means of a transfer function [9]∆ h ( τ, k ) ≡ T h ( τ, k )∆ h, inf ( k ) , T h ( τ, k ) ≡ (cid:18) a k a ( τ ) (cid:19) , (2.10)which characterizes the expansion history between the moment of horizon re-entry τ = τ k of a given mode k , defined as a k H k ≡ k with a k ≡ a ( τ k ), H k ≡ H ( τ k ), and a later moment τ > τ k . For the power spectrum today we will use the notation T h ( k ) ≡ T h ( τ , k ). Note thatthe factor in Eq. (2.10) is simply due to averaging over harmonic oscillations of the modesdeep inside the horizon.If we assume that immediately after inflation, the Universe became radiation dominated(RD) with equation of state w = 1 /
3, the resulting present-day GW energy density spectrumis (quasi-)scale invariant, for the frequency range corresponding to the modes crossing theHubble radius during RD. Setting n t = 0 and averaging over oscillations, the amplitude of– 5 –he plateau characterizing the energy density spectrum today isΩ (0)GW (cid:12)(cid:12)(cid:12) plateau (cid:39) G k Ω (0)rad π (cid:18) H inf m p (cid:19) (cid:39) · − (cid:18) H inf H max (cid:19) , (2.11)where we have used k = a k H k and introduced the RD transfer function [9] T h ( k ) (cid:39) (cid:18) a k a (cid:19) (cid:39) G k Ω (0)rad (cid:18) a H a k H k (cid:19) , G k ≡ (cid:18) g ∗ ,k g ∗ , (cid:19) (cid:18) g s, g s,k (cid:19) / . (2.12)In the rhs of Eq. (2.11) we have used Ω (0)rad (cid:39) · − , g s, (cid:39) . g ∗ , = 3 .
36, and g s,k (cid:39) g ∗ ,k (cid:39) .
75 (so that G k (cid:39) . g s,k , g ∗ ,k equal to theStandard Model (SM) degrees of freedom ( dof ) before the electroweak symmetry breaking,and hence independent of k . In reality the number of SM relativistic dof change with thescale, but we postpone the discussion of this spectral distortion to Section 4.3.2. For thetime being we simply consider an identical suppression for all the modes as G k ∼ . plateau of the inflationary GW (quasi-)scaleinvariant energy density spectrum today, corresponding to the modes that crossed the horizonduring RD. However, if after inflation there is a transient period of evolution with EoS w (cid:54) = 1 /
3, before RD is established, the resulting GW energy density spectrum today will nolonger remain scale-invariant. As we will see next, the spectrum today will actually consistof two parts: a high-frequency branch, corresponding to the modes that crossed the horizonduring the transient epoch, and a (quasi-)scale invariant branch corresponding to the modesthat crossed the horizon during RD . Let us consider now that there is a period in the early Universe, spanning from the end ofinflation till the onset of RD, with EoS w (cid:54) = 1 / V ∝ φ n ,an effective (oscillation averaged) EoS emerges as w (cid:39) ( n − / ( n + 1) [15]. For n = 2( V ∝ φ ) we obtain a RD period with w (cid:39) /
3, while for n = 1 ( V ∝ φ ) we obtain instead amatter dominated (MD) era with w (cid:39)
0. Interestingly, for n ≥ stiff dominated (SD) period with EoS w > / − / < w <
1. Even though it is common to assume that 0 ≤ ¯ w ≤ /
3, thereis a priori no reason (theoretical or observational) to exclude the stiff case 1 / < w < V (cid:29) K , where V and K are the inflaton potential and kinetic energy densities.Inflation cannot be sustained however if the potential drops to V < K/
2. Furthermore, if a There is yet another part of the spectrum, corresponding to modes that crossed the Hubble radius aftermatter-radiation equality, which behaves as Ω ∝ /k . This corresponds to very small frequencies today f (cid:46) − Hz, and hence we will not be concerned with such low frequency end of the spectrum, as it onlyaffects the CMB and it cannot be probed by direct-detection GW experiments. – 6 –eature in the inflaton potential allows its value V to drop much below the kinetic energy K ,the EoS can become stiff after inflation, w = ( K − V ) / ( K + V ) > / V (cid:29) K during inflation to some small value V (cid:28) K after inflation. Thetransition itself would actually trigger the end of inflation, leading to a post-inflationary EoS w (cid:39) −O ( V /K ). In general we expect that the EoS can approach unity from below, but neverachieve w = 1 exactly, as this would require an exactly flat direction with V = 0. A naturalscenario where inflation is followed by a KD phase is that of Quintessential-Inflation [25],where the inflaton potential V ( φ ) is engineered so that the necessary transition occurs atthe end of inflation, and the potential is also adjusted to describe the observed dark energyas a quintessence field, see e.g. [26–33] for different proposals. As mentioned before, a stiffperiod can be also engineered through the oscillations of the inflaton with potential V ∝ φ n , n ≥
3. In this case we note however that the stiff period cannot be sustained for very long,as self-resonant effects lead eventually to a fragmentation of the coherent oscillating inflatoncondensate [53].Considering the presence of an SD period before the onset of RD is actually not onlytheoretically well motivated but also phenomenologically interesting. As we have mentionedin the Introduction, and as we will show in detail, an SD period with equation of state w > / ρ tot = ρ ∗ exp {− (cid:82) daa ( w ( a ) + 1) } , with ρ ∗ = 3 m p H ∗ theinitial energy density at the end of inflation. We will distinguish the Hubble rate at the endof inflation H ∗ from that during inflation H inf throughout the text, except in Section 4.2and 4.3 where we take H inf (cid:39) H ∗ in deriving the numerical results. In general the EoS isdetermined by the inflaton potential and is a function of time. However we expect it tochange only adiabatically during SD, and in any case we can always describe the scaling ofthe energy density in terms of an effective (logarithmic-averaged) value of the EoS duringthe stiff period, 3 (cid:82) daa (1 + w ( a )) ≡ w ) log( a/a ∗ ), so that ρ tot = ρ ∗ ( a/a ∗ ) − w +1) .Let us consider therefore an SD period between the end of inflation and the onset ofRD, with effective EoS w ≡ w s > / dof , according to the standard hot Big Bangpicture. We sketch the different epochs of the expansion history we consider in Figure 1.From now on, the subscripts * , RD , and BBN , stand for “evaluation at” or “corresponding to”the end of inflation/beginning of SD epoch, end of SD epoch/beginning of RD period, andonset of BBN, respectively. In order to not sabotage the success of BBN, a minimal requisitethat we need to impose over the assumed expansion history is that the stiff epoch must endbefore the beginning of BBN, i.e. τ RD < τ BBN .In this section we will solve (2.4) for the aforementioned cosmological scenario, propagatethe solution to the present era when it can be detected, and express the result in terms of theGW energy density spectrum today (2.9). Our focus is mainly on the modes that re-enter– 7 – inf τ BBN τ RD now τ eq w < − w = w S w = w = 0 w = − inflation SD RD MD ΛD τ Figure 1 . ΛCDM+inflation expansion history with a stiff epoch. the horizon during the stiff phase (SD-reentering modes), as they constitute the part of theGW spectrum that is enhanced and hence potentially detectable. The expression for thepresent-day GW energy spectrum depends on the assumption on how the transition from theSD to RD era takes place. There are two natural cases to consider:1.
Instant transition.
Here the transition from SD to RD is modeled as ’instantaneous’,i.e. it occurs in a very short time interval compared to the Hubble timescale at themoment of the transition. The scale factor a ( τ ) and Hubble rate H ( τ ) are continuousduring the transition, but we consider a sudden jump in the effective EoS from w = w s > / w = 1 /
3. For example, this could happen if the stiff fluid decays into radiationat some point (which marks the end of the SD epoch) by some process characterizedby a timescale much shorter than the instantaneous Hubble time at that moment.2.
Smooth transition.
Here the SD spoch is driven by a fluid component (typicallyformed by the inflaton field itself) which dominates the energy budget of the Universeand has a stiff equation of state w s > /
3, but there is also a relatively small amount ofradiation present at the end of inflation / onset of SD. The energy densities of the twocomponents scale freely as the Universe expands, i.e. we assume no interaction betweenthe two sectors. Due to the different scalings, as ∝ a − for radiation and as ∝ a − w s ) for the stiff fluid, the radiation component eventually dominates the energy budget ofthe Universe. The transition from SD to RD then occurs smoothly, over a few Hubbletimes, as the energy density of radiation gradually overtakes that of the stiff fluid.Since the time-dependence of the scale factor around the SD-to-RD transition is differentin the two cases, the GWs that have entered the horizon during SD will evolve differently.In turn, this means that the GW energy spectrum for the modes that re-entered the horizonbefore and around the SD-to-RD transition will differ in the two cases. The GW energyspectrum can be computed fully analytically in the instant transition case. However, inthe smooth transition case we can only compute analytically the asymptotic high and lowfrequency branches of the GW spectrum, corresponding to the modes that entered the horizondeep inside SD and RD, respectively. From an observational point of view, this is not aproblem, as only the high frequency branch of the spectrum can be potentially probed byGW detectors. For completeness, in any case, we will provide a numerical computation ofthe full GW spectrum in the smooth transition case.– 8 – .1 Instant transition We consider in this subsection the expansion history shown in Figure 1, assuming the tran-sition from the SD to RD epoch occurs instantaneously, that is, in a time much shorter thanthe instantaneous Hubble time at the moment of transition. We consider the scale factorand the Hubble rate around the transition as continuous smooth functions, but we modelthe EoS as a discrete function¯ w = w s Θ( τ RD − τ ) + 13 Θ( τ − τ RD ) , (3.1)where Θ( x ) is the step-function. The energy density of the background is therefore continu-ous, and scales as ρ tot = ρ ∗ ( a/a ∗ ) − w s ) , τ ≤ τ RD (Stiff Domination) ρ RD ( a/a RD ) − , τ ≥ τ RD (Radiation Domination) (3.2)where ρ RD ≡ ρ ∗ ( a RD /a ∗ ) − w s ) . From here the scale factor during each period can then besolved exactly as a ( τ ) = a ∗ (cid:20) a ∗ H ∗ α s ( τ − τ ∗ ) (cid:21) α s = a ∗ (cid:18) a ∗ H ∗ α s (cid:19) α s [˜ τ s ( τ )] α s , τ ∗ ≤ τ ≤ τ RD (3.3) a ( τ ) = a RD [1 + a RD H RD ( τ − τ RD )] = a H RD ˜ τ r ( τ ) , τ RD ≤ τ (cid:28) τ eq (3.4)and, correspondingly, the Hubble rate aH ( τ ) ≡ d log adτ as aH ( τ ) = a ∗ H ∗ a ∗ H ∗ α s ( τ − τ ∗ ) = α s ˜ τ s ( τ ) , τ ∗ ≤ τ ≤ τ RD (3.5) aH ( τ ) = a RD H RD a RD H RD ( τ − τ RD ) = 1˜ τ r ( τ ) , τ RD ≤ τ (cid:28) τ eq , (3.6)where τ eq denotes the time at the radiation-matter equality, and we have defined α s ≡
21 + 3 w s . (3.7)For convenience, we have introduced the time variables ˜ τ s and ˜ τ r via a ∗ H ∗ α s ˜ τ s ( τ ) ≡ a ∗ H ∗ α s ( τ − τ ∗ ) , a RD H RD ˜ τ r ( τ ) ≡ a RD H RD ( τ − τ RD ) . (3.8)From now on, and without loss of generality, we fix a ∗ = 1 and τ ∗ = 0 at the end of inflation(onset of SD).During the stiff epoch the GW equation of motion (2.4) reads, using (3.5), h (cid:48)(cid:48) k + 2 α s ˜ τ s h (cid:48) k + k h k = 0 , (3.9)where (cid:48) denotes derivatives with respect to ˜ τ s (since d ˜ τ s = dτ , we do not distinguish betweenthe derivatives with respect to ˜ τ s or τ ). Requiring that the tensor mode function must match– 9 –he inflationary spectrum h k ( τ ) = h inf k and h (cid:48) k ( τ ) = 0 in the superhorizon limit kτ (cid:28)
1, thesolution to (3.9) during the stiff period when τ ∗ ≤ τ ≤ τ RD is h (stiff) k ( τ ) = Γ (cid:18) α s + 12 (cid:19) (cid:18) α s y (cid:19) α s − J α s − ( α s y ) h inf k , y ≡ kaH = k ˜ τ s ( τ ) α s , (3.10)where J ν ( x ) is the Bessel function of the first kind. Using the small argument limit of theBessel function, J ν ( x ) (cid:39) x ν / (2 ν Γ( ν +1)) for x (cid:28)
1, we obtain h (stiff) k ( τ ) (cid:39) h inf k when k (cid:28) aH ,as it should. The superscript (stiff) indicates that the solution applies in the stiff epoch. Wewill use analogous notations in what follows.Using (3.6), the GW equation of motion (2.4) during the the RD epoch reads h (cid:48)(cid:48) k + 2˜ τ r h (cid:48) k + k h k = 0 , (3.11)where this time (cid:48) denotes derivatives with respect to ˜ τ r (again simply because d ˜ τ r = dτ ).The solution during the RD era to (3.11) τ RD ≤ τ (cid:28) τ eq , is h (rad) k ( τ ) = 1 √ y (cid:104) A ( k ) J ( y ) + B ( k ) Y ( y ) (cid:105) , y ≡ kaH = k ˜ τ r ( τ ) , (3.12)where J / ( y ) , Y / ( y ) are Bessel functions of the first and second kind and the superscript (rad) indicates that the solution applies during RD. At τ = τ RD this solution must match withsolution (3.10) [and hence simultaneously with h inf k if the mode is superhorizon]. Continuityof the tensor modes and their derivatives requires h (stiff) k ( τ RD ) = h (rad) k ( τ RD ) , h (cid:48) (stiff) k ( τ RD ) = h (cid:48) (rad) k ( τ RD ) (3.13)from which we get (cid:40) A ( k ) B ( k ) (cid:41) = 1 √ (cid:18) α s (cid:19) α s Γ (cid:18) α s + 12 (cid:19) κ − α h inf k × (cid:40) a ( κ ) b ( κ ) (cid:41) , (3.14)where κ ≡ kk RD , (3.15)with k RD ≡ a RD H RD , and a ( κ ) = √ α s (cid:32) J α s − ( α s κ ) Y ( κ ) − J α s + ( α s κ ) Y ( κ ) J ( κ ) Y ( κ ) − J ( κ ) Y ( κ ) (cid:33) , (3.16) b ( κ ) = √ α s (cid:32) J α s + ( α s κ ) J ( κ ) − J α s − ( α s κ ) J ( κ ) J ( κ ) Y ( κ ) − J ( κ ) Y ( κ ) (cid:33) . (3.17)The sub-horizon kτ (cid:29) τ (cid:29) τ RD limit of (3.12) is h (rad) k ( τ ) = (cid:114) πy [ A ( k ) sin y − B ( k ) cos y ] , (3.18)– 10 –here we used the large argument expansion of Bessel functions J ν ( x (cid:29) (cid:39) (cid:113) πx sin( x − δ ν ), Y ν ( x (cid:29) (cid:39) (cid:113) πx cos( x − δ ν ) and ˜ τ ≈ τ . Substituting (3.14) into (3.18), squaring, andaveraging over mode oscillations, we obtain | h (rad) k ( τ ) | = 12 πy (cid:18) α s (cid:19) α s Γ (cid:18) α s + 12 (cid:19) κ − α s ) W ( κ ) | h inf k | (3.19)= (cid:18) a RD a ( τ ) (cid:19) π (cid:18) α s (cid:19) α s Γ (cid:18) α s + 12 (cid:19) κ − α s W ( κ ) | h inf k | , (3.20)where we have used y = κ ( a/a RD ) in the second line, and defined W ( κ ) ≡ a ( κ ) + b ( κ ) = πα s κ (cid:20)(cid:16) κJ α s + ( κ ) − J α s − ( κ ) (cid:17) + κ J α s − ( κ ) (cid:21) . (3.21)It can be shown that subhorizon modes always scale as h k ∝ a − regardless of how thescale factor a evolves with time. Thus, even though the second expression of | h (rad) k ( τ ) | c.f. Eq. (3.20), was derived in the RD epoch, once its time evolution (i.e. the damping ofthe tensors due to the expansion of the Universe) is written as ∝ a − ( τ ) the expressionremains valid for the subsequent epochs as well. This is, however, only true for modes thatbecame subhorizon before the moment of matter-radiation equality at τ = τ eq , i.e. for modes k (cid:29) k eq ≡ a eq H eq .Building the present-day tensor power spectrum ∆ h ( τ , k ) = k π | h (rad) k ( τ ) | with (3.20)and plugging this into (2.9), leads to the present-day energy spectrum for the modes k (cid:29) k eq re-entering the horizon during the SD or RD epochs,Ω (0)GW ( k ) ≡ k ∆ h ( τ , k )12 a H = a k πa H (cid:18) α s (cid:19) α s Γ (cid:18) α s + 12 (cid:19) κ − α s W ( κ ) ∆ h, inf ( k )= (cid:18) a RD a (cid:19) (cid:18) H RD H (cid:19) π (cid:18) H inf m p (cid:19) Γ ( α s + 1 / − α s ) α α s s Γ (3 / W ( κ ) κ − α s ) , (3.22)where in the last step we have introduced inflationary tensor power spectrum (2.6) (with n t = 0), and used κ ≡ k/k RD , k RD = a RD H RD , and π = 4Γ (3 / ρ rad ( τ RD ) = ρ crit ( τ RD ) = 3 m p H . This and the scaling law of radiationenergy density implies (cid:18) a RD a (cid:19) (cid:18) H RD H (cid:19) = 8 πGρ rad ( τ )3 H = Ω (0)rad (cid:18) g ∗ ,k g ∗ , (cid:19) (cid:18) g s, g s,k (cid:19) / . (3.23)Plugging (3.23) into Eq. (3.22), using Eq. (2.11) for the inflationary plateau, and expressingthe result as a function of present-day frequencies f = k/ (2 πa ), we finally obtainΩ (0)GW ( f ) = Ω (0)GW (cid:12)(cid:12)(cid:12) plateau × W ( f /f RD ) × A s (cid:18) ff RD (cid:19) − α s ) , (3.24) As we will see later on, an analogous relation in the smooth transition case differs by a factor of 2. – 11 –here f RD ≡ k RD / (2 πa ) is the frequency corresponding to the horizon scale at the onsetof RD, k RD = a RD H RD , W ( x ) is the window function defined in Eq. (3.21), and we haveintroduced the constant A s ≡ Γ ( α s + 1 / − α s ) α α s s Γ (3 / , (3.25)which ranges as 1 < A s < /π (cid:39) .
27 for 1 / < w s <
1. The window function W ( x ) variessmoothly around the frequencies f ∼ f RD , and its asymptotic limits at large frequencies f (cid:29) f RD (corresponding to modes crossing the horizon during SD) and small frequencies f (cid:28) f RD (corresponding to modes crossing the horizon during RD), determine the asymptoticbehaviour of the energy density spectrum. In particular we obtain W ( f /f RD (cid:28) −→ A − (cid:18) ff RD (cid:19) − − α s ) , W ( f /f RD (cid:29) −→ , (3.26)and hence Ω (0)GW ( f ) (cid:39) Ω (0)GW (cid:12)(cid:12)(cid:12) plateau × , f (cid:28) f RD A s (cid:16) ff RD (cid:17) − α s ) , f (cid:29) f RD . (3.27)What matters from the point of view of detection prospects of this signal is the fact that thehigh-frequency branch of the spectrum rises with frequency, exhibiting a significant blue tiltfor a stiff EoS w S > / n t ≡ d log Ω (0)GW d log f = 2(1 − α s ) = 2 (cid:18) w S − w S + 1 (cid:19) > , (3.28)which approaches unity n t −→ w S −→
1. It is precisely this large tilt thatleads us to consider the ability of GW detectors to measure this signal: as we will discusslater, a significant fraction of the parameter space characterizing the shape of the spectrum, { w S , f RD , H inf } leads to the high-frequency part of the spectrum being above the sensitivityof LISA and LIGO at their corresponding key frequencies.The window function characterizes, in a sense, the ’interpolation’ around the centralfrequencies f ∼ f RD , of the two asymptotic regimes at large f (cid:29) f RD and small f (cid:28) f RD frequencies. In the next section, we will actually compute numerically the exact frequencydependence of the window function W ( f ) when the SD-to-RD transition is not modeledas instantaneous, but rather a smooth transition resulting from the gradual dominationof the energy budget in the Universe of an initially small radiation component. From anobservational point of view, given the current upper bound on the amplitude of the smallfrequency plateau [c.f. Eq. (2.11)], the actual frequency dependence of W ( f ) around f ∼ f RD is irrelevant, as it cannot be observed by direct detection experiments. Nevertheless, as theexpansion history of the Universe changes depending on the expansion rate assumed aroundthe SD-to-RD transition, the high frequency branch f (cid:29) f RD of the GW spectrum willexperience a slightly different expansion history once the corresponding modes cross insidethe horizon. As we will show next, this translates into a correction of the normalizationconstant A s characterizing the rising high-frequency branch of the spectrum.Note that the late-time acceleration correction factor (Ω m / Ω Λ ) pointed out in [54] doesnot apply to our derivation. The correction factor is only needed if we adopt the commonly– 12 –sed fitting formula which can be traced back to Eq. (16) of [55]. This fitting formula isobtained by solving (2.4) numerically in the presence of matter and radiation, but withoutcosmological constant, and then fitting the resulting transfer function. The fitted transferfunction needs to be corrected by the factor (Ω m / Ω Λ ) simply because the cosmologicalconstant was not taken into account in getting the transfer function (see e.g. Section 3 of[55] or the Appendix of [56]). Now, in our derivation all the information about the expansionhistory between the horizon reentry τ k and the present time τ , including that of the latecosmological-constant dominance, is contained in the ratio a /a k . Thus, it would be incorrectto include the said correction factor. If we set the scale factor today a = 1 and specified agiven mode k , the correction factor is just the ratio between the scale factor at the moment ofhorizon crossing a k for a universe with cosmological constant and that for a universe withoutcosmological constant. Since we are always considering a complete universe with cosmologicalconstant included, there is nothing to correct for in our derivation. Let us consider now a situation where at the end of inflation τ = τ ∗ , the total energy density ρ ∗ ≡ m p H ∗ is split between a dominant stiff fluid with energy density ρ ∗ stiff = (1 − (cid:15) ) ρ ∗ , (cid:15) (cid:28) ρ ∗ rad = (cid:15)ρ ∗ (cid:28) ρ ∗ stiff . We assumethat the stiff fluid has an equation of state parameter w s > /
3, which we take as constantfor simplicity. The expansion of the Universe is then driven by the energy densities of thetwo fluids, each of which scale freely, as ρ rad ∝ a − and as ρ stiff ∝ a − w S ) . The outlineof the expansion history will roughly follow the sketch shown in Figure 1, but this time, thetransition from SD to RD is slow and smooth, instead of ’instantaneous’. While the stifffluid dominates, ρ stiff (cid:29) ρ rad , the equation of state of the Universe remains approximatelyconstant and equal to w s . Since ρ stiff scales down faster than ρ rad , there is always a moment,which we denote as τ = τ RD , at which ρ rad ( τ = τ RD ) = ρ stiff ( τ = τ RD ). This point marksthe end of the SD epoch and the beginning of the RD epoch, although the expansion historyaround this moment is neither purely SD nor RD, but rather dictated by a mixture of thetwo fluids. In a few Hubble times after τ = τ RD , the radiation component becomes theenergy-dominating fluid, and from then on the expansion history follows that of standardRD embedded in the usual ΛCDM scenario.In such smooth SD-to-RD transition, the scale factor cannot be solved analytically,let alone the GW equation of motion. It is, however, possible to work out analytically theblue-tilted high-frequency branch of the GW spectrum corresponding to modes entering thehorizon deep inside the SD epoch at τ (cid:28) τ RD , when the equation of state of the Universe w (cid:39) w s is approximately constant. Fortunately, these are the modes that can actually beprobed by GW detectors. Far before the SD-to-RD transition takes place, the differencebetween the previously considered instant transition and the presently considered smoothtransition is not yet apparent and the tensor mode function is given by the same expressionas in the instant transition case (3.10). In order to avoid having to deal with the part of theexpansion history close to the SD to RD transition where the evolution of the scale factoris not analytically solvable, we employ the trick we used earlier, namely rewriting the tensorpower spectrum in terms of the scale factor ∝ a − . Once we do that, the resulting expressionwill be valid in all the subsequent epochs.In order to proceed, we need first to obtain the value of a RD in the two fluid approach.The condition ρ rad ( τ RD ) = ρ stiff ( τ RD ) at τ RD implies ρ stiff ( τ RD ) /ρ stiff ( τ ∗ ) = H / H ∗ and thescaling of the energy density of the stiff fluid gives ρ stiff ( τ RD ) /ρ stiff ( τ ∗ ) = ( a RD /a ∗ ) − w s +1) .– 13 –ogether, they yield a RD a ∗ = (cid:18) / H ∗ H RD (cid:19) α s1+ α s . (3.29)Taking the sub-horizon limit of expression (3.10), squaring it, and averaging over oscillations,we arrive at (cid:12)(cid:12)(cid:12) h (stiff) k (cid:29) k RD ( τ ) (cid:12)(cid:12)(cid:12) = 12 π Γ (cid:18) α s + 12 (cid:19) (cid:18) k ˜ τ s (cid:19) α s (cid:12)(cid:12)(cid:12) h inf k (cid:12)(cid:12)(cid:12) , = 12 π Γ (cid:18) α s + 12 (cid:19) (cid:18) α (cid:19) α s κ − α s (cid:18) a RD a ( τ ) (cid:19) (cid:12)(cid:12)(cid:12) h inf k (cid:12)(cid:12)(cid:12) , (3.30)where we recall that κ ≡ k/k RD = f /f RD , a ∗ H ∗ ˜ τ s ( τ ) = α s + a ∗ H ∗ ( τ − τ ∗ ) [c.f. (3.8)], and inthe second line we have used the scale factor a ( τ ) = a α s ∗ H α s ∗ α − α s s [˜ τ s ( τ )] α s deep inside SDduring τ ∗ ≤ τ (cid:28) τ RD [c.f. (3.3)], together with (3.29) and k RD ≡ a RD H RD .Now that the solution is expressed in terms of the scale factor, it remains valid in all thesubsequent epochs, and we can omit the superscript (stiff) . Building the present-day tensorpower spectrum ∆ h ( τ , k ) = k π | h k (cid:29) k RD ( τ ) | with (3.20), and plugging this into (2.9), leadsto the present-day energy spectrum for the modes k (cid:29) k RD re-entering the horizon duringthe SD,Ω (0)GW ( f (cid:29) f RD ) ≡ k ∆ h ( τ , k )12 a H = a k πa H (cid:18) α (cid:19) α s Γ (cid:18) α s + 12 (cid:19) κ − α s ∆ h, inf ( k )= (cid:18) a RD a (cid:19) (cid:18) H RD H (cid:19) π (cid:18) H inf m p (cid:19) Γ ( α s + 1 / − α s α α s s Γ (3 / κ − α s ) = Ω (0)GW (cid:12)(cid:12)(cid:12) plateau × ˜ A s (cid:18) ff RD (cid:19) − α s ) , ˜ A s ≡ − α s A s , (3.31)where in the second step we have introduced the inflationary tensor power spectrum (2.6)(with n t = 0) and used κ ≡ k/k RD , k RD = a RD H RD and π = 4Γ (3 / κ = f /f RD , the definition of A s [c.f. Eq. (3.25)] and of the inflationary plateau [c.f. (2.11)], and the fact that in a smooth transition ρ rad ( τ RD ) = ρ crit ( τ RD ) / H / πG which implies (cid:18) a RD a (cid:19) (cid:18) H RD H (cid:19) = 2 ρ rad ( τ )3 m p H = 2 Ω (0)rad (cid:18) g ∗ ,k g ∗ , (cid:19) (cid:18) g s, g s,k (cid:19) / . (3.32)We notice that in (3.32) there is an extra factor of 2 compared to the analogous expression(3.23) for the instant transition case. As before, in the final expression of Eq. (3.31) weabsorbed the effects due to the changes in the relativistic dof into Ω (0)GW (cid:12)(cid:12)(cid:12) plateau .If we compare the expression of the high-frequency branch of the GW energy spectrumwe just obtained in the smooth transition case with its instant transition counterpart (3.27),we see that the normalization constant is now a factor 2 − α s larger, i.e. a factor that rangesfrom 1 (if w s = 1 /
3) to √ w s → w s (cid:39)
1, the amplitudeof the observable high-frequency branch of the spectrum is actually ∼
40% larger than in theinstant transition case. For comparison we plot in Figure 2 the present GW energy densitypower spectrum obtained in the instant SD-to-RD transition model, c.f. Eq. (3.24), together– 14 – .01 0.10 1 10 10010 - - - - f / f RD h Ω G W ( ) ( f ) Figure 2 . Comparison among different forms of the present-day GW energy density spectra. For theinstant transition case we show the full oscillatory solution (blue solid line) computed with Eqs. (3.33)-(3.34), together with its oscillation envelope (black dashed line) and oscillation averaged analyticalspectrum (black solid line) computed with Eq. (3.24). For the the smooth transition case we plot thefull oscillatory form (red solid line) based on numerical evaluation of Eqs. (3.35)-(3.38), and for its highfrequency branch the oscillation envelope (green dashed line) and the oscillation averaged analyticalspectrum (green solid line) obtained with Eq. (3.31). For each type of transition the jump betweenthe oscillation envelope (dashed lines) and oscillation average (solid lines) is a factor 2 exactly, asthe tensor mode functions exhibit exact harmonic oscillations at deep inside sub-horizon scales. Thedifference in amplitude in the high frequency branch between the smooth and the instant transitioncases, correspond to a factor 2 − α s (cid:39) .
41, as obtained for w s = 0 . with the high frequency branch obtained in the smooth transition case, c.f. Eq. (3.31). Forcompleteness, we also plot the GW spectrum without averaging over oscillations. In theinstant transition this corresponds toΩ (0)GW ( f ) = Ω (0)GW (cid:12)(cid:12)(cid:12) plateau × A s κ − α s ) × W osc ( κ ) , κ ≡ kk RD = ff RD (3.33) W osc ( κ ) ≡ πy (cid:104) a ( κ ) J ( y ) + b ( κ ) Y ( y ) (cid:105) , y ≡ kaH . (3.34)In a smooth transition we need to obtain the spectrum fully numerically asΩ (0)GW ( f ) = Ω (0)GW (cid:12)(cid:12)(cid:12) plateau × A s κ − α s ) × W num ( κ ) , κ ≡ kk RD = ff RD (3.35) W num ( κ ) ≡ πκ (cid:18) a ( τ ) a RD (cid:19) |H k ( τ ) | , (3.36)with a ( τ ) and H k ( τ ) the solution to the differential equations a (cid:48) ( τ ) = a ( τ ) H ∗ (cid:16) (1 − δ ) a ( τ ) − w s ) + δ a ( τ ) − (cid:17) / , a ( τ ∗ ) = 1 (3.37) H (cid:48)(cid:48) k ( τ ) + 2 a (cid:48) ( τ ) a ( τ ) H (cid:48) k ( τ ) + k H k ( τ ) = 0 , (cid:26) H k ( τ = (cid:15)/k ) = 1 , H k ( τ = (cid:15)/k ) = 0 , , (3.38)– 15 –here (cid:15) (cid:28) δ ≡ ρ ∗ rad /ρ ∗ is the initial fraction of the radiationenergy density. We observe that if we average the expression of W osc ( κ ) from Eq. (3.34) overmode oscillations, we recover the expression for W ( κ ) from Eq. (3.21), as it should. Regardless of how we model the SD-to-RD transition, the overall shape of the GW energyspectrum today h Ω GW ( τ , f ) can be characterized by three parameters: the inflationaryHubble rate H inf , the equation of state parameter deep inside the SD epoch w s , and thefrequency today f RD corresponding to the horizon scale at the transition from SD to RD.The spectra exhibit in all cases a low frequency plateau at f (cid:28) f RD , and a power law branch ∝ f − α s ) at large frequencies f (cid:29) f RD . In the top and left-bottom panels in Figure 3,we plot the oscillation averaged spectra in the instant transition case, evaluated at differentvalues of the parameters { w s , H inf , f RD } . As can be seen in these panels, H inf controls thelevel of the plateau (and hence the overall amplitude of the full spectrum, see the left-toppanel), whereas w s determines the slope of the high-frequency part, and f RD determines thelocation of the “elbow” where the plateau and the blue-tilted parts are connected. In thementioned panels in Figure 3, we also plot the power-law sensitivity curves for stochastic GWbackgrounds of LISA and the O1, O2 and O5 runs of advanced LIGO. As it is evident fromthe figure panels, there are values of the parameters for which we expect the signal to beclearly observable by LIGO and LISA. In the following, we will first determine the parameterspace { w s , H inf , f RD } compatible with a detection and later subtract from it the region thatis incompatible with current upper bounds on stochastic GW backgrounds, as set by BBNand CMB constraints.Had we also plot in Figure 3 the high-frequency branch of the GW oscillation-averagedenergy spectra in the smooth-transition case, the difference between the instant- and smooth-transition spectra would be barely noticeable to the eye (given the logarithmic scales in-volved). For convenience, whenever we plot a GW spectrum against frequency today, wewill use the instant transition case, as only then we have an analytic expression for the fullspectral range, all the way from small f (cid:28) f RD to large f (cid:29) f RD frequencies, c.f. Eq. (3.24).For the analysis we present in the rest of this section, however, we could use either case, aswe only need to evaluate the high-frequency branch of the spectrum, for which we have theanalytic expressions for both instant- and smooth- transitions, recall Eqs. (3.27) and (3.31).As an instantaneous SD-to-RD transition can at best be an approximation, we consider thesmooth-transition case as more realistic. For the following parameter analysis we will usetherefore the smooth transition case.To start with, we limit our parameter-space scan within the following direct boundsimposed on each of the three parameters: the non-detection of primordial B-modes in theCMB puts an upper bound on the energy scale of inflation H inf ≤ H max (cid:39) . × GeV,c.f. Eq. (2.8); the EoS of an SD epoch must lie within the range 1 / < w s < f RD ≥ f BBN , where f BBN is the red-shifted frequency today corresponding to the horizon scale at the onset of– 16 – - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - f [ Hz ] h Ω G W Figure 3 . Top-left, top-right, bottom-left panels : GW energy spectra in the instant transition case(black solid lines), together with the LISA sensitivity curve (blue solid line), and LIGO sensitivitycurves (green solid lines). We also indicate the BBN upper bound (red dotted line), which we introducelater in Sect. 4.2. In the top-left panel we fix w s = 0 . f RD = 10 − Hz and plot the GW energyspectra for H inf = 10 GeV, 10 GeV, 10 GeV. In the top-right panel we fix H inf = 10 GeV and w s = 0 . f RD = 10 − Hz, 10 − Hz, 10 − Hz. In the bottom-leftpanel we fix H inf = 10 GeV and f RD = 10 − Hz and plot the GW energy spectra for w s = 0 .
5, 0 . Bottom-right panel : GW energy spectra in the instant transition case taking into account changesin the number of relativistic dof , depending on whether the SD-to-RD transition takes place before(red dashed line) or after (blue solid line) the QCD phase transition. For comparison we show thecorresponding spectra (gray solid lines) without correcting for changes in the number of dof . Wepostpone the discussion on this to Sect. 4.3.2.
BBN, f BBN ≡ π a BBN a H BBN
GeV × . · Hz (4.1)= Ω (0)rad π (cid:114) H Hz H BBN
GeV × . · Hz (cid:39) . · − Hz , (4.2)where we have used Ω (0)rad (cid:39) . · − , g s, BBN = g s, (cid:39) . g ∗ , BBN = g ∗ , = 3 . H (cid:39) . × . · − Hz, and H BBN = π ( g ∗ , BBN / / ( T /m p ), with T BBN = 0 .
001 GeV.
While expressing the GW energy spectrum in terms of f RD yields a neat expression, it is moreuseful from a physical point of view to characterize the point of the SD-to-RD transition in– 17 –erms of an energy scale, namely the temperature T RD of the radiation dof at that moment .For a smooth SD-to-RD transition T RD ≡ (cid:18) g ∗ , RD π ρ rad ( τ RD ) (cid:19) / = (cid:18) g ∗ , RD π m p H (cid:19) / = G − / g − / ∗ , RD × √ π √ (cid:113) Ω (0)rad (cid:114) m p H × . · − (cid:18) f RD Hz (cid:19) GeV (cid:39) G − / g − / ∗ , RD × . · (cid:18) f RD Hz (cid:19) GeV (4.3)where we have used where we have used f RD = 2 πa RD k RD , k RD = a RD H RD , (3.32), andin the last line Ω (0)rad (cid:39) × − and H (cid:39) . rhs of Eq. (4.3) would simply be multiplied by √
2. Eq. (4.3) informs usthat the energy scales of, e.g. the electroweak phase transition ∼
100 GeV, BBN ∼ ∼ . ∼ − Hz, ∼ − Hz, and ∼ − Hz, respectively. It also says that the typical frequenciesthat LISA and LIGO can probe, namely f LISA ∼ . f LIGO ∼
10 Hz, correspond toenergy scales ∼ GeV and ∼ GeV, respectively. The requirement that SD ends beforeBBN starts limits the possible values of the energy scale to T RD (cid:38) H inf and fixed- f RD (or T RD ) slices of the parameter-spaceregions, probe-able by LISA and the O2 run of aLIGO. For the time being, in the presentanalysis, and until the end of Section 4.2, we fix n t = 0 and the number of dof to g ∗ , RD = g s, RD (cid:39) .
75. In other words, the effects of a red-tilt in the GW spectrum and of changesin the number of relativistic dof in the RD epoch are not included in these plots. We willpresent corrections to our analysis due to the inclusion of these effects in Sections 4.3.1 and4.3.2, respectively. In light of the left panels in Figure 4, we find that in order to detect theGW background by LISA, the values of the parameters must lie in the following ranges10 GeV (cid:46) H inf ≤ . × GeV , . (cid:46) w s < , − Hz (cid:46) f RD < . × − Hz , − GeV (cid:46) T RD < . × GeV , which are understood as follows. If a parameter, say H inf , lies in the above specified range,then there are values of the other two parameters, i.e. w s and f RD ( T RD ), that give rise toGW signals that are detectable by LISA.Let us note that as we scan the parameters H inf , w s , and f RD from somewhere outside theprobe-able region by LISA, the GW spectrum first intersects the LISA sensitivity curve close The radiation components in our case are likely to be in thermal equilibrium before the SD-to-RD tran-sition, but this is not guaranteed. To avoid this problem, we could define instead an energy scale associatedto f RD , as the 4 th root of the radiation energy density E RD = ρ rad ( τ RD ) / at τ = τ RD . However, since inthe hot Big Bang picture it is customary to parametrize the energy scale of the RD era by the temperatureof the thermal radiation, we will stick to this convention. We will simply characterize the energy scale by atemperature parameter T RD , related to the radiation energy density as ρ rad ( τ RD ) ≡ π g ∗ , RD T , regardless ofwhether the relativistic degrees of freedom are in thermal equilibrium. For all cases of (observational) interest,we expect that the SM degrees of freedom in the radiation bath will have reached thermal equilibrium beforethe transition time τ = τ RD . – 18 –ISA LIGO O2 - - - - - - - - - - f RD ( Hz ) ω s T RD ( GeV ) H inf = H max H inf = H max / H inf = H max / H inf = H max / H inf = H max / H inf = H max / - - - - - - - - - - f RD ( Hz ) ω s T RD ( GeV ) H inf = H max H inf = H max / H inf = H max / H inf = H max / H inf = H max / H inf = H max / H inf ( GeV ) ω s f RD = f BBN f RD = f BBN f RD = f BBN f RD = f BBN f RD = f BBN f RD = f BBN H inf ( GeV ) ω s f RD = f BBN f RD = f BBN f RD = f BBN f RD = f BBN f RD = f BBN f RD = f BBN
Figure 4 . Colored regions represent the parameter-space probe-able by LISA (left column) and LIGOO2 run (right column). In the top row we show fixed- H inf slices and in the bottom row fixed- f RD slices. The vertical dashed line indicates f BBN in the top panels, and H max in the bottom panels. Ineach figure, different slices of the parameter space overlay one another, with darker ones placed upperin the stack. Beneath each slice of darker color, there are always hidden lighter colored regions. to the minimum point of the curve, i.e. where the sensitivity is of LISA is maximum (recallthe top and left-bottom panels in Figure 3). Since the LISA sensitivity curve is steep enougharound its minimum, the first intersection point of a GW spectrum lies always very close to thetip of the LISA sensitivity curve (we have checked explicitly that this is always true by varyingthe parameters around their boundary values that separate detection from non-detection).Hence, based on this observation, the detectable-undetectable boundary in the parameterspace is given roughly by the condition h Ω GW ( f LISA ) (cid:38) h Ω (LISA)GW , where f LISA (cid:39) . h Ω (LISA)GW (cid:39) − is the amplitudeof the LISA power law sensitivity curve at such minimum. Using this, we arrive to thefollowing condition for ’detectability’ with LISA in the smooth-transition case, (cid:18) H inf . × GeV (cid:19) (cid:18) . f RD (cid:19) w s − w s+1 (cid:38) . A / ∼ . − . , (4.4)where one should substitute ˜ A s −→ A s for the instant transition case. This rough analyticbound is useful for making quick estimates.– 19 –ISA LIGO O2 ω s Figure 5 . Parameter-space regions probe-able by LISA (left panel) and LIGO O2 (right panel),sampling gradually the value of w s . The color coding in both figures is the same, so it is clear thatLISA can probe some fraction of the parameter space inaccessible to LIGO: in particular the regionwith sufficiently low values of w s , represented by the green-ish colors in the upper left corner of theleft panel, but absent in the right panel. Vertical and horizontal dashed lines represent f BBN and H max , respectively. On the other hand, as it is evident from the panels in Figure 3, whenever a GW spectrumcrosses through the LISA sensitivity curve, it is possible that it will also cross through theLIGO sensitivity curve, depending on the parameters. Recently, the results from cross-correlation analysis on data from Advanced LIGO’s second observing run (O2), combinedwith the results of the first observing run (O1), showed no evidence for the presence of astochastic background in LIGO [57] (the analysis actually focused on a search for powerlaw GW spectra just like in our case). This implies that the parameter space accessible toaLIGO O2 should be subtracted to the parameter space probe-able by LISA that we justquantified above. The parameter space that would have led LIGO O2 to a detection of thisGW background, lie in the following ranges (see right panels in Figure 4)2 . × GeV (cid:46) H inf < . × GeV , . (cid:46) w s < , − Hz (cid:46) f RD < . × − Hz , − GeV (cid:46) T RD < . . As can see there is still margin for a potential detection of this background by LISA, asthe parameter space regions probed by LIGO O2 are in general smaller than those probedby LISA. If we were to obtain the parameter space expected to be probe-able by LIGO O5(using the projected power law sensitivity curve), we would still reach a similar conclusion. Ingeneral, LIGO’s ability to probe the GW background studied here leads to a looser minimum H inf , looser maximum f RD , and tighter minimum w s , than a detection with LISA. Thisis understandable since, as displayed in the top- and left-bottom panels of Figure 3, thesensitivity curves of either LIGO O2 or LIGO O5 are a factor of ∼ O (10 − ) − O (10 − ) lesssensitive than that of LISA, but probe frequencies a factor of ∼ higher than LISA. Sincethe logarithmic slope of our the high-frequency observable branch of our GW background is– 20 –onstrained to be d log ρ GW d log f = w − / w +1 / <
1, it is understandable that there is always a fractionof the parameter space that can lead to a detection in LISA, but not in LIGO (note that thisnever happens in the opposite direction). This is quantified in Figure 5, where we plot withthe same color coding the isocurves w s ( H inf , f RD ) = const , as we vary H inf and f RD , bothfor LISA (left panel) and LIGO O2 (panel). In the figure we can clearly see that the regionwith sufficiently low values of w s , represented by the green-ish colors in the upper left cornerof the LISA panel, is absent in the LIGO panel .The next logical step would be to quantify the ’reduced’ observable parameter space byLISA, after we subtract the parameter space region ruled out by the present non-detectionof a stochastic GW background by LIGO O1+O2. Instead, we will address first in the nextsubsection other upper bounds on stochastic GW backgrounds, due to constraints on theexpansion rate during BBN and CMB decoupling. We will find that such upper boundsare actually stronger (for the problem at hand) than the current upper bounds from LIGOO1+O2 [57]. Gravitational waves contribute to the energy budget of relativistic species in the Universe.Overly abundant GWs may alter too much the expansion rate of the Universe during BBN,which, in turn, may affect the resulting abundances of light elements. To avoid sabotagingthe success of BBN, it is required that [23] (cid:90) f ∗ f BBN h Ω GW ( τ , f ) dff ≤ . × − ∆ N ν (cid:39) . × − , (4.5)where ∆ N ν parametrizes the extra amount of radiation beyond the SM dof [58], wherewe have used the CMB constraint on the number of extra relativistic species ∆ N ν (cid:46) . C.L. [59]. In Eq. (4.5), the lower bound f BBN is the frequency today corresponds to themode crossing the horizon at the onset of BBN, c.f. Eq. (4.2), whereas the upper bound f ∗ isthe frequency today corresponds to the mode re-entering the horizon at the end of inflation(hence it is the high-frequency end of the GW spectrum). The BBN constraint above is anintegral constraint, i.e. a non-local constraint in the frequency space. To get a local, slightlyweaker constraint out of it, we can simply perform the integration in (4.5). As our signalis simply a power-law at high frequencies Ω GW ∝ f − α ) , the integral is dominated by itsupper bound, and we obtain h Ω GW ( f ∗ ) (cid:46) (1 − α s ) · . × − , (4.6)where the lhs should be evaluated using Eq. (3.27) for the instant-transition case, or withEq. (3.31) for the smooth-transition case. Thus, the BBN constraint tends to rule out caseswith a high inflationary scale H inf , large EoS w s , and low transition frequencies f RD (equiv-alently low energy scales T RD ); these parameters would yield an amplitude of h Ω GW ( f ∗ )that is well above the rhs of Eq. (4.6). The BBN bound may therefore rule-out the highestpossible GW signals that the detectors may probe at their respective frequency range. Had we used the LIGO O5 expected sensitivity, the small − w s region probe-able by LISA but inaccessibleto LIGO O5 would still exist, but with a reduced volume of the parameter space. The contribution from extra radiation during any stage of evolution of the Universe is typicallyparametrized in terms of an effective deviation ∆ N ν from the number of SM neutrino species N ν = 3 . – 21 – - - - - - - - - - - f RD ( Hz ) ω s T RD ( GeV ) H inf = H max H inf = H max / H inf = H max / H inf = H max / H inf = H max / H inf = H max / H inf ( GeV ) ω s f RD = f BBN f RD = f BBN f RD = f BBN f RD = f BBN f RD = f BBN f RD = f BBN
Figure 6 . Colored regions represent the parameter space compatible with the requirement of nothaving too much GWs at the onset of BBN (to obtain these contours we have used the smoothtransition GW spectrum Eq. (3.31) evaluated at f ∗ , assuming g ∗ ,RD = g s,RD = 106 .
75 at the SD-to-RD transition). In the left panel we show fixed- H inf slices of the parameter space, whereas in theright panel fixed- f RD slices. In each figure, different slices of the parameter space overlay one another,with darker ones placed lower in the stack. Beneath each slice of darker color, there are always hiddenlighter colored regions. In order to obtain the BBN bound on the parameter space { H inf , w s , f RD } , we need tocompute the ratio of f ∗ to f RD . This will tell us about the maximum value of h Ω GW ( f ) andwhether it exceeds the BBN bound. By construction, we have f ∗ f RD = k ∗ k RD = a ∗ H ∗ a RD H RD . (4.7)To proceed, we express a ∗ /a RD in terms of H ∗ , H RD , and α s using (3.3) and (3.29) inthe smooth transition case. Then, we express H RD in terms of f RD , Ω (0)rad , and H using f RD = k RD / (2 πa ), k RD = a RD H RD , and (3.32) in the smooth-transition case (alternatively(3.23) in the instant-transition case). The result is f ∗ f RD = (cid:16) G RD Ω (0)rad (cid:17) / H H ∗ π f α s , abrupt transition2 (cid:16) − α s1+ α s (cid:17) (cid:16) G RD Ω (0)rad (cid:17) / H H ∗ π f α s , smooth transition (4.8)The relative strength of the high-frequency end of the GW spectra in the two cases is thusΩ (smooth)GW ( τ , f (smooth) ∗ )Ω (abrupt)GW ( τ , f (abrupt) ∗ ) = 2 (cid:16) − α s1+ α s (cid:17) , (4.9)where the ratio ranges from 1 for w s = 1 / / (cid:39) .
59 for w s = 1.For simplicity, from now on we will make the assumption that the Hubble parameterduring inflation H inf is approximately equal to that at the end of inflation H ∗ , i.e. H inf (cid:39) H ∗ .– 22 –sing (4.8), we can obtain the BBN bound (4.6) in terms of H inf , w s , and f RD , and work outthe parameter space region compatible with the BBN bound. The latter is shown in Figure 6for the smooth transition case. If we exclude the region incompatible with the BBN boundfrom the parameter space region probe-able by LISA shown in Figure 4, only a relativelysmall region remains, see left panel of Figure 7. On the other hand, the parameter spaceregion probe-able by current LIGO O2, or even by the will-be LIGO O5 expected sensitivityis completely wiped out when we consider the BBN bound. In other words, there is nopossible value of H inf , w s , or f RD , that can lead to a detection in LIGO without violatingthe BBN bound. This implies that the inflationary background of GWs we are studying hereis expected to be undetectable by current or future runs of LIGO. To put it differently, ifa stochastic signal with power law spectrum was to be discovered by LIGO in the comingyears, it should not be identified with the inflationary GW background in the presence of anSD era.Interestingly, we see in the left panel of Figure 7 that the BBN bound appears to ruleout the region w > .
56, essentially independently of f RD and H inf . To understand this, let usrecall that the parameter space region probe-able by LISA can be described to a good approx-imation by (4.4). In terms of this approximate bound, the cut shown in the Figure 7 is theintersection between the LISA bound saturation surface h Ω GW ( τ , f LISA ; H inf , f RD , w s ) = h Ω (LISA)GW , and the BBN bound saturation surface, h Ω GW ( τ , f ∗ ; H inf , f RD , w s ) = h Ω (BBN)GW .Incidentally, both h Ω GW ( τ , f LISA ) and h Ω GW ( τ , f ∗ ) depend on exactly the same combi-nation of H inf and f RD , namely H /f − α s )RD . Note that this is only true because, as statedearlier, we assumed H inf (cid:39) H ∗ . Consequently, we can combine the two equations in such away that the dependence on H /f − α s )RD drops, leaving behind an equation for α s (or ω s ).This explains why the BBN bound cuts in the left panel in Figure 7 appears as a straightline, independent of H inf and f RD .A similar bound on the amount of GWs can also be obtained from constraints on theHubble rate at CMB decoupling, as this can be used to infer an upper bound on extraradiation components [60–62]. This translates to an upper bound on the amount of GWs,which actually extends to a greater frequency range than the BBN bound, down to f (cid:46) − Hz. In Ref. [60] two cases were identified, depending on the initial conditions ofthe GWs. In the first case, labeled as the adiabatic initial condition, the GW background isassumed to have perturbations imprinted on its energy density following the same distributionas all other components in the Universe. In the second case, labeled as the homogeneousinitial condition, the GW background is not perturbed and the curvature perturbation isthe one of the standard adiabatic case. We view this second option for initial conditionsas more justified, since it applies to most of the known cosmological GW backgrounds, andcertainly to the background studied in this paper. The most recent analysis on this [62]only analyzes the case of adiabatic initial conditions. Extrapolating this result to the caseof GWs with homogeneous initial conditions, Ref. [23] concludes that the constraint shouldyield approximately a bound as h Ω GW ( τ , f ) (cid:46) × − , (4.10)which is a factor ∼ | Probe-able region compatible with CMB - - - - - - - f RD ( Hz ) ω s T RD ( GeV ) H inf = H max H inf = H max H inf = H max H inf = H max H inf = H max H inf = H max - - - - - - - f RD ( Hz ) ω s T RD ( GeV ) H inf = H max H inf = H max H inf = H max H inf = H max Figure 7 . Remaining parameter space region probe-able by LISA after removing the part that isincompatible with the upper bounds on GW stochastic backgrounds. In the left panel we show fixed- H inf slices of the parameter space compatible with the BBN constraint Eq. (4.5), whereas in theright panel we show the analogous plot but imposing the more restrictive CMB constraint Eq. (4.10).If we vary H inf continuously, instead of discretely as done here, the spiky feature of the remainingparameter space region would become infinitely dense, and the resulting figures on the left and rightpanels would have clean, horizontal cuts at w s ≈ .
56 and w s ≈ .
53 respectively, both which areshown above in dashed lines. the parameter space compatible with a detection by LISA and satisfying at the same timeEq. (4.10) is smaller than the parameter space compatible with a detection by LISA whilesatisfying Eq. (4.5) [compare the size of the areas of the parameter regions in the right panelof Figure 7 with respect to those in the left panel]. Identical reasoning as before explains aswell the straight line cut in the right panel in Figure 7 which rules out the region w s (cid:38) . H inf and f RD . The CMB bound is however not as robust as the BBNbound, because Eq. (4.10) is based on a extrapolation of the actual CMB constraint, basedon an adiabatic initial condition for the GW background. We therefore prefer to stick onlythe BBN bound (4.6) in what follows.From the left panel of Figure 7 it becomes clear that the initially large regions ofparameter space compatible with a signal detection by LISA, recall left panels in Figure 4,has shrunk into a small region where the values of the parameters lie in the following ranges6 . × GeV (cid:46) H inf < . × GeV , . (cid:46) w s < . , − Hz (cid:46) f RD (cid:46) . · − Hz , (cid:46) T RD (cid:46)
215 MeV . subsec:BBNandCMBbounds There are other physical effects capable of distorting the shapeof the GW energy spectrum that we have not taken into account in the analysis so far. Theseeffects have not been included in the main analysis because they induce only small deviationsin the results, and depend on some unknown aspects of the physics at high energies. In thissection we study the impact of adding a small red-tilt in the GW spectrum, the spectral– 24 –istortion due to changes in the number of relativistic dof during RD, and an anisotropicstress due to the presence of free-streaming relativistic neutrinos. Due to the fact that the inflationary phase cannot be perfectly de Sitter , we should considerthe effect of a red-tilt in the spectrum of GWs. This is expected e.g. in basic slow-roll infla-tionary scenarios. We can easily incorporate such tilt into our analysis, by simply multiplyingthe previous GW energy spectra with a factorΩ GW ( τ , f ) → Ω GW ( τ , f ) × (cid:18) ff p (cid:19) n t , (4.11)where f p is a pivot frequency (related to CMB scales). The spectral tilt is constrained as − n t (cid:46) .
008 [1, 2], so it changes the tensor spectrum only very gradually with scale. Asmentioned before, in the absence of running of the spectral index, i.e. if dn t d log k = 0, theamplitude of the tensor spectrum only falls by a factor ∼ (10 ) − . ∼ . e ) ∼
26 orders of magnitude separating the CMB scales from the Hubble radius at theend of inflation. This is what justifies our assumption that the low-frequency branch f (cid:28) f RD of the spectrum is scale-invariant. To quantify the effect in our analysis due to adding a tiltin the spectrum, we choose the same pivot scale used in the Planck series of papers, givenby a present day wavenumber k p a (cid:39) .
002 Mpc − , (4.12)which corresponds to a frequency today as f p = 12 π k p a (cid:39) . × − Hz . (4.13)The left panel of Figure 8 shows the remaining parameter space region that is probe-ableby LISA when imposing a red spectral tilt n t = − .
008 in the GW spectrum. Due to theweakening of the GW spectrum, we expect that slightly lower values of H inf , higher values of f RD , and smaller values of w s to be no longer probe-able by LISA. At the same time, we alsoexpect the BBN bound to loosen, allowing higher values of w s to be probe-able. Such shiftingeffects of the parameter space regions are clearly visible in the left panel of Figure 8, wherethe probe-able regions in the presence of tilt (colored areas, delimited by solid lines) aredisplaced upwards (i.e. to higher values of w s ) and towards the left (i.e. to smaller values of f RD ), when compared to the probe-able regions in the absence of tilt (empty areas, delimitedby dashed lines). Overall, the probe-able regions still represent very small islands of theparameter space. As expected, the presence of a tilt (even in the maximum case as we chose)results in only a small correction. dof during RD While the energy density of the GWs always scales as a − , the total energy of the universeduring the RD period at temperature Tρ tot = π g ∗ T , (4.14)– 25 – - - - - - - - f RD ( Hz ) ω s T RD ( GeV ) H inf = H max H inf = H max H inf = H max H inf = H max H inf = H max H inf = H max - - - - - - - f RD ( Hz ) ω s T RD ( GeV ) H inf = H max H inf = H max H inf = H max H inf = H max H inf = H max H inf = H max Figure 8 . The remaining parameter space probe-able by LISA, accounting for the presence of ared spectral tilt n t = − .
008 (left panel), and accounting for the effect of changes in the number ofrelativistic dof (right panel). In both panels, we show for comparison as empty regions delimited bydashed line contours the probe-able regions before we considered the presence of a spectral tilt orcorrected by the pre-factor C ∆ g ∗ . may evolve in more complicated ways due to possible changes in the effective number g ∗ of relativistic degrees of freedom ( dof ). In general, the temperature evolves according toentropy conservation s = 2 π g s T ∝ a − ⇒ T ∝ g − / s a − , (4.15)with g s the number of entropic dof . Thus, the radiation energy density evolves as ρ tot ∝ g ∗ g − / s a − ∼ g − / ∗ a − , (4.16)where in the last equality we have taken into account the fact that energy and entropy dof sare approximately equal. If there is a sudden change in the number of dof , say, from g ∗ to g ∗ , ρ tot would change by a factor of ∼ ( g ∗ /g ∗ ) / , and consequently Ω GW ≡ ( dρ GW /d ln f ) /ρ tot will change by a factor of ∼ ( g ∗ /g ∗ ) / . Note that the change in Ω GW only affects themodes that are already sub-horizon before the change g ∗ −→ g ∗ takes place; those that aresuper-horizon remain frozen. After a series of changes in the relativistic dof , the net changein Ω GW for a particular mode k depends only on the value at the moment of horizon re-entry g ∗ ,k ≡ g ∗ ( k = aH ). As a rule of thumb, modes that re-enter the horizon earlier get a largersuppression, since the temperature and hence g ∗ ,k were larger then. The dependence on thenumber of relativistic dof is encoded in the pre-factor G k [c.f. Eq. 2.12] in the amplitude ofthe GW spectrum [recall Eq. (2.11)], which we rewrite here for convenience,Ω GW ( τ , k ) ∝ G k , G k ≡ (cid:18) g ∗ ,k g ∗ , (cid:19) (cid:18) g s, g s,k (cid:19) . (4.17)As far as the Standard Model is concerned , there are only two events that lead to sig-nificant changes (in terms of percentage) in the relativistic dof that we need to pay attentionto. They are the QCD phase transition at around the temperature T QCD ∼
200 MeV, andthe electron-positron annihilation at a temperature T e + e − = 0 . Changes in beyond the SM dof can be also probed by LISA [63], but here we focus on the SM dof only. – 26 – ∗ ,k = g s,k = 106 .
75, which is valid whenever all the SM particles are relativistic and dominatethe energy budget of the Universe. We evaluated the plateau amplitude Eq. (2.11) using suchvalues for the number of dof , which lead to a suppression G k (cid:39) .
39. After the QCD phasetransition, however, the number of dof drop to g ∗ ,k = g s,k = 10 .
75, so that the suppressionis smaller, with G k ∼ .
83. After electron-positron annihilation and neutrino decoupling, thenumber of dof drop to their final values g ∗ ,k = g ∗ , = 3 .
36 and g s,k = g s, (cid:39) .
91, so thatthere is no suppression any more as G k = 1 (the neutrinos are treated as if they were masslessbecause they are relativistic during the radiation era, and give a negligible contribution to thetotal energy density during the matter era). Therefore, to account for the effect of changesin the relativistic dof in the RD epoch, we simply need to modify the main result as followsΩ GW ( τ , k ) → C ∆ g ∗ ( k )Ω GW ( τ , k ) , C ∆ g ∗ ( k ) ≡ G k ( g ∗ ,k , g s,k ) G k (106 . , . . (4.18)A simple approximation for the correction C ∆ g ∗ ( k ) is to consider it as a piece-wise constantfunction, discontinuous at every moment the number of dof changes significantly, i.e. C ∆ g ∗ ( k ) = . , k > k QCD . , k QCD > k > k e + e − . , k < k e + e − , (4.19)with k QCD , k e + e − the comoving horizon scales at the QCD phase transition and electron-position annihilation, respectively. Here we have implicitly assumed that the QCD phasetransition happened during the RD epoch, i.e. T RD > T QCD (equivalently k RD > k QCD ). Inthe bottom-left panel of Figure 3, we show an example of a GW spectrum (red dashed line)taking into account the correction (4.19). For comparison we show in the same panel (graysolid line) the GW spectrum for the same parameters but without the correction C ∆ g ∗ .If the change in g ∗ happens during the SD epoch, the correction needs however to bemodified. During the SD period, a jump in the value of g ∗ in the (sub-dominant) radiationsector does not entail a jump in Ω GW . Therefore, if the QCD phase transition happenedduring the SD epoch , i.e. if T RD < T QCD , the correction factor is then changed to C ∆ g ∗ ( k ) = (cid:40) .
15 for k > k e + e − .
59 for k < k e + e − , (4.20)which reflects that in this case, only changes in the number of dof taking place during theRD, impact on the final form of Ω GW ( τ , k ). In the bottom-left panel of Figure 3, we showan example of a GW spectrum (blue solid line) taking into account the correction (4.20),and for comparison we also show the GW spectrum for the same parameters but without thecorrection (gray solid line). The pre-factors C ∆ g ∗ (cid:38) g RD ∼ .
75. In these cases, the overall amplitudeof the GW energy spectra for the frequency range of interest ( f (cid:29) f e + e − ∼ − Hz) is inreality a factor C ∆ g ∗ ∼ We always assume that electron-positron annihilation occurs in the RD epoch, as the minimum reheatingtemperature required for a consistent cosmology is T RH (cid:38) T BBN ∼ MeV. – 27 – ∆ g ∗ ( k ) factor taken into account. As the probe-able temperature scales correspond all tosmall values T (cid:46) T QCD , it is important that we apply the prescription given by Eq. (4.20),and not by (4.19). This represents therefore a relatively important correction to the GWspectra with such low f RD frequencies, as we were extra-suppressing their amplitude (inthe frequency range of interest) by a factor ∼ .
15. The impact of this in our parameterconstraint analysis, is that the remaining probe-able regions by LISA are displaced to theright (towards larger values of f RD ) and are slightly wider in w s , see colored regions in theright panel of Figure 8, which now extend the probe-able region down to w s (cid:39) .
46. Theprobe-able parameter space ranges change now to10 GeV (cid:46) H inf ≤ . × GeV , . (cid:46) w s < . , − Hz (cid:46) f RD < . × − Hz , (cid:46) T RD <
150 MeV . Overall, this effect is still a small correction, which does not change the fact that aftersubtraction of the parameter space region incompatible with the BBN bound, the remainingparameter space probe-able by LISA is still a very small island of the full parameter space.Finally we note that correcting our previous GW spectra by a factor ∼ .
15 still does notchange the fact that there is no surviving parameter space compatible with having a signaldetection by LIGO O2 or O5, without violating the BBN constraint.
Free-streaming decoupled relativistic particles moving along their geodesics back-react to thespacetime in an anisotropic way, which manifests itself as an anisotropic contribution to thestress-energy tensor [10, 64]. An extra term appears then in the right hand side of the equationof motion of tensor perturbations (2.4). If the free-streaming species in question are stable attime scales much longer than the instantaneous Hubble time, their energy density ρ FS makeup an approximately constant fraction f FS of the total energy density of the universe duringRD. In that case, their anisotropic contribution to the stress energy tensor can be calculatedexplicitly [10] Π k = − ρ FS ( τ ) (cid:90) ττ ν dec dτ (cid:48) (cid:26) j [ k ( τ − τ (cid:48) )] k ( τ − τ (cid:48) ) (cid:27) h (cid:48) k ( τ (cid:48) ) , (4.21)so that (2.4) becomes an integro-differential equation h (cid:48)(cid:48) k + 2 a (cid:48) a h (cid:48) k + k h k = − FS (cid:18) a (cid:48) a (cid:19) (cid:90) ττ ν dec dτ (cid:48) (cid:26) j [ k ( τ − τ (cid:48) )] k ( τ − τ (cid:48) ) (cid:27) h (cid:48) k ( τ (cid:48) ) . (4.22)The fact that the kernel on the right hand side can be written as a sum of spherical Besselfunctions suggests that there are solutions that can be expressed as a sum of spherical Besselfunctions. Keeping the first five non-vanishing terms of such sum provides a solution withan error of less than 0 . | A | <
1, relativeto that in the absence of free-streaming particles. The damping factor | A | ranges from 1 to0.35, corresponding to Ω FS = 0 and Ω FS = 1, respectively.This damping effect is however local in frequency space. In the case of the SM neutrinos,for example, the damping only applies to modes with the present frequency range f eq < – 28 – < f ν dec , i.e. those that crossed the horizon after neutrino decoupling at T (cid:46) In this paper we have studied how a period characterized by a stiff equation of state 1 / 1, spanning from the end of inflation to the onset of RD, impacts on the GW backgroundfrom inflation. Due to this SD period, the GW energy density spectrum acquires a blue tilt n t ≡ d log ρ GW d log f = 2 ( w s − / w s +1 / in the frequency regime f (cid:29) f RD associated to the modes thatre-enter the horizon during the stiff epoch (here f RD is the frequency today corresponding tothe horizon scale at the time of the SD-to-RD transition). As a result, the energy spectrumin the high frequency domain f (cid:29) f RD is considerably amplified relative to the (quasi-)scaleinvariant part in the low frequency domain f (cid:28) f RD , corresponding to the modes re-enteringthe horizon during RD.To obtain an exact expression for the spectrum around the scales reentering during theSD-to-RD transition, we have done a matching in the frequency space between the (highfrequency) modes that reenter during SD, and the (low frequency) modes that reenter duringRD. This requires the use of the exact time-dependence of the scale factor, Hubble rate andequation of state of the Universe, and a careful treatment of the evolution of the modescrossing the horizon along the transition itself. We have obtained the exact transfer functionof the GW energy density spectrum around the frequencies f ∼ f RD , both numerically and(whenever possible) analytically, considering both ’instant’ and smooth modelings of theSD-to-RD transition, see Eqs. (3.21) and (3.36). We find that the high frequency branchof the spectrum in the smooth case is a factor 2 − α s larger than in the instant case, where α s = 2 / (1 + 3 w s ), recall Eq. (3.31). We consider a smooth transition as more realistic, sincein this case all background quantities (the scale factor, the Hubble rate and the equation ofstate) change smoothly and continuously along the SD-to-RD transition. Because of this wedecided to perform our parameter analysis only using the GW spectra obtained for a smoothtransition case. Using the GW spectrum from the instant-transition case, which assumesa sudden jump in the equation of state, leads in any case to only marginal changes in theanalysis results.The shape of the GW spectrum is controlled by w s , f RD , and the energy scale of inflation H inf . We have determined the parameter space compatible with a detection of this signal byAdvanced LIGO and LISA. See Figures 4 and 5, which exhibit that, a priori , a large region ofthe parameter space is potentially observable. Consistency with upper bounds on stochasticGW backgrounds, due to constraints on the expansion rate during BBN and CMB decoupling,rules out however a significant fraction of the would-be observable parameter space. We findthat the signal becomes completely inaccessible to Advanced LIGO, independently of theparameters. In other words, there is no solution in the parameter space { w s , f RD , H inf } compatible with a detection at LIGO, while not violating at the same time the BBN/CMBconstraints [c.f. Eqs. (4.5), (4.10)]. This is independent of whether we consider the current runO2 or the projected run O5. In the case of LISA, a small region of parameter space remainsstill probe-able. This is depicted in the plots of Figure 7, which show clearly that the initially– 29 –arge portion of parameter space compatible with a signal detection by LISA, shrinks intoa very small region. This region corresponds to scenarios with a large inflationary scale H inf (cid:38) . H max , transitioning into a RD stage at a low temperature 1 MeV (cid:46) T RD (cid:46) − Hz (cid:46) f RD (cid:46) . · − Hz), and where the stiff EoS is confined within a narrow rangeof small values, 0 . (cid:46) w s (cid:46) . 56 [c.f. the end of Section 4.2].We have also studied the impact of adding a small red tilt into the GW spectrum. Wefind that slightly lower values of H inf , higher values for f RD , and small values of w s , are nolonger probe-able by LISA, see left panel of Figure 8. Overall, the probe-able regions stillrepresent a small island of the parameter space. The inclusion of a constant tilt (even in themaximum case n t = − . 008 allowed by CMB observations) represents therefore only a verysmall correction. We have considered as well the spectral correction in GW spectrum dueto the inclusion of changes in the number of relativistic dof during RD. If the temperatureof the radiation fluid component at the time of the SD-to-RD transition is large enough, T RD > T QCD (cid:39) 200 MeV, this does not affect the parameter constraint analysis (as the GWsignal does not change within the observable frequency range). However, for low reheatingtemperatures T RD (cid:46) T QCD (as actually required by the surviving probe-able region of param-eter space, for which T RD (cid:46) 150 MeV), this effect yields a correction in the GW spectrumamplitude of a factor ∼ f RD , and become slightly wider in the range of w s , see the rightpanel in Figure 8. The probe-able region now extends down to w s (cid:39) . 46. Overall, this effectis still a small correction, which does not change the fact that the parameter space regionprobe-able by LISA still represents a very small island of the full parameter space. We notethat neither of these corrections, adding a small tilt or changing the number of dof , changethe fact that there is no surviving parameter space compatible with having a signal detectedby advanced LIGO O2/O5 while not violating the BBN constraint. Finally we note that theexpected anisotropic stress due to the presence of the free-streaming relativistic SM neutrinosdoes not affect the GW signal in the frequency range that is relevant to our analysis.A stiff epoch is a crucial aspect in various early Universe scenarios, for instance the Quintessential inflation [25–33], where inflation is followed by an epoch dominated by thekinetic energy of the inflaton, with the potential adjusted to describe the observed darkenergy today. Various mechanisms for excitation of the radiation component that eventuallywill dominate the energy budget in these scenarios have been proposed and worked out indetail [65–69] such that the duration of the SD phase is controlled by the relevant parametersinvolved in each construction. One of such mechanisms is the original gravitational reheating model [34, 35], where the Universe is reheated by the decay products of inflationary spectatorfields excited during or towards the end of inflation. Basic implementations of this idea havebeen shown however to be inconsistent with BBN/CMB constraints [36], unless unnatural ad hoc constructions with many fields and identical tuned properties are accepted. Variantsof the idea that avoid the previous inconsistency have been proposed, like the Higgs-reheating scenario [37], where the Standard Model (SM) Higgs is a spectator field with a non-minimalcoupling to gravity, and the Universe is reheated into SM relativistic species (decay productsof the Higgs) during the SD period. The same mechanism can actually be realized withgeneric self-interacting scalar fields other than the SM Higgs [38], see [39] for a recent re-analysis of the idea. In general, models with blue-tilted inflationary GW spectrum due to thepresence of an SD epoch have been studied in various contexts in the past, see e.g. [7–9, 11, 17–– 30 –1, 23, 40, 41]. Unfortunately, the small island of parameter space compatible with a detectionby LISA, H inf (cid:38) . H max , 1 MeV (cid:46) T RD (cid:46) 150 MeV, and 0 . (cid:46) w (cid:46) . 56, does notseem particularly appealing/suitable from the model building perspective: these parameterscorrespond to scenarios where the SD epoch spans for many decades in energy scale, from theend of inflation till the onset of RD. From the point of view of field theoretical constructionsit seems actually difficult to obtain an equation of state within the narrow allowed range.For whichever scenario resting upon an early SD period after inflation, BBN sets stringentconstraints (recall Figures 6 and 7): if f RD (cid:46) − Hz (typically corresponding to a long SDperiod), the only observable window by LISA is the small island of parameter space quantifiedjust above; if f RD (cid:38) − Hz, the GW signal is then simply not observable. Whicheverthe scenario studied, this simple frequency domain rule must be taken into considerationwhen assessing the observability of the inflationary GW background distorted because of thepresence of a period with stiff equation of state. In light of our analysis, we find it veryunlikely that this GW signal will be detected in the future.Let us note that in our analysis we have asumed Eq. (2.7) for the spectral tilt of thetensor modes excited during inflation. This is a relation expected whenever inflation iswell described by a quasi-de Sitter background, as we assume in this work. In alternativemodels, n t could be very different. However, except for very concrete inflationary set-up’s,the inflationary tensor spectrum rarely exhibits a large blue tilt for the modes excited duringinflation. In our analysis we want to constrain the expansion rate of the Universe afterinflation, so we avoid on purpose introducing n t as a free parameter, as this would representa degenerate constraint between the post-inflationary expansion rate and the freedom tochose n t during inflation. We rather focused on the ability of detectors to probe the post-inflationary expansion rate alone, assuming a natural implementation of inflation. Becauseof this, we have not considered standard constraints due to pulsar timing arrays (PTA) in ouranalaysis: PTA cannot set constraints in our case, as BBN imposes a minimum SD-to-RDfrequency of the order of f RD = f BBN ∼ − Hz. This means that even in the highestpossible slope of the tensor spectrum, when h Ω GW ∝ f for w (cid:39) 1, the tensor backgroundcan only be at most h Ω GW = h Ω (0) GW | (plateau) × ˜ A s × ( f PTA /f BBN ) ∼ − at the typicalfrequencies of PTA experiments f PTA ∼ − Hz. Such large tilt are however forbiddendue to the BBN constraint, as shown in our analysis, c.f. Fig. 7. Therefore, at the PTAfrequencies, taking into accout the BBN constraint for the lowest values of f RD , the valueof GW background at f PTA is always below any of the projected sensitivities expected forpresent and future PTA experiments, in particular h Ω GW ( f PTA ) (cid:46) − . For constraintson GW signals assuming that n t can be a free variable, see e.g. [42–52].As a final remark, let us mention that if we consider an inflationary tensor tilt withrunning, or a smooth transition from inflation into the SD epoch, there are possible waysto reduce the tension with the BBN/CMB bounds, and hence to enlarge the probe-ableparameter space of the GW signal. In realistic inflationary models, typically the deviationfrom slow-roll becomes more noticeable towards the end of inflation, naturally leading to arunning in the tensor tilt. This will reduce further the amplitude of the GW modes in thehigh frequency end of the spectrum, but this becomes a model dependent computation. Forstandard single field inflation monomial potentials, one may obtain a reduction factor of thespectral energy amplitude at the highest frequency mode of the order of (cid:46) . 1. Furthermore,the transition into an SD era is also a model dependent mechanism, and if the transitionwere smooth, say lasting for a few Hubble times since the end of inflation, this would alsoreduce the amplitude at the high frequency modes of the spectrum (as the earliest modes– 31 –eentering the horizon would not evolve in an SD background yet, or the initial stiffness ofthe EoS would be softer, i.e. closer to 1 / Acknowledgements . We thank Benjamin Stefanek for commments on the manuscript.The work of DGF was supported partially by the ERC-AdG-2015 grant 694896 and partiallyby the Swiss National Science Foundation (SNSF). References [1] Planck collaboration, Y. Akrami et al., Planck 2018 results. X. Constraints on inflation , .[2] BICEP2, Keck Array collaboration, P. A. R. Ade et al., BICEP2 / Keck Array x:Constraints on Primordial Gravitational Waves using Planck, WMAP, and New BICEP2/KeckObservations through the 2015 Season , Phys. Rev. Lett. (2018) 221301, [ ].[3] L. P. Grishchuk, Amplification of gravitational waves in an istropic universe , Sov. Phys. JETP (1975) 409–415.[4] A. A. Starobinsky, Spectrum of relict gravitational radiation and the early state of the universe , JETP Lett. (1979) 682–685.[5] V. A. Rubakov, M. V. Sazhin and A. V. Veryaskin, Graviton Creation in the InflationaryUniverse and the Grand Unification Scale , Phys. Lett. (1982) 189–192.[6] R. Fabbri and M. d. Pollock, The Effect of Primordially Produced Gravitons upon theAnisotropy of the Cosmological Microwave Background Radiation , Phys. Lett. (1983)445–448.[7] M. Giovannini, Gravitational waves constraints on postinflationary phases stiffer thanradiation , Phys. Rev. D58 (1998) 083504, [ hep-ph/9806329 ].[8] M. Giovannini, Production and detection of relic gravitons in quintessential inflationary models , Phys. Rev. D60 (1999) 123511, [ astro-ph/9903004 ].[9] L. A. Boyle and P. J. Steinhardt, Probing the early universe with inflationary gravitationalwaves , Phys. Rev. D77 (2008) 063504, [ astro-ph/0512014 ].[10] Y. Watanabe and E. Komatsu, Improved Calculation of the Primordial Gravitational WaveSpectrum in the Standard Model , Phys. Rev. D73 (2006) 123515, [ astro-ph/0604176 ].[11] L. A. Boyle and A. Buonanno, Relating gravitational wave constraints from primordialnucleosynthesis, pulsar timing, laser interferometers, and the CMB: Implications for the earlyUniverse , Phys. Rev. D78 (2008) 043531, [ ].[12] S. Kuroyanagi, T. Chiba and N. Sugiyama, Precision calculations of the gravitational wavebackground spectrum from inflation , Phys. Rev. D79 (2009) 103501, [ ].[13] S. Kuroyanagi, T. Chiba and N. Sugiyama, Prospects for Direct Detection of InflationaryGravitational Waves by Next Generation Interferometric Detectors , Phys. Rev. D83 (2011)043514, [ ].[14] S. Kuroyanagi, K. Nakayama and J. Yokoyama, Prospects of determination of reheatingtemperature after inflation by DECIGO , PTEP (2015) 013E02, [ ].[15] M. S. Turner, Coherent Scalar Field Oscillations in an Expanding Universe , Phys. Rev. D28 (1983) 1243. – 32 – 16] K. D. Lozanov and M. A. Amin, Equation of State and Duration to Radiation Dominationafter Inflation , Phys. Rev. Lett. (2017) 061301, [ ].[17] A. Riazuelo and J.-P. Uzan, Quintessence and gravitational waves , Phys. Rev. D62 (2000)083506, [ astro-ph/0004156 ].[18] V. Sahni, M. Sami and T. Souradeep, Relic gravity waves from brane world inflation , Phys.Rev. D65 (2002) 023518, [ gr-qc/0105121 ].[19] H. Tashiro, T. Chiba and M. Sasaki, Reheating after quintessential inflation and gravitationalwaves , Class. Quant. Grav. (2004) 1761–1772, [ gr-qc/0307068 ].[20] M. Giovannini, Stochastic backgrounds of relic gravitons, T Λ CDM paradigm and the stiff ages , Phys. Lett. B668 (2008) 44–50, [ ].[21] M. Giovannini, Thermal history of the plasma and high-frequency gravitons , Class. Quant.Grav. (2009) 045004, [ ].[22] Y. Cui, M. Lewicki, D. E. Morrissey and J. D. Wells, Cosmic Archaeology with GravitationalWaves from Cosmic Strings , Phys. Rev. D97 (2018) 123505, [ ].[23] C. Caprini and D. G. Figueroa, Cosmological Backgrounds of Gravitational Waves , Class.Quant. Grav. (2018) 163001, [ ].[24] Y. Cui, M. Lewicki, D. E. Morrissey and J. D. Wells, Probing the pre-BBN universe withgravitational waves from cosmic strings , JHEP (2019) 081, [ ].[25] P. J. E. Peebles and A. Vilenkin, Quintessential inflation , Phys. Rev. D59 (1999) 063505,[ astro-ph/9810509 ].[26] M. Peloso and F. Rosati, On the construction of quintessential inflation models , JHEP (1999) 026, [ hep-ph/9908271 ].[27] G. Huey and J. E. Lidsey, Inflation, brane worlds and quintessence , Phys. Lett. B514 (2001)217–225, [ astro-ph/0104006 ].[28] A. S. Majumdar, From brane assisted inflation to quintessence through a single scalar field , Phys. Rev. D64 (2001) 083503, [ astro-ph/0105518 ].[29] K. Dimopoulos and J. W. F. Valle, Modeling quintessential inflation , Astropart. Phys. (2002) 287–306, [ astro-ph/0111417 ].[30] C. Wetterich, Variable gravity Universe , Phys. Rev. D89 (2014) 024005, [ ].[31] C. Wetterich, Inflation, quintessence, and the origin of mass , Nucl. Phys. B897 (2015)111–178, [ ].[32] M. W. Hossain, R. Myrzakulov, M. Sami and E. N. Saridakis, Variable gravity: A suitableframework for quintessential inflation , Phys. Rev. D90 (2014) 023512, [ ].[33] J. Rubio and C. Wetterich, Emergent scale symmetry: Connecting inflation and dark energy , Phys. Rev. D96 (2017) 063509, [ ].[34] L. H. Ford, Gravitational Particle Creation and Inflation , Phys. Rev. D35 (1987) 2955.[35] B. Spokoiny, Deflationary universe scenario , Phys. Lett. B315 (1993) 40–45, [ gr-qc/9306008 ].[36] D. G. Figueroa and E. H. Tanin, Inconsistency of an inflationary sector coupled onlygravitationally , .[37] D. G. Figueroa and C. T. Byrnes, The Standard Model Higgs as the origin of the hot Big Bang , Phys. Lett. B767 (2017) 272–277, [ ].[38] K. Dimopoulos and T. Markkanen, Non-minimal gravitational reheating during kination , JCAP (2018) 021, [ ].[39] T. Opferkuch, P. Schwaller and B. A. Stefanek, Ricci Reheating , . – 33 – 40] M. Giovannini, Stochastic backgrounds of relic gravitons: a theoretical appraisal , PMC Phys. A4 (2010) 1, [ ].[41] B. Li, P. R. Shapiro and T. Rindler-Daller, Bose-Einstein-condensed scalar field dark matterand the gravitational wave background from inflation: new cosmological constraints and itsdetectability by LIGO , Phys. Rev. D96 (2017) 063505, [ ].[42] W. Zhao, Constraint on the early Universe by relic gravitational waves: From pulsar timingobservations , Phys. Rev. D83 (2011) 104021, [ ].[43] W. Zhao, Y. Zhang, X.-P. You and Z.-H. Zhu, Constraints of relic gravitational waves by pulsartiming arrays: Forecasts for the FAST and SKA projects , Phys. Rev. D87 (2013) 124012,[ ].[44] S. Kuroyanagi, T. Takahashi and S. Yokoyama, Blue-tilted Tensor Spectrum and ThermalHistory of the Universe , JCAP (2015) 003, [ ].[45] R. Jinno, T. Moroi and T. Takahashi, Studying Inflation with Future Space-Based GravitationalWave Detectors , JCAP (2014) 006, [ ].[46] L. Lentati et al., European Pulsar Timing Array Limits On An Isotropic StochasticGravitational-Wave Background , Mon. Not. Roy. Astron. Soc. (2015) 2576–2598,[ ].[47] P. D. Lasky et al., Gravitational-wave cosmology across 29 decades in frequency , Phys. Rev. X6 (2016) 011035, [ ].[48] NANOGrav collaboration, Z. Arzoumanian et al., The NANOGrav Nine-year Data Set:Limits on the Isotropic Stochastic Gravitational Wave Background , Astrophys. J. (2016)13, [ ].[49] X.-J. Liu, W. Zhao, Y. Zhang and Z.-H. Zhu, Detecting Relic Gravitational Waves by PulsarTiming Arrays: Effects of Cosmic Phase Transitions and Relativistic Free-Streaming Gases , Phys. Rev. D93 (2016) 024031, [ ].[50] S. Kuroyanagi, T. Chiba and T. Takahashi, Probing the Universe through the StochasticGravitational Wave Background , JCAP (2018) 038, [ ].[51] F. D’Eramo and K. Schmitz, Imprint of a scalar era on the primordial spectrum ofgravitational waves , .[52] N. Bernal and F. Hajkarim, Primordial Gravitational Waves in Non-standard Cosmologies , .[53] K. D. Lozanov and M. A. Amin, Self-resonance after inflation: oscillons, transients andradiation domination , Phys. Rev. D97 (2018) 023533, [ ].[54] Y. Zhang, Y. Yuan, W. Zhao and Y.-T. Chen, Relic gravitational waves in the acceleratingUniverse , Class. Quant. Grav. (2005) 1383–1394, [ astro-ph/0501329 ].[55] M. S. Turner, M. J. White and J. E. Lidsey, Tensor perturbations in inflationary models as aprobe of cosmology , Phys. Rev. D48 (1993) 4613–4622, [ astro-ph/9306029 ].[56] S. Kuroyanagi, C. Gordon, J. Silk and N. Sugiyama, Forecast Constraints on Inflation fromCombined CMB and Gravitational Wave Direct Detection Experiments , Phys. Rev. D81 (2010)083524, [ ].[57] LIGO Scientific, Virgo collaboration, B. P. Abbott et al., A search for the isotropicstochastic background using data from Advanced LIGO’s second observing run , .[58] G. Mangano, G. Miele, S. Pastor, T. Pinto, O. Pisanti and P. D. Serpico, Relic neutrinodecoupling including flavor oscillations , Nucl. Phys. B729 (2005) 221–234, [ hep-ph/0506164 ].[59] R. H. Cyburt, B. D. Fields, K. A. Olive and T.-H. Yeh, Big Bang Nucleosynthesis: 2015 , Rev.Mod. Phys. (2016) 015004, [ ]. – 34 – 60] K. M. Smith, W. Hu and M. Kaplinghat, Cosmological Information from Lensed CMB PowerSpectra , Phys. Rev. D74 (2006) 123002, [ astro-ph/0607315 ].[61] I. Sendra and T. L. Smith, Improved limits on short-wavelength gravitational waves from thecosmic microwave background , Phys. Rev. D85 (2012) 123002, [ ].[62] L. Pagano, L. Salvati and A. Melchiorri, New constraints on primordial gravitational wavesfrom Planck 2015 , Phys. Lett. B760 (2016) 823–825, [ ].[63] R. R. Caldwell, T. L. Smith and D. G. E. Walker, Using a Primordial Gravitational WaveBackground to Illuminate New Physics , .[64] S. Weinberg, Damping of tensor modes in cosmology , Phys. Rev. D69 (2004) 023503,[ astro-ph/0306304 ].[65] G. N. Felder, L. Kofman and A. D. Linde, Inflation and preheating in NO models , Phys. Rev. D60 (1999) 103505, [ hep-ph/9903350 ].[66] J. De Haro and L. Arest Sal, Reheating constraints in quintessential inflation , Phys. Rev. D95 (2017) 123501, [ ].[67] J. de Haro and L. Arest Sal, Reheating via Gravitational Particle Production in Simple Modelsof Quintessence or CDM Inflation , Galaxies (2017) 78, [ ].[68] J. Haro, W. Yang and S. Pan, Reheating in quintessential inflation via gravitational productionof heavy massive particles: A detailed analysis , JCAP (2019) 023, [ ].[69] J. de Haro, S. Pan and L. Arest Sal, Understanding gravitational particle production inquintessential inflation , ..