{'fraction_non_alphanumeric': 0.041773635864484245, 'fraction_numerical': 0.012257557095582234, 'mean_word_length': 4.505105792010457, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Precision cosmology and the stiff-amplified gravitational-wave background from inflation: NANOGrav, Advanced LIGO-Virgo and the Hubble tension

The recent NANOGrav finding of a common-spectrum process has invited interpretations as possible evidence of a primordial stochastic gravitational-wave background (SGWB) stronger than predicted by standard inflation+LCDM. Such an SGWB would contribute an extra radiation component to the background Universe which may affect its expansion history. As such, it may help alleviate the current Hubble tension, a novel connection between gravitational waves and cosmology. We demonstrate this by considering a cosmological model, the"standard inflation + stiff amplification"scenario, with two components added to the LCDM model: a stiff component (w=1) and the primordial SGWB. Previously, we showed that even for standard inflation, the SGWB may be detectable at the high frequencies probed by laser interferometers, if it is amplified by a possible early stiff era after reheating. Models that boost the SGWB enough to cause significant backreaction, however, must still preserve the well-measured radiation-matter equality, as precision cosmology demands. For that, we calculate the fully-coupled evolution of the SGWB and expansion history, sampling parameter space (tensor-to-scalar ratio, reheating temperature and temperature at stiff-to-radiation equality). We then perform a joint analysis of the NANOGrav results and latest upper bounds from Planck, big bang nucleosynthesis and Advanced LIGO-Virgo, to constrain the model. The resulting blue-tilted, stiff-amplified SGWB is still too small to explain the NANOGrav results. However, if someday, Advanced LIGO-Virgo detects the SGWB, our model can explain it within standard inflation (without requiring an initial blue tilt). Meanwhile, this model may bring current high-z measurements of the Hubble constant within 3.4 sigma of the low-z measurements by SH0ES (from 4.4 sigma) and within 2.6 sigma of those by H0LiCOW (from 3.1 sigma), reducing the tension.


Introduction
A stochastic gravitational-wave background (SGWB) from primordial tensor fluctuations is generically produced in the inflationary paradigm [1][2][3]. Once deemed too small to be detected, this primordial SGWB is now possibly within reach of various experiments, from the cosmic microwave background (CMB) to gravitational-wave (GW) interferometers, over a wide range of frequencies [4,5]. It may even contribute significantly enough to the energy content of the Universe as to affect the expansion history, with possible observable consequences beyond its direct detection [e.g., 6,7]. These direct and indirect probes of the primordial SGWB can, therefore, potentially reveal important information about inflation and other physical processes in the early Universe, which are otherwise poorly understood [e.g., 8,9].
In fact, even before inflation was proposed, Grishchuk realized that in an expanding Universe, significant parametric amplification can occur, not only for classical gravitational waves (GWs), but even for quantum fluctuations of the vacuum [10,11]. It requires that (1) modes spend time outside the Hubble radius (i.e., the background Universe expands more rapidly than GWs vary in time), when (2) the Universe is not radiation-dominated (RD). When both conditions are met, GWs, or tensor fluctuations, will be amplified relative to the "adiabatic" solution (for which h ∝ 1/a for modes always well-inside the Hubble radius).
The inflationary paradigm [12][13][14] naturally provides such a period that enables parametric amplification and production of macroscopic GWs from initial vacuum fluctuations. When tensor modes are stretched well outside the Hubble radius, their amplitudes become time-independent, or "frozen" [1,15]. These amplitudes define the primordial tensor power spectrum. For standard single-field, slow-roll inflation, their distribution is nearly Gaussian, with a nearly scale-invariant power spectrum that satisfies the consistency relation [16,17]. After inflation ends, tensor modes start to reenter the Hubble radius and each, thereafter, evolves according to the adiabatic solution, redshifting like radiation. Together, all modes that reentered and remain inside the Hubble radius constitute the primordial SGWB.
The parametric amplification regime for a given mode spans its Hubble exit and reentry [18]. 1 While all modes of interest exit during inflation, different modes can reenter during post-inflationary eras with different equations of state (EoS). This actually leads to another kind of amplification/attenuation of the primordial SGWB, as we describe below, in terms of the departure of the amplitudes of modes at a given time after their reentries, relative to those if the EoS of the Universe during their reentries were radiation-like (w ≡p/ρ = 1/3).
Our observed Universe must undergo a standard RD era which begins no later than big bang nucleosynthesis (BBN) and ends at radiation-matter equality. For nearly scale-invariant initial conditions, the contribution to the present-day SGWB energy spectrum, Ω GW (f ) ≡ d Ω GW /d ln f , by modes that reentered during this RD era is nearly frequency-independent. This results in a spectrum with a long "plateau" [19]. In what follows, we shall henceforth use the term amplification to refer, not to the parametric amplification effect described above, but rather to the amplification of a mode at a given time after it reenters, relative to this plateau associated with Hubble reentries that take place during the RD era. For modes of longer wavelengths that reenter during the matter-dominated (MD) era (w = 0) which follows the RD era, Ω GW (f ) is amplified relative to the plateau, since the time-dependence of the Hubble parameter then differs from that of an RD Universe [20,21]. On the other hand, for modes with short enough wavelength to reenter before BBN, the possibility exists for amplification, too, since the expansion history or, equivalently, the EoS of the Universe during this period are poorly constrained and may also depart from w = 1/3. In fact, for these modes of higher frequencies, Giovannini [22] considered the interesting case in which Ω GW (f ) is amplified relative to the plateau by an early phase whose EoS is stiffer than radiation (i.e., w > 1/3). This possibility has subsequently been studied by many authors [7,[23][24][25][26][27][28][29][30].
In a previous paper [7, hereafter LSR17], we investigated this amplification effect in the particular context of complex scalar field dark matter (SFDM) made up of charged ultralight bosons [31]. If all cosmological dark matter consists of SFDM (the ΛSFDM model), the Universe would be dominated at early times by the stiff phase of SFDM (w SF = 1), before the standard RD era. The stiff phase of a scalar field is also known as the "kination" phase [32], since the energy density of the SFDM is dominated by the kinetic energy. LSR17 showed that this early stiff-SFDM-dominated era (w = 1; "stiff era" for short) indeed amplifies the high-frequency part of Ω GW (f ) relative to the plateau value. This amplified SGWB may contribute a non-negligible radiation component to the total energy density, which boosts the expansion rate during the RD era. Meanwhile, this same effect results in a blue tilt in Ω GW (f ), which may even make direct detection of the SGWB possible at high frequencies by current laser interferometer experiments, e.g., Advanced LIGO-Virgo [33,34] and LISA [35]. Therefore, the stiff-era amplification effect (henceforth, "stiff amplification") encourages multi-wavelength search of the primordial SGWB using different GW probes [e.g., 36,37]. In this paper, we again focus on the stiff-amplified primordial SGWB, in a more general context, not limited to that involving SFDM.
Besides the CMB and laser interferometers, pulsar-timing array (PTA) observations can probe the SGWB by searching for correlated timing deviations in millisecond pulsars induced by the SGWB [38,39]. Recently, the North American Nanohertz Observatory for Gravitational Waves (NANOGrav) reported strong evidence for a stochastic common-spectrum process in their 12.5 yr pulsar-timing data set [40] with a high amplitude (h c ∼ 10 −15 at f yr = 1 yr −1 ). Though it has not been confirmed as an SGWB detection yet, many interpretations in this direction have flourished since then. Possible SGWB sources include a cosmic population of supermassive black hole binaries [41], cosmic strings [42][43][44], phase transitions [45][46][47][48], the primordial SGWB with a large initial blue tilt from non-standard inflationary scenarios (relaxing the consistency relation) [49,50] and others [e.g., 51,52].
In this paper, we are, however, interested in the secondary blue tilt in the primordial SGWB produced by stiff amplification, within the standard inflationary scenario which preserves the consistency relation (henceforth, the "standard inflation + stiff amplification" scenario). The case of stiff amplification (w = 1) is the one that maximizes the possible secondary blue tilt that results for modes that reentered when the EoS of the Universe has w > 1/3. Thus, the first part of this paper is dedicated to the question of whether stiffamplified primordial SGWB can explain the high common-spectrum amplitude reported by NANOGrav. To this end, we consider a cosmological model with two additional components to the base-ΛCDM model [53]: a stiff component and the primordial SGWB. In our model, when inflation ends, there is an extended phase of reheating with a matter-like EoS (w = 0) [54,55]. When reheating ends, the Universe is assumed to be dominated by the stiff component and remains so until the onset of the RD era. In order to constrain our model parameters, we perform a joint analysis of the latest observational results from the CMB, BBN, NANOGrav and Advanced LIGO-Virgo's third observing run (O3) [56].
Our analysis has a novel feature: we self-consistently include the backreaction of the SGWB on the background expansion rate, as we did in LSR17. Although noted before [e.g., 22,57], this backreaction effect is unfortunately often neglected when modelling the SGWB. Nevertheless, as stated above, the stiff-amplified SGWB can contribute a non-negligible (percent level) radiation component to the total energy density during the RD era. This will, in return, not only affect the evolution of tensor modes, but also other observables, e.g., radiation-matter equality and the CMB damping tail. A precise analysis of the primordial SGWB ought to account for its coupling with the background expansion history, therefore.
In the meantime, the well-known Hubble tension [e.g., 58, 59] also motivates our treatment. The present-day Hubble constant, H 0 , now measured at better than 3% precision by several experiments, shows a discrepancy (> 3σ) between its value measured by the CMB [53] and that by the distance ladder or time delays of lensed quasars in the nearby Universe [60,61]. With respect to the aforementioned radiation-matter equality, one way to alleviate the Hubble tension is to exploit the H 0 − N eff degeneracy: the redshift of this equality can be kept constant by increasing the value of H 0 and the effective number of relativistic species at the same time [62][63][64]. Our model implements this H 0 − N eff degeneracy, boosting H 0 in accordance with the additional radiation-like SGWB contribution, while the coupled evolution of the Hubble parameter and tensor modes is properly taken into account. Thus, the second part of this paper is dedicated to the implication of current constraints on the primordial SGWB for the Hubble tension. We investigate the extent to which the stiff-amplified SGWB can bring the value of H 0 from the CMB into agreement with those from local measurements.
The paper is organized as follows. In section 2, we demonstrate the stiff amplification effect on the primordial SGWB and introduce our model. In section 3, we discuss all current measurements and upper bounds on the primordial SGWB, for each of several probes in turn. In section 4, we combine these probes in a joint analysis and derive the constraints on the "standard inflation + stiff amplification" scenario that result. The implication of these results for the Hubble tension is explored in section 5. We conclude in section 6.

Stiff amplification of the primordial SGWB
In this paper, we consider the primordial tensor perturbations with respect to a flat FLRW background metric, so the short-wave, weak-field limit is apparently satisfied for the GWs described by these tensor modes (see appendix A). We can write down the perturbed metric in the transverse-traceless gauge (the "TT gauge") [10,57], In section 2.1, we review the basic equations concerning the primordial SGWB from inflation, and its amplification by a post-inflationary stiff era. In section 2.2, we present our cosmological model for the "standard inflation + stiff amplification" scenario, which selfconsistently includes the stiff-amplified primordial SGWB.

Basic equations
Primordial tensor perturbations can be expanded in Fourier space [e.g., 65], where f is the comoving frequency,k is a unit vector, k ≡ 2πfk/c, and P ij are the polarization tensors for the + and × states. In our convention, f,k)) * due to the reality of h ij . When a mode is well-inside the Hubble radius, h P (t, f,k) ∝ e −2πif η /a, where η is the conformal time, dη ≡ dt/a. 2 For an isotropic, stationary and Gaussian SGWB, the most straightforward observable is the two-point correlation function. In Fourier space, it is defined as where δ S 2 D is the Dirac function on the two-sphere and S h (f ) is the one-sided power spectral density of the SGWB [67]. S h (f ) is related to the characteristic amplitude/strain of the SGWB, h c (f ), and the (dimensionless) tensor power spectrum, The primordial SGWB is characterized by its power spectrum at an initial time, . Standard singlefield, slow-roll inflation predicts a nearly scale-invariant initial power spectrum for tensor modes, Here A t is the tensor amplitude, r ≡ A t /A s defines the tensorto-scalar ratio, and the tensor spectral index satisfies the consistency relation, n t = −r/8. Following the convention of Planck, the pivot scale is chosen as k * = 2πf * /c = 0.05 Mpc −1 [68]. As for the transfer function, its evolution follows from the wave equation (A.2), where the overdot denotes the derivative with respect to the cosmic time, t, and we have omitted the arguments (t, f ) for brevity. For any mode that has undergone inflation, its amplitude is frozen while it is well-outside the Hubble radius, so that ∆ 2 h (t, f ) ∆ 2 h,i (f ), T h 1 andṪ h 0. After Hubble reentry, the transfer function for all modes asymptotically evolves as T h ∝ 1/a (the adiabatic solution). However, their relative amplitudes (frequency dependence) at a given time is subject to the EoS of the Universe at the reentry of each mode, allowing for stiff amplification. We will recap this effect in what follows.
For a given mode of frequency f , its sub-Hubble solution for Ω GW (f ) can be approximated as (using eq. [A. 6]) where H ≡ȧ/a is the Hubble parameter and a f H(a f ) ≡ 2πf defines the scale factor at its Hubble reentry (cf. eq. [58] in LSR17). As mentioned in the introduction, modes that reentered during the RD era correspond to a plateau in Ω GW (f ) for standard inflation. Alternatively, an early stiff era gives rise to a blue-tilted spectral shape in Ω GW (f ). Such a stiff era is proposed by a variety of physical mechanisms [e.g., 69,70] and many of them involve a scalar field dominated by its kinetic energy [23,32,[71][72][73][74]. In LSR17, the stiff era is due to the stiff phase of SFDM, interposed between reheating and the RD era. We illustrate this stiff phase by an example ΛSFDM universe in appendix B. For a mode that reentered during the stiff era, we showed in section III.B.3 of LSR17 that its Hubble reentry happens later than it would if the Universe were RD all the time. In other words, a f,stiff > a f,rad for the value of a f appearing in eq. (2.4). Therefore, eq. (2.4) shows that for such a mode, the value of Ω GW (a, f ) in the sub-Hubble limit (a a f ) is greater than it would be if the Hubble reentry happened in the RD era (i.e., the plateau value). In this way, the primordial SGWB is amplified for modes reenter during the stiff era, relative to the plateau.

Stiff-amplified SGWB: self-consistent model for precision cosmology
Stiff amplification causes a secondary blue tilt in Ω GW (a, f ) evaluated at late times when all modes of interest are in the sub-Hubble limit. Whereas any pre-RD era with an EoS stiffer than radiation would generically lead to a blue tilt (whose spectral index depends on the EoS), we will only consider a stiff era (w = 1) in this paper, in order to maximize the amplification. Then, for modes that reentered during the stiff era, On the other hand, an extended period of reheating, with a matter-like EoS (w = 0), precedes the stiff era, as mentioned above. For modes that reentered during reheating, Ω GW (f ) ∝ f −2 . Therefore, in the "standard inflation + stiff amplification" scenario, the combined effect of reheating and the stiff era introduces an excess in the spectrum of Ω GW (f ) relative to the plateau associated with the standard RD era, which appears as a triangle (in logarithmic scales; see LSR17). This triangle peaks at f re , which corresponds to the mode that reentered at the end of reheating, characterized by T re , the reheating temperature.
To account for stiff amplification, we consider a cosmological model which contains a stiff component (w s = 1) and the primordial SGWB, in addition to all the base-ΛCDM components. 3  , derived from the numerical solutions to the dynamical system described in the text for those parameters, are listed here as well. These observables will be discussed in section 3.
the stiff component. Thereafter, the energy density of the stiff component evolves as ρ s ∝ a −6 and dominates the total energy density of the early Universe between the end of reheating and the end of the stiff era. The latter endpoint is defined as the moment of equality between the energy density of the stiff component and that of the radiation components, parameterized by the temperature at this equality, T sr . Therefore, apart from the base-ΛCDM parameters, our model has three parameters: r, T re and T sr . As we shall describe below, the model requires us to solve a set of coupled, integro-differential equations for each set of model parameters.
To solve for tensor transfer functions, we apply the dynamical system approach. For a given mode with comoving frequency f , the following dynamical variables can be defined: Apparently, T h = y f /e ζ f . The wave equation (2.3) can then be rearranged into the following dynamical system: 6a) where the prime denotes the derivative with respect to the number of e-foldings, N ≡ ln a (dN = H dt), and where Ω GW and Π GW are defined in eqs. (A.4) and (A.5), and Ω m , Ω r and Ω s are the energy fractions of matter (CDM+baryons), radiation (photons+massless neutrinos) and the stiff component, respectively. Apparently, σ is related to the EoS of the Universe by σ = 1 + w. Therefore, the evolution of each tensor mode is coupled to the expansion history of the background Universe via σ.
For illustrative purposes, we have solved these equations above for several example models, with parameters listed in table 1. While we will say more about our numerical method Matter Λ a eq a sr BBN below, it is useful to present these examples first, in order to anticipate the general behavior of the solutions in the discussion which follows. The left panel of figure 1 shows the energy density evolution of each component for Model I. The time evolution of σ is illustrated in the right panel of figure 1, for Model I in table 1. In order for amplification to take place, σ = 4/3 is required. As mentioned in the introduction, when stiff amplification (σ = 2) of the primordial SGWB occurs, the coupling between the radiation-like SGWB and the background metric may cause significant backreaction from the SGWB on the Hubble parameter.
To account self-consistently for this backreaction, we must solve the coupled dynamical system of eqs. (2.6) and (2.7) for each frequency, for any given set of model parameters, (r, T re , T sr ). Our method of solution is described in appendix C. Ordinarily, the solution of these coupled equations would be subject to boundary conditions at the present, fixed by the observational values adopted for Ω m,0 and H 0 (where Ω Λ,0 = 1 − Ω m,0 for a flat FLRW Universe). However, observations of the CMB and baryon acoustic oscillations (BAO) also fix the value of the redshift of radiation-matter equality, z eq , to an exquisite precision [e.g., 53].
Since the SGWB adds an extra radiation component to the background energy density, we must ensure that our solution yields the observed z eq , despite this. In so doing, we encounter the degeneracy between the value of H 0 measured by the CMB and BAO and the boost to the radiation energy density by the SGWB, allowed by the requirement that z eq is fixed. 4 As we shall show, by the end of the stiff era, the contribution of the SGWB to the background energy density reaches an asymptotic value, relative to that of the other radiation components. This asymptotic ρ GW (which thereafter evolves as ρ GW ∝ a −4 ) can be represented by a constant value of ∆N eff , the effective number of extra relativistic species. ∆N eff ≡ N eff − N eff,0 , where N eff,0 = 3.046 corresponds to three Standard Model neutrinos [75]. As a result, we are able to utilize the H 0 − N eff degeneracy for which the value of z eq is preserved, to determine the boundary conditions in our solutions, as follows. While H 0 and N eff can both vary, we keep Ω m,0 and Ω r,0 + Ω GW,0 fixed, thus fixing z eq and z mΛ (the redshift of matter-Λ equality). In our model, then, the H 0 − N eff degeneracy is stated as where H 0,ΛCDM is the value in the base-ΛCDM model (for which ∆N eff = 0). Thus, our model can actually help alleviate the Hubble tension by boosting H 0 . Figure 2 illustrates the necessity of our treatment for the backreaction of the SGWB. The left panel shows that for Model III in table 1, the extra radiation due to the stiffamplified SGWB can indeed cause a ≈ 3% increase in the Hubble parameter during the RD era (∆N eff ≈ 0.37). The right panel shows that the backreaction of this SGWB would lead to a shift of z eq more than 6σ away from its value measured by Planck using the base-ΛCDM model [53], if Ω m,0 and H 0 were both fixed at the ΛCDM best-fit values. In short, precision cosmology requires that the simultaneous backreaction of the primordial SGWB on the background expansion history be self-consistently taken into account throughout its evolution. We have confirmed that our treatment meets this requirement with a precision ∼ 10 −3 (cf. appendix C), for all viable model parameters, (r, T re , T sr ).
Our results for the present-day SGWB energy spectra, Ω GW (f ), for the example models in table 1, are shown in figure 3. These models are chosen so as to illustrate the dependence of the spectral shape on the model parameters. To begin with, they all share a plateau in Ω GW (f ) of the same height since they assume the same value of r. Models I -V all display the blue tilt and the triangle-shaped spectrum at high frequencies due to stiff amplification. 5 Model I has the lowest value of T sr and thus the highest amplitude at f yr ≡ 1 yr −1 , the . Present-day SGWB energy spectra for all example models in table 1. Vertical dashed lines indicate the representative frequencies of each SGWB probe, f * , f yr and f LIGO for the CMB, PTA and LIGO-Virgo, respectively. The 5%-95% confidence interval of the common-spectrum amplitude reported by the 12 yr NANOGrav results (labeled "NG12") is displayed [40], along with the 95% upper limit from BKP [68] and that from Advanced LIGO-Virgo O3 [56].
reference frequency for PTAs. This example also shows that when the end of the stiff era slightly overlaps BBN, the value of N eff at BBN can be different from that at late times. Models I and II have different values of T re and T sr but the same "area" under the triangle in Ω GW (f ) (as if the same triangle slides along the plateau), which manifests itself in the equal values of N eff for these two models at late times. Models II and III have the same value of T sr , so their stiff eras end at the same time.  The cross-sectional planar slices of the 3-D space of model parameters chosen in each view are colorcoded, and the grey region represents parameters entirely excluded by the observational constraints. The vertical black line indicates the 95% CL upper limit on r from BKP [68]. The magenta bars on the color box (labeled by "NG12") indicate the 5%-95% confidence interval of the common-spectrum amplitude, h c (f yr ), from the 12.5 yr NANOGrav results [40].

CMB temperature and polarization
The primordial SGWB can leave an observable imprint on the CMB temperature and polarization anisotropy [76][77][78]. In particular, detection of the CMB B-mode polarization around ∼ 100 would be a convincing signature of the primordial SGWB. Currently, BICEP2/Keck Array+Planck (BKP) only provides an upper bound on the tensor-to-scalar ratio, r < 0.061 at 95% confidence level (CL) [68]. This upper bound directly applies to our model, too, since the stiff era does not affect long-wavelength modes that reentered around recombination. In the future, CMB-S4 experiments will continuously seek to measure the primordial SGWB from inflation [79].

NANOGrav results
PTA observations measure the times of arrival ("ToAs") of radio pulses from millisecond pulsars. Those ToAs can be modulated by an SGWB permeating the spacetime between the pulsar and the earth. In fact, the existence of an SGWB would be manifested in the timing-residual cross-power spectral density (cf. eq. [2] in [40]) as a time-correlated, commonspectrum stochastic process across all pulsar-earth pairs, with quadrupolar spatial correlations between pulsars (i.e., the Hellings & Downs curve [80]). In PTA analysis, the characteristic strain of an SGWB is usually modeled as a power law, NANOGrav recently discovered a time-correlated, stochastic process with a common amplitude and spectral index in their 12.5 yr data set. However, there is little evidence for quadrupolar spatial correlations in this common-spectrum process, required to identify it with an SGWB. Hence, the NANOGrav results are still inconclusive with regard to GW detection, and, in the meantime, have yet to be confirmed by other PTAs. Nevertheless, despite its uncertainty, this reported common-spectrum process has incited many attempts to explain it in terms of the SGWB.
In our model, the present amplitude of h c , or, equivalently, Ω GW (f ), at frequencies near f yr can be higher than in the ΛCDM model only if the corresponding modes have experienced Present-day energy density fraction per logarithmic frequency of the stiff-amplified primordial SGWB, Ω GW (f ), at the reference frequency f LIGO = 25 Hz in our model, presented as the three-view projections with respect to the model parameters, (r, T re , T sr ). The cross-sectional planar slice of the 3-D space of parameters shown in each view is color-coded, and the grey region is entirely excluded by the observational constraints. The vertical black line indicates the 95% CL upper limit on r from BKP [68]. The magenta curves indicate the 95% CL upper limit on Ω ref, LIGO from the Advanced LIGO-Virgo O3 results [56]. stiff amplification. For example, as shown in figure 3, these modes lie within the blue-tilted part of the SGWB spectrum for Model I, so their amplitudes at f yr are higher than the ΛCDM-plateau value. Here, we sample our model parameters, (r, T re , T sr ), throughout the entire parameter space, to calculate the value of h c (f yr ) of the primordial SGWB for all model parameters of interest. Our results are shown in figure 4 (with the T sr axis upside-down in all figures in this paper that show the results for different models, spanning the range of parameter space). 6 We compare our results with the A CP posterior reported by NANOGrav [40], for which the 5%-95% confidence interval is 1.75 − 3.83 × 10 −15 in the case of the bluetilted spectral slope predicted for the stiff-amplified SGWB (α = −1/2). 7 Figure 4 shows that the amplitude of the SGWB at f yr in the "standard inflation + stiff amplification" scenario, as constrained by other observations, is too small to explain the common-spectrum process in the 12.5 yr NANOGrav data set. 8

Advanced LIGO-Virgo
Laser interferometers like the Advanced LIGO-Virgo network can directly detect SGWBs by cross-correlating data from different detectors [e.g., 82]. Recently, the LIGO Scientific Collaboration and Virgo Collaboration published results of a search for an isotropic SGWB 6 Throughout our analysis, we do not sample the grey region in parameter space displayed in figures 4 -7 for computational efficiency, because models in this region result in too much extra radiation energy density from the stiff-amplified SGWB, and are thus firmly excluded by late-Universe N eff bounds (cf. section 3.4).
7 These values we quote are different from those in the fiducial model in NANOGrav's analysis, because the latter assumes α = −2/3, as expected for the SGWB from unresolved mergers of supermassive black-hole binaries. 8 Our result here that the amplitude of the stiff-amplified SGWB spectrum at fyr is constrained to be far below the NANOGrav results is qualitatively consistent with the argument in [81], based upon applying the BBN constraint (which we shall discuss in section 3.4) to limit how late the stiff era can end. While the model in [81] differs from ours (e.g., it posits a stiff era that immediately follows inflation, with no standard reheating process), this reflects the fact that the example GW spectra in [81], computed numerically, share the spectral feature of ours for modes whose Hubble reentry occur during the stiff era, with a blue tilt of ΩGW(f ) ∝ f . using data from their first three observing runs (O1, O2 and O3) [56]. While the crosscorrelation spectrum from data does not show evidence for an SGWB signal, a new upper limit is placed on the present-day SGWB energy spectrum, modeled as a power law, The reference frequency is chosen to be f LIGO = 25 Hz.
We again calculate the value of Ω ref, LIGO in our model, sampling the model parameters (r, T re , T sr ), as shown in figure 5. Since the stiff-amplified SGWB in our model has a triangleshaped spectrum (i.e., Ω GW (f ) is a broken power law), it does not always have a fixed spectral index across the LIGO-Virgo frequency range, 20-1726 Hz. Therefore, we compare our results with the marginalized 95% CL upper limit from the O3 analysis, Ω ref, LIGO < 6.6 × 10 −9 , obtained by integration over α LIGO . Figure 5 displays this upper limit.

Integral bounds: BBN, CMB+BAO
The primordial SGWB can also be searched by indirect probes, e.g., light element abundances from BBN, the CMB, and large-scale structure of the Universe. These cosmological probes provide what is known as integral bounds on the SGWB, since the observables in each probe are (indirectly) affected by the integration of Ω GW (f ) over a wide range of frequencies. In the following, we will examine all such current probes, classifying them according to the epoch in the expansion history of the Universe to which each probe is sensitive.
Early-Universe cosmology: big bang nucleosynthesis. Standard BBN predicts certain relic abundances for light elements like helium-4 and deuterium (see [83] for a brief review). These abundances are sensitive to the cosmology of the background Universe during BBN (when T ∼ 10 9 K), in particular the baryon-to-photon ratio and the expansion rate then. Therefore, one can infer related cosmological parameters, namely the baryon density, Ω b,0 h 2 (where we use h here to mean the Hubble constant in units of 100 km s −1 Mpc −1 ), and the effective number of relativistic species at that time, N eff, BBN , by combining observations of the primordial 4 He and D abundances with theoretical BBN calculations [e.g., 84]. We, henceforth, use N eff, BBN to denote its value during BBN, in order to distinguish it from the value in the late Universe, N eff, late (which affects the CMB and BAO). We note that, in our discussion in section 2.2 of the asymptotic ∆N eff associated with ρ GW , we were actually referring to this latter N eff, late . By contrast, the value of N eff, BBN in our model may have contributions from both the primordial SGWB and the stiff component, the latter because it increases the expansion rate of the Universe relative to the rate for a standard RD Universe with three neutrino species, even though it does not, itself, evolve like a radiation-like component. As a result, the constraint on N eff, BBN can be translated into constraints on the sum of the stiff-amplified primordial SGWB and the stiff component (rather than on the SGWB alone) in our model, and thus on the model parameters.
In this paper, we quote the 95% CL bounds on N eff, BBN , marginalized over Ω b,0 h 2 , obtained from combining measurements of the primordial 4 He mass fraction, Y P , and the primordial deuterium abundance, (D/H) P . For the Y P measurement, our baseline is from the data compilation of [85] (A15), but we also quote the bounds from [86] (I14). For the (D/H) P measurement, we reference the results from [87] (C14). 9 The combined observational bounds on N eff, BBN are presented as follows: The discrepancy between them is due to the moderate tension between the Y P measurements from A15 and I14, which is still under debate. It is worth noting that the lower bound from I14+C14 slightly disfavors the standard value N eff,0 = 3.046.
We have calculated the value of ∆N eff, BBN in our model for each choice of model parameters. The results are shown in figure 6, where the 95% CL upper limits from each of the two combined observational constraints presented in eq. (3.1) are also displayed.
Late-Universe cosmology: radiation-matter equality and CMB damping tail. Extra radiation components that survive through the late Universe, like the stiff-amplified primordial SGWB, may significantly contribute to the value of N eff, late . In our model, unlike the case of N eff, BBN discussed above, the SGWB is the only contribution to ∆N eff, late , since the stiff component is negligible in the late Universe. Therefore, constraints on N eff, late from the CMB and BAO can be translated into constraints on the stiff-amplified primordial SGWB, and thus on our model parameters. Before we discuss these late-Universe constraints, it is important to (1) point out that a general analysis should distinguish them from the BBN constraints above (which are only concerned with the early-Universe cosmology), and (2) clarify which kind of analysis is suitable for this purpose.
Currently, the CMB+BAO data sets are compatible with primordial element abundance data inasmuch as they can be jointly fitted by an extension to the base-ΛCDM model with a fixed, time-independent N eff which is common to both the early epoch of BBN and the late times probed by the CMB and BAO observations (i.e. by assuming ∆N eff, BBN = ∆N eff, late , henceforth "ΛCDM+N eff ") [53]. While many analyses are based on this assumption [e.g., 64], the concordance between BBN and the late-Universe cosmology, however, does not demand Figure 7. Effective number of extra relativistic species in the late Universe, ∆N eff, late , in our model, presented as the three-view projections with respect to the model parameters, (r, T re , T sr ). The crosssectional planar slices of the 3-D space of model parameters chosen in each view are color-coded, and the grey region is entirely excluded by the observational constraints. The vertical black line indicates the 95% CL upper limit on r from BKP [68]. The magenta curves indicate the 95% CL upper limit on ∆N eff, late from combining the CMB and BAO data with a prior for Y P (eq. [3.2a], our baseline). The cyan curves indicate the 95% CL upper limit from CMB+BAO only (eq. [3.2b]). that these two values of N eff be equal. As a matter of fact, viable cosmological models like ours allow for different values of N eff at BBN and at late times (cf. Model I in table 1). In light of such models, it is therefore more general and favorable to constrain N eff, BBN and N eff, late separately. For the latter, the analysis should only involve physical processes that directly determine the late-Universe observables, independent of BBN, in a clean way. We carefully examine those physical processes in the following: (1) The physical size of the sound horizon depends on the duration of the RD era. Thus, the angular scales of the sound horizon measured by both the CMB and BAO (at different redshifts) are sensitive to the value of z eq .
(2) The early Integrated Sachs-Wolfe (ISW) effect refers to the enhancement of the CMB temperature anisotropies due to the time-variation of gravitational potentials after recombination, when the Universe was still not yet fully transitioned from RD to MD. In particular, the relative heights of the first three peaks in the CMB temperature power spectrum are sensitive to the value of z eq via the early ISW effect [e.g., 89].
(3) On even smaller scales, the CMB temperature anisotropies are damped by photon diffusion, an effect known as Silk damping. The slope of the damping tail in the CMB power spectrum reflects the amount of Silk damping [63]. Since it depends on both the expansion rate at recombination and the number density of free electrons, the CMB damping tail measurements are thus sensitive to both N eff, late and Y P .
The z eq measurements based on the first two physical processes above are still subject to the H 0 −N eff degeneracy, as described in section 2.2. This degeneracy can, however, be broken by the CMB damping tail measurements, because they provide additional information that enables constraining N eff, late on its own. Therefore, one can constrain N eff, late independently of information involving any early-Universe process (e.g., a BBN determination), by fitting the CMB+BAO data with an extended ΛCDM model which allows both N eff, late and Y P to vary freely (henceforth, "ΛCDM+N eff +Y P "). In this paper, we quote the N eff, late constraints from such an analysis provide by the Planck 2018 results (CMB+BAO). 10 Optionally, the Y P value from helium abundance measurements can be additionally combined as a prior in the analysis. This provides a tighter constraint on N eff, late while still independent of BBN [53]. We choose this constraint as our baseline in the paper.
The N eff, late constraints from Planck are quoted as follows:

Results: joint constraints on standard inflation + stiff amplification
In this section, we combine the constraints on our model parameters, (r, T re , T sr ), from all the probes of the primordial SGWB described in the previous section, to obtain the joint constraints on the "standard inflation + stiff amplification" scenario. We summarize these joint constraints in figure 8, which we shall now discuss.
As already described in section 3.2 and shown by figure 4, the amplitude of h c (f yr ) for the stiff-amplified SGWB in our model, while enhanced with respect to that without stiff amplification, in the base-ΛCDM model, is still much lower than the common-spectrum amplitude from the NANOGrav 12 yr data set. This is due to the fact that the model parameters required to amplify the primordial SGWB at this frequency to the level of the NANOGrav signal would result in excessively large values of ∆N eff, late , well above its current 95% CL upper bound from observations. In fact, the combined constraints from all other probes indicate that the difference in h c (f yr ) between our model predictions and the NANOGrav results is more than two orders of magnitude. Therefore, the "standard inflation + stiff amplification" scenario cannot explain the common-spectrum process reported by NANOGrav. If the latter were indeed due to an SGWB, astrophysical sources are more likely to be the origin.
In that case, as long as we remain below the NANOGrav results, our model can still be constrained by other probes of the primordial SGWB, e.g., laser interferometric experiments and indirect probes. The joint constraints on our model parameters from these bounds are shown in figure 8. It displays three-dimensional views of the 95% CL constraints in this parameter space, first as required to satisfy each observational constraint separately, for the constraints from the O3 data of the Advanced LIGO-Virgo, the N eff, BBN measurements, and the N eff, late measurements, respectively, 11 and then with a view of the range of parameters allowed by all three of those constraints -the 95% CL joint constraint. Unfortunately, since not all the likelihood data from these measurements are publicly available, we cannot yet perform a full Bayesian joint analysis to obtain the posteriors for our model parameters. Instead, our joint analysis here simply combines the 95% CL constraints on our model parameters from Figure 8. Three-dimensional view of the constraints on the "standard inflation + stiff amplification" scenario in its parameter space, (r, T re , T sr ). Panels (a)-(c): constraints from the LIGO-Virgo O3 results, the N eff, BBN measurements and the N eff, late measurements, respectively. For each probe, the constraints are visualized as the isosurface and contours for the corresponding 95% CL upper limit. In panels (b) and (c), the isosurfaces are from our baseline constraints on N eff, BBN and N eff, late , eqs. (3.1a) and (3.2a), respectively. In panels (a)-(c), the color-coded planar cross-sections are the same as those in the leftmost panels of figures 5-7, respectively. In each panel, the three thick magenta curves are the same as those magenta curves shown in the three views of the corresponding figure (among figures 5-7), and the grey vertical plane (with black borders in the figure) indicates the 95% CL upper limit on r from BKP [68]. Panel  Table 2.
Overall 95% CL allowed range of our model parameters (r, T re , T sr ), described as the 95% CL lower limit on T sr for given values of r and T re (i.e., the top surface of the allowed region in panel [d] of figure 8). The overall allowed range has r < 0.061 (95% CL) from BKP [68]. each probe, i.e., combining the isosurfaces in panels (a)-(c) of figure 8. The resulting 95% CL allowed range of (r, T re , T sr ) is indicated by the light green volume in panel (d).
To describe the features of this overall allowed range, we first note that there must be a lower bound on T re to allow BBN to occur, T re 4 MeV, and there is an upper bound on r from the CMB, r < 0.061 (95% CL). For fixed r and T re , a lower value of T sr (i.e., larger Ω s,0 ) implies longer duration of the stiff era and thus higher degree of stiff amplification of the primordial SGWB. Therefore, there must be a lower bound on T sr for given values of r and T re . As a matter of fact, this bound is described by the top surface of the allowed region in panel (d) of figure 8, the 95% CL lower limit on T sr (as a reminder, the T sr axis is upside-down is this figure).
Using this description, we find that the parameter range allowed by the joint constraints can be characterized by dividing it into three regimes, according to which observational probe of our model dominates in each regime. They can be roughly parameterized by the range of T re . The results are laid out in table 2 and labeled in panel (d) of figure 8. We can describe these three regimes as follows, itemized by the regime number: (i) This regime has the lowest values of T re . The dominant constraint on T sr is from N eff, BBN (cf. the discussion on the N eff, BBN constraint in section 3.4, where we explain that in our model, not only the SGWB but also the stiff component can contribute to N eff, BBN significantly). In fact, the lower limit on T sr in this regime is roughly manifested as a horizontal plane in the parameter space, insensitive to the values of r and T re . It reflects the fact that in this regime, the SGWB is unimportant and our model is mainly constrained by the requirement that the stiff-to-radiation transition must finish early enough, so that the stiff component alone, for which ρ s ∝ a −6 , does not boost the expansion rate during BBN beyond the observational constraints.
(ii) This regime has the intermediate range of T re . The dominant probe is the LIGO-Virgo measurements, since the frequency of the peak in Ω GW (f ) in our model, f re -which corresponds to the mode that reentered the Hubble radius at the end of reheating (the beginning of the stiff era) -is around f LIGO = 25 Hz then. The lower limit on T sr for given values of r and T re in this regime is manifested as the "waterfall" surface shown in panel (d) of figure 8.
(iii) This regime has the highest values of T re and the dominant constraint is from N eff, late , which amounts to an upper bound on the area under the triangle in Ω GW (f ) (cf. figure 3). Correspondingly, the lower limit on T sr for given r and T re in this regime roughly appears as a plane in the parameter space, whose equation is specified in table 2.
By contrast with our model predictions for the PTA frequency range, which fall far short of the reported NANOGrav signal for our allowed model parameters, there is no such gap between our model and the current upper limits from all other observational probes. As a result, if the NANOGrav signal holds up over time and is confirmed, then these other probes (i.e., LIGO-Virgo, N eff, BBN and N eff, late ) may be more likely to detect the "standard inflation + stiff amplification" scenario than PTA. For example, for the results presented here, within the range of model parameters allowed by the joint analysis, any future detection by LIGO-Virgo of the primordial SGWB, consistent with its current O3 95% CL upper limit, can be explained by our model. If so, then our model would provide an explanation within standard inflation which does not require an initial spectral tilt.
Each of these probes will be continually upgraded in the future for its sensitivity (e.g., CMB-S4 [79], LIGO A+ 12 ), and a comparison among those sensitivities would be required to determine which probe will provide the first evidence of this early-Universe scenario. 13 If detection results for multiple probes, all consistent with the predictions of "standard inflation + stiff amplification", this would then constitute smoking-gun evidence in favor of the model.

Implications for the Hubble tension
New observational results continue to increase the significance of the tension between H 0 measurements from CMB+BAO and those from the nearby Universe. This has motivated growing interest in alternatives to the base-ΛCDM model [e.g., 61]. Here we examine the possibility of reconciling/alleviating the Hubble tension in our "standard inflation + stiff amplification" scenario, by the presence of the stiff-amplified primordial SGWB as an extra radiation component. In the presence of this extra radiation component, current measurements of z eq (for a precision ∼ 0.6%) then drive us to incorporate the H 0 − N eff degeneracy [64] in our analysis, expressed as eq. (2.8) in our model. We note that the H 0 − N eff degeneracy is concerned with N eff, late , the late-Universe value of N eff , independent of N eff, BBN .
The resulting impact of our model on the Hubble tension is illustrated in figure 9. , BAO+BBN [64] and CMB+BAO [53]. The last two measurements are both based on the "ΛCDM+N eff " model, and both involve a BBN calculation which assumes N eff, BBN = N eff, late . In contrast, our analysis is free from these assumptions, as described in section 3.4. Since the H 0 − N eff degeneracy is only concerned with N eff, late , we adopt the BBN-independent upper bound on N eff, late from the CMB+BAO analysis provided by Planck (eq. [3.2]), based on the "ΛCDM+N eff +Y P " model [53]. Therefore, combining eq. (2.8) with the Planck fits, we obtain the corresponding upper limits on H 0 for our model, as follows: with Y P prior without Y P prior Figure 9. Hubble tension implications of the "standard inflation + stiff amplification" scenario. The tension is illustrated by the discrepancy between the H 0 measurements from the nearby Universe, including SH0ES [60] and H0LiCOW [61], and those from the CMB and BAO, including BAO+BBN [64] and CMB+BAO [53]. The latter two analyses are both based on the ΛCDM+N eff model, and the solid ellipses are their 95% CL contours, respectively. All the shaded regions represent the 68% CL ranges for the respective measurements. With respect to our model, the H 0 − N eff degeneracy relation, eq.   figure 9. They occur for model parameters corresponding to the current 95% CL upper limits of N eff, late , and thus to the boundary surface for regime (iii) described above (cf. figure 8 and table 2).
Meanwhile, for these current observational limits on N eff, late , our model may be able to reduce the discrepancy between measurements of H 0 by CMB+BAO and SH0ES from 4.4σ to 3.8σ for the baseline upper limit (eq. [3.2a]), and that between CMB+BAO and the H0LiCOW measurement from 3.1σ to 2.8σ. In addition, if we take the more relaxed upper limit, eq. (3.2b), our model may further bring the H 0 from CMB+BAO to within 3.4σ of the SH0ES measurement, and within 2.6σ of the H0LiCOW measurement.

Summary and discussion
The recent NANOGrav 12 yr results have invited many interpretation attempts involving a primordial stochastic gravitational-wave background. All those attempts are based upon non-standard early-Universe scenarios which can predict stronger amplitudes in the PTA frequency range than that from standard inflation + ΛCDM. While it is understood by many authors that such an SGWB contributes an extra radiation component to the background Universe, which may therefore affect its expansion history, we here investigate the possibility that this extra radiation may then help alleviate the current Hubble tension, thus drawing a novel connection between gravitational waves and cosmology.
We demonstrate this by considering a cosmological model, the "standard inflation + stiff amplification" scenario, with two components added to the base-ΛCDM model: a stiff component and the primordial SGWB. In our model, an early stiff era (with w = 1) arises, when the stiff component dominates, between the end of an extended period of reheating and the standard RD era. Unlike some other suggestions to explain the NANOGrav signal by postulating nonstandard inflation with a blue-tilted primordial spectrum for the tensor modes responsible for the SGWB, our model does not require us to depart from standard inflation so as to tilt the primordial spectrum. In our case, the primordial spectrum is nearly untilted, but a secondary blue tilt results, instead, from the stiff amplification caused by the stiff era. Clarifying its distinction from parametric amplification, we revisit this stiff-amplification effect on the primordial SGWB under the general scenario considered here, which has three parameters, (r, T re , T sr ). The secondary blue tilt in the stiff-amplified primordial SGWB is manifested in its present-day energy spectrum as Ω GW (f ) ∝ f , for the frequency range of modes that reentered the Hubble radius during the stiff era. In this paper, we address the questions of whether such a blue tilt may explain the NANOGrav results and to what extent the stiff-amplified primordial SGWB can reduce the Hubble tension. Along the way, we make predictions for other direct and indirect observables of the SGWB, as well.
In doing so, we develop a new method to include the backreaction of the SGWB on the background expansion rate self-consistently. In fact, we point out that any GW analysis based on a model that can significantly boost the SGWB, like ours, must account for its backreaction on the background Universe, in order to preserve the well-measured redshift of radiation-matter equality as precision cosmology demands. For that, we solve the fullycoupled dynamical system of the SGWB and expansion history. We update its boundary conditions by boosting the Hubble constant in accordance with the extra radiation associated with the SGWB. In so doing, we utilize the H 0 − N eff degeneracy which preserves z eq .
We then sample the three-dimensional parameter space, (r, T re , T sr ), to perform a joint analysis of the NANOGrav results and the latest upper bounds from Planck, big bang nu-cleosynthesis and Advanced LIGO-Virgo, to constrain the model. We find that the resulting blue-tilted, stiff-amplified SGWB is still too small to explain the common-spectrum amplitude reported by NANOGrav (by at least two orders of magnitude), when constrained by current upper limits of the other observables: Ω ref, LIGO , N eff, BBN and N eff, late . The latter together provide joint constraints on our model parameters. We find that the parameter range allowed by the joint constraints can be characterized by dividing it into three regimes, according to which observational probe of our model dominates in each regime.
While we have shown that, for its allowed parameters, the maximum amplitude of the predicted primordial SGWB for the "standard inflation + stiff amplification" scenario in the PTA frequency range is far smaller than the reported NANOGrav signal, there is no such gap between our model predictions and the current upper limits on the SGWB from all other observational probes. In the future, therefore, even if the NANOGrav signal is confirmed and too large to be the primordial SGWB predicted here, these other probes (i.e., LIGO-Virgo, N eff, BBN and N eff, late ) may still be able to detect our predicted SGWB. For example, for model parameters which satisfy the constraints derived here from our joint analysis, any future detection of the primordial SGWB by LIGO-Virgo that is consistent with the current O3 95% CL upper limit can be explained by our model. If so, then this would be an explanation within standard inflation which does not require an initial spectral tilt.
As the sensitivity of these probes increases in the future (e.g., CMB-S4, LIGO A+), the chances of detecting the primordial SGWB will increase, as well. Which probe is likely to provide the first evidence of this early-Universe scenario will depend upon the relative improvements among their sensitivities over time. If detections occur for multiple probes, each consistent with the predictions of "standard inflation + stiff amplification", then this would be smoking-gun evidence in favor of the model.
With regard to the Hubble tension, we have shown that the "standard inflation + stiff amplification" scenario may reduce the discrepancy between the measurement of H 0 by CMB+BAO for the baseline upper limit (eq. [3.2a]) and that by SH0ES, from 4.4σ to 3.8σ, and that by H0LiCOW, from 3.1σ to 2.8σ. Moreover, according to our analysis, if we take the more relaxed upper limit, eq. (3.2b), instead, then our model can bring the value of H 0 derived from CMB+BAO to within 3.4σ of its value from the SH0ES measurement, and within 2.6σ of its value from the H0LiCOW measurement. Hence, our results demonstrate that, while existing attempts to reconcile the Hubble tension often appeal to an extra relic radiation component, the primordial SGWB as generally present in the current cosmological paradigm can provide a favorable candidate for such extra radiation which may at least partially reduce the tension. In fact, given the unknown expansion history of the Universe between the end of inflation and BBN, non-negligible extra radiation is a natural consequence of the primordial SGWB produced within the standard inflationary paradigm, caused by stiff amplification, without further new ingredients, e.g., dark photons or early dark energy.
Our results for the "standard inflation + stiff amplification" scenario can also be contrasted with those for models recently proposed to explain the NANOGrav results as the primordial SGWB from nonstandard inflation. In these nonstandard inflation models, the consistency relation between r and the spectral index n t is relaxed, to allow the primordial spectrum to have a large initial blue tilt [e.g., 49,50]. In our model, by contrast, the inflationary consistency relation is obeyed, and, therefore, the primordial spectrum is nearly flat. We have shown that the SGWB from standard inflation must be well below the reported NANOGrav amplitude, even after stiff amplification, while the nonstandard inflationary models can match the level of that amplitude only by postulating a large enough initial blue tilt.
As a result, the latter models often need to place limits on the impact of reheating, either by restricting T re to be small (i.e., 10 3 GeV) or else invoking some nonstandard process; otherwise, the blue tilt required to match the NANOGrav amplitude would be so large as to violate the current BBN constraints. Since our stiff-amplified SGWB does not rise to the level of the NANOGrav results anyway, our model is not restricted to low values of T re .
A Primordial SGWB: the short-wave, weak-field limit The classical description of GWs is based on a clean separation of perturbations from the background metric [66,92,93]. When such a separation occurs in spatial dimensions, this is the "short-wave" limit: all the wavelengths of interest are much less than the typical curvature radius of the background. If perturbations are additionally small (the "weak-field" limit), one can expand the Ricci tensor around its background value. In this expansion, the first-order term, R (1) µν , yields the equation of motion of the GWs (the wave equation) and the second-order term, R (2) µν , yields their effective stress-energy tensor. The latter is defined as whereḡ µν is the background metric, R (2) ≡ḡ µν R (2) µν and . . . 3D denotes the spatial average over a scale greater than all modes of interest.
Primordial tensor perturbations over a flat FLRW background naturally satisfy the short-wave, weak-field limit. Thus, in the TT gauge, the wave equation for these perturbations takes the following standard form, 14 Also, T GW µν turns out to be homogeneous since it is a spatially-averaged quantity by definition. For tensor fluctuations produced by inflation, the spatial average is equal to the ensemble average . . . according to the ergodic theorem. Thus, primordial GWs from inflation constitute a stochastic background.
The energy density and pressure of this SGWB can be explicitly written as Therefore, moving into Fourier space, the dimensionless energy and pressure spectra can be written as (A.4) 14 We consider no sources with anisotropic stress in this paper. The inverse relations are apparently For modes well-inside the Hubble radius (2πf /aH 1), the high-frequency limit is satisfied in addition to the short-wave limit. In this case, the adiabatic solution for plane waves reads T h ∝ cos (2π if η)/a. It oscillates much more rapidly than the Universe expands. Thus, only time-averaged values over several oscillations, . . . t , can be measurable in practice. We then haveṪ 2 h (2πf /a) 2 T 2 h (the time-averaging notation for sub-Hubble solutions is omitted throughout the paper for brevity). This yields showing that sub-Hubble modes indeed evolve like radiation, w(a, f ) = 1/3. The last equality above also implies the following relation in the sub-Hubble limit: B Illustrative example of an early stiff era: ΛSFDM universe Figure 10 shows the energy density evolution of all the components in an example ΛSFDM universe, for which the cosmological dark matter is composed of ultralight (m ∼ 10 −22 eV/c 2 ) bosonic particles in Bose-Einstein condensate [31,LSR17]. This dark matter model is described by a complex scalar field, thus known as scalar field dark matter (SFDM). Complex SFDM is a variant of the fuzzy dark matter (which is otherwise described by a real scalar field). It generically undergoes a kination or stiff phase as the earliest stage of its dynamical evolution. As a result, the primordial SGWB from inflation is subject to stiff amplification and may then causes significant backreaction on the background Universe. The example model shown in figure 10 has a repulsive quartic self-interaction, which causes the radiation-like phase of SFDM as shown in panel (a). The radiation-like SFDM contributes yet another extra radiation component to the critical density of the Universe, manifested as the corresponding plateau in panel (c). Later on, SFDM transitions into the matter-like phase and becomes dark matter, responsible for cosmological structure formation. C Numerical scheme: approximation by model with constant ∆N eff As described in section 2.2, we must solve the coupled dynamical system of eqs. (2.6) and (2.7) for each frequency, for each set of model parameters, (r, T re , T sr ). 15 It is necessary, in fact, to solve the system of equations for all frequencies at once, since it is a set of integro-differential equations, in which Ω GW and Π GW in eq. (2.7) are both integrals over all frequencies. An iterative solution is required, in that case, since the integrated quantities at a given time are not known until the solution is known for each frequency and at all times. Moreover, the finite-difference scheme must contend with the requirement of resolving the high-frequency oscillatory behavior of the solution in time, which requires many small steps. Even if we only integrate the dynamical system exactly during the Hubble reentry for each mode (for about 10 e-foldings) and stitch that solution with its analytical super-Hubble and sub-Hubble asymptotes, the total solution to the system consisting of all frequencies can be costly for a single set of model parameters alone. In addition, in order to constrain our model parameters by comparing the solutions for different parameters with observational constraints, we must sample a large grid of representative points in the three-dimensional parameter space, (r, T re , T sr ), and find the solution to the coupled equations for each point.
Faced with these computational challenges, we have developed an efficient numerical scheme to solve the coupled system. It takes advantage of the fact that in the exact solution for cases with significant stiff amplification, the contribution of the SGWB to the background energy density reaches an asymptotic value by the end of the stiff era, i.e., a constant fraction of that of other radiation components. As a result, ρ GW can be represented in terms of a constant value of ∆N eff . During the stiff era, itself, the expansion history is not sensitive to this value, except that it determines when the stiff era ends and the Universe becomes RD. As such, if we knew what that asymptotic value of ∆N eff was going to be, we could approximate the entire expansion history and therefore the evolution of σ (cf. eq. [2.7]) quite well, by adopting this asymptotic ∆N eff and assuming that it is constant from the beginning of the integration forward in time. Since the evolution of σ solely determines the transfer of tensor modes, as eq. (2.6) implies, this method would yield a good approximate solution to the coupled equations. Unfortunately, we do not know this asymptotic ∆N eff in advance of solving the coupled equations. However, we can approach this value iteratively, if we have an initial guess for the value of ∆N eff . In the end, this approximation enables us to produce computationally-efficient solutions, with negligible differences from the exact solutions.
In what follows, we describe this approximate model for the treatment of the backreaction of the SGWB (mentioned in section 2.2 with other detail) and justify our use of it. First, we explain why the value of ∆N eff due to ρ GW is asymptotically constant in cases with significant stiff amplification. The degree of stiff amplification depends on the duration of the stiff era. If the latter spans more than a few e-foldings, the frequency range of modes that reentered the Hubble radius during the stiff era will extend more than a few orders of magnitude, accordingly. In this case, combining eq. (A.5) and the spectrum of the stiff-amplified SGWB (Ω GW (f ) ∝ f ), one can shown that ρ GW must be dominated by high-frequency modes that reentered at the beginning of the stiff era (cf. figure 3). Since the tensor modes within a fixed frequency range evolve like radiation in the sub-Hubble limit (cf. eq. [A.6]), the overall stiff-amplified primordial SGWB can therefore be well approximated by a radiation component with a constant ∆N eff , shortly after the onset of the stiff era. As a result, we can replace eq. (2.7) in the exact dynamical system by the following approximation: where Ω er is the energy fraction of this extra radiation component. The simultaneous coupling between the primordial SGWB and the background Universe can thus be approximated by the system of eqs. (2.6) and (C.1), which is more efficient to solve than the exact system. We solve this approximate system iteratively with an update on its present-day boundary conditions for each iteration. Particularly, we calculate the asymptotic value of ∆N eff associated with ρ GW from the last iteration, and then update the value of H 0 as part of the boundary conditions for the next iteration, using the H 0 − N eff degeneracy relation, eq. (2.8). As a reminder, the sum Ω r,0 + Ω er,0 that appears as a parameter in eq. (C.1) is fixed in our treatment so as to fix z eq (cf. section 2.2). It needs no update between iterations, therefore. We adopt the following convergence criteria for the iterative scheme that the fractional difference between the asymptotic values of ∆N eff from consecutive iterations is less than 10 −3 . Fortunately, only a few iterations are required in order to converge or else to reach the conclusion that the adopted model parameters be excluded (i.e., the grey region in figures 4-7).
Justification of our numerical scheme is illustrated in figure 11, using Model III in table 1 as an example. The left panel demonstrates the consistency between the exact model (with the primordial SGWB) and the corresponding approximate model (with its converged value of the asymptotic ∆N eff ≈ 0.37), in terms of σ. It shows that their relative difference in σ is less than 10 −5 throughout the expansion history. The right panel displays the evolution of ∆N eff in both models. It confirms that whenever the backreaction is important, that is, during the RD era, the value of ∆N eff from the approximate model agrees with that from the exact model. In summary, the self-consistent expansion history for the exact model can be faithfully mimicked by that from the computationally-efficient, approximate model.