# Inflationary Gravitational Waves as a probe of the unknown post-inflationary
primordial Universe

Athul K. Soman  Swagat S. Mishra  Mohammed Shafi  Soumen Basak

###### Abstract

One of the key predictions of the standard inflationary paradigm is the
quantum mechanical generation of the transverse and traceless tensor
fluctuations due to the rapid accelerated expansion of space, which later
constitute a stochastic background of primordial gravitational waves (GWs).
The amplitude of the (nearly) scale-invariant inflationary tensor power
spectrum at large scales provides us with crucial information about the energy
scale of inflation in the case of the minimal inflaton coupling to gravity.
Furthermore, the spectral energy density, $\Omega_{{}_{\rm GW}}(f)$, of the
GWs at sufficiently small scales (or, large frequencies $f$) serves as an
important observational probe of post-inflationary primordial dynamics. In
fact, the small-scale spectral tilt, $n_{{}_{\rm GW}}=\frac{{\rm
d}\log{\Omega_{{}_{\rm GW}}}}{{\rm d}\log{f}}$, of the spectral energy density
of GWs is sensitive to the (unknown) post-inflationary equation of state
(EoS), $w$, of the universe; with a softer EoS ($w<1/3$) leading to a red
tilt: $n_{{}_{\rm GW}}<0$, while a stiffer EoS ($w>1/3$) resulting in a blue
tilt: $n_{{}_{\rm GW}}>0$. The post-inflationary dynamics, however, is
generically expected to be quite complex, potentially involving a number of
distinct phases. Hence, in this work, we discuss the possibility of multiple
sharp transitions, namely $w_{1}\to w_{2}\to w_{3}\to...\to w_{n}$, in the EoS
of the post-inflationary universe and compute the corresponding spectral
energy density of the inflationary GWs. We explicitly determine the region of
the parameter space $\\{{w_{1},\,w_{2},\,w_{3},...,w_{n}\\}}$ which leads to a
potentially detectable signal in the upcoming GW detectors, without violating
the current constraints.

##  1 Introduction

Cosmic inflation [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14] is currently
the leading paradigm in providing natural initial conditions for the hot Big
Bang phase [6, 15, 16, 17]. In the simplest scenario, inflation is assumed to
be sourced by a single (real) scalar field $\phi$, called the inflaton field,
which rolls slowly down its potential $V(\phi)$, and couples minimally to the
Ricci scalar [6, 8, 9, 14]. Of particular importance is the spectrum of
(nearly) scale-invariant scalar perturbations on super-Hubble scales,
generated via quantum fluctuations due to the rapid accelerated expansion of
space during inflation [18, 19, 20, 21, 6]. Upon their Hubble-entry, the
scalar curvature fluctuations induce the temperature and density
inhomogeneities in the primordial plasma, which later grow via gravitational
instability, leading to the formation of the large-scale structure (LSS) in
the universe [6, 16, 22, 23, 24]. The latest Cosmic Microwave Background (CMB)
observations [25, 26, 27, 28] provide strong support for the single field
slow-roll inflationary paradigm.

Another important aspect of the inflationary dynamics is the quantum
mechanical generation of a nearly scale-invariant spectrum of transverse and
traceless tensor fluctuations [29] on super-Hubble scales. The inflationary
tensor fluctuations, upon their Hubble-entry in the post-inflationary epochs,
constitute a stochastic background of primordial gravitational waves (GWs)
[30, 31, 32]. Both scalar and tensor primordial fluctuations leave distinct
imprints on the CMB, which can be parameterised in terms of two sets of
inflationary observables, namely, the amplitudes and the spectral indices of
the primordial power spectra. In particular, the parameter space comprising of
the scalar spectral index $n_{{}_{S}}-1$ and the tensor-to-scalar ratio $r$ is
quite illustrative in ruling out a large class of inflationary models [33, 8,
26, 27, 34, 35]. The latest CMB observations [26, 28] impose constraints on
the amplitude of primordial GWs, expressed as an upper bound on the tensor-to-
scalar ratio, namely, $r<0.036$ at $95\%$ confidence, which already disfavours
a large class of single field inflationary models, see Ref. [35]. Moreover,
the inflationary tensor fluctuations may potentially be important in
distinguishing between different competing frameworks of the pre-hot Big Bang
universe [36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52,
53, 54, 55, 56, 57] observationally.

Since these GWs originate during the pre-hot Big Bang universe and propagate
towards the present epoch, they encode important information about the
physical processes during inflation as well as the post-inflationary dynamics
of the universe. In particular, given that they are generated due to the
exponential expansion during inflation, the amplitude of the spectral energy
density $\Omega_{{}_{\rm GW}}(f)$ of the primordial GWs at large length scales
(or small frequencies $f$) is related to the energy scale (Hubble scale $H$)
of inflation [30, 31, 58, 34, 59]. Additionally, the spectral index of
$\Omega_{{}_{\rm GW}}(f)$, defined by $n_{{}_{\rm
GW}}=\frac{\mathrm{d}\log{\Omega_{{}_{\rm GW}}}}{\mathrm{d}\log{f}}$, is
sensitive to the equation of state (EoS) of the post-inflationary evolution
[31]. In particular, the spectral tilt $n_{{}_{\rm GW}}$ at small scales (or
high frequencies $f$) is determined by the (unknown) post-inflationary EoS $w$
of the primordial universe prior to the commencement of the Big Bang
Nucleosynthesis (BBN) [31, 34]. A softer EoS ($w<1/3$) leads to a red tilt:
$n_{{}_{\rm GW}}<0$, while a stiffer EoS ($w>1/3$) results in a blue tilt:
$n_{{}_{\rm GW}}>0$. A combination of the amplitude and tilt of
$\Omega_{{}_{\rm GW}}(f)$ across a broad range of frequency bands would then
enable us to probe the dynamics of the early universe at energy scales far
beyond the reach of the most powerful terrestrial particle accelerators [31,
58, 60, 61, 62, 63, 34, 59]. Hence the search for the primordial GWs is often
regarded as one of the key future targets in Cosmology.

A number of GW detectors are already operational at present, while more
sensitive and advanced future detectors, both ground-based and space-based, at
different frequency ranges are under construction. The existing (and upcoming)
GW detectors can be primarily classified into three categories, namely, (i)
CMB B-mode polarisation probes operating in the ultra-low frequency range
($10^{-18}\,{\rm Hz}\lesssim f\lesssim 10^{-16}\,{\rm Hz}$), (ii) pulsar
timing arrays (PTAs) sensitive to GWs in the low frequency range
($10^{-9}\,{\rm Hz}\lesssim f\lesssim 10^{-7}\,{\rm Hz}$), and (iii) atomic
and laser interferometers operating in the relatively high frequency range
$(10^{-5}\,{\rm Hz}\lesssim f\lesssim 10^{3}\,{\rm Hz})$ [64]. Amongst the
ultra-low frequency GW detectors, there are CMB probes such as the BICEP/Keck
Array [65], LiteBIRD [66], and POLARBEAR/Simons Array [67], which are designed
to measure the B-mode polarisation of the CMB photons, induced by the
quadrupolar primordial GWs [68, 69]. The recent data release by the PTA
observations, such as the NANOGrav collaboration [70, 71, 72], the EPTA/InPTA
[73, 74, 75] and IPTA [76] provide tentative evidence for the detection of a
background of stochastic GWs at low frequencies ($\sim{\cal
O}\left(10^{1}\right)\,{\rm nHz}$ range), whose primary origin, whether
astrophysical or primordial, is yet to be established [77, 78]. The currently
operational interferometers such as the advanced Laser Interferometer
Gravitational-wave Observatory (aLIGO) [79, 80], Virgo [81], KAGRA [82] and
the upcoming/proposed interferometers, such as the Laser Interferometer Space
Antenna (LISA) [83, 84, 85], Cosmic Explorer (CE) [86], Einstein Telescope
(ET) [87], Big Bang Observer (BBO) [88], and the DECi-hertz Interferometer
Gravitational wave Observatory (DECIGO) [89, 90], will be able to probe the
(relatively) large frequency GWs from the early universe associated with the
small-scale inflationary dynamics, as well as the post-inflationary expansion
history of the universe, prior to the epoch of Quantum Chromodynamics (QCD)
phase transition.

During the period of single field slow-roll inflation, the energy density of
the universe is primarily comprised of the homogeneous inflaton condensate
$\phi$, which rolls down its potential $V(\phi)$ in a slow terminal speed,
thereby driving the near-exponential expansion of space. If the inflaton
couplings to other external fields are low enough, then particle production
during inflation can be neglected, since the rapid accelerated expansion of
space quickly dilutes away the decay products. After the end of inflation,
however, the inflaton field begins to oscillate around the minimum of its
potential [91] and begins to transfer its energy to the particles of external
offspring fields that it is coupled to via parametric resonance [92, 93, 94,
95, 96], and/or to the inflaton inhomogeneities $\delta\varphi$ via self-
resonance (when attractive self-interactions are present in $V(\phi)$) [97,
98, 96, 99, 100]. In certain classes of inflation, where the inflaton
potential does not feature a stable minimum, rather it decays to zero at large
field values, particle production via resonance is absent. However, the
inflaton can still decay either perturbatively [101] or via the process of
instant preheating [102]. Nevertheless, the offspring fields (and other decay
products) undergo further interaction, scattering and thermalization; as a
result, the universe eventually transitions to the radiation-dominated thermal
plasma phase (or, the hot Big bang phase) before the commencement of BBN. This
transient phase between the end of inflation and the beginning of the thermal
radiation domination is referred to as the epoch of reheating, which is
supposed to be the origin of all primordial matter and energy in the universe
[94, 96, 14].

A number of comprehensive studies conduced in the past three decades [92, 94,
93, 95, 96, 103, 104, 105, 106, 107, 108, 109] demonstrate that the physical
processes taking place during reheating is highly complex, and potentially
non-linear. However, in most cases, the reheating dynamics can be broadly
divided into three stages: (i) preheating, (ii) backreaction, and (iii)
thermalization [92, 94, 110, 14]. The initial stage of preheating usually
exhibits rapid and efficient particle production via broad parametric
resonance [92, 94, 109]. Eventually, the backreaction of the produced
particles leads to the complete/partial fragmentation of the coherent inflaton
condensate. Additionally, cosmological redshifting drives the momenta of
produced particles away from the broad resonance band, thereby shutting the
resonance down and quenching the resonant particle production. This phase is
succeeded by a long and gradual perturbative decay of the inflaton, as well as
the re-scattering of decay products. Hence, the final stage of thermalization
is usually the longest of the three stages, which eventually leads to the hot
Big Bang phase, as discussed before. If strong attractive self-interaction of
the inflaton field is present, then it can lead to the formation of quasi-
solitonic structures like oscillons, resulting in the production of GWs [111,
112, 113, 97, 98, 96, 99, 100]. Depending upon the type of inflaton potential
and the inflaton couplings to external fields, the duration and EoS of each of
these phases can be different.

Despite the profusion of theoretical and numerical work in the field, the
epoch of reheating remains observationally inaccessible at present. In
general, there exist two distinct types of observational probes of the early
universe physics, namely (1) primordial relics such as ultra long-lived
solitons (e.g. strings, domain walls, oscillons), and (2) propagating
primordial messengers (signals) such as primordial GWs. Furthermore, these
primordial GWs may either originate from the physical processes during
reheating or from an earlier epoch, such as inflation. In this work, we focus
on the GW signal which gets generated during inflation and propagates towards
us, encoding the unknown primordial dynamics of reheating after inflation. In
particular, the duration of and the equation of state (EoS) during reheating
can be inferred from the spectral energy density of primordial GWs at higher
frequencies [34, 110, 59]. Furthermore, the GW spectrum is also influenced by
a number of other factors, such as the anisotropic stress from the free-
streaming relativistic neutrinos [114] and the variation in the effective
number of relativistic degrees of freedom throughout the history of the
universe [115, 116, 117, 116]. Since our interest lies with in the gross
features of the GW spectrum, we will not consider these factors in our
computation.

There has been a number of interesting papers in the recent literature [118,
119, 34, 120, 121, 59, 122, 123, 124, 125, 126, 127, 128, 129] investigating
the possibility of probing the epoch of reheating via inflationary GWs, where
the kinematics during reheating is parameterised by the expansion rate during
($\sim$ at the end of) inflation $\,H_{\rm inf}$, the temperature achieved at
the end of reheating $T_{\rm r*}$, the duration and the EoS $w_{{}_{\rm re}}$
during reheating. However, a majority of the papers has primarily assumed a
single (average) EoS, $w_{{}_{\rm re}}$, throughout the reheating history111
However, see refs. [121, 59] for earlier work involving multiple primordial
equation of states.. As mentioned earlier, the dynamics of reheating features
a number of complex non-linear phases, and hence, the universe is likely to
have undergone multiple transitions through different equations of state
during reheating [130, 104]. Additionally, the universe might have undergone
multiple phase-transitions [131, 132, 133, 134] after the completion of
reheating in the radiation dominated epoch, leading to a significant deviation
of the EoS from $w=1/3$.

Consequently, in this paper, we carry out a thorough investigation of the
spectrum of primordial stochastic GWs by considering the unknown post-
inflationary history of the universe to have undergone multiple phases of
piece-wise (nearly) constant EoS, namely, $w_{1}\to w_{2}\to w_{3}\to...\to
w_{n}$. For the sake of simplicity in calculation, we further assume that the
transition between any two successive epochs is sharp enough to be modelled as
an instantaneous transition. Under the aforementioned assumptions, we first
solve the evolution equation for the Fourier mode functions of the tensor
fluctuations in different epochs of constant $w_{i}$, and then use the Israel
junction matching conditions [135, 31, 136] to determine the full tensor mode
functions. Using the thus obtained analytical expression for the mode
functions, we finally compute the present-day spectral energy density of the
(first-order) inflationary GWs. We highlight the subspace of the parameter
space of EoS $\\{w_{1},\,w_{2},\,...\,w_{n}\\}$, and the corresponding
duration of each epoch of constant $w_{i}$, satisfying two conditions: (i)
yielding a potentially detectable background of stochastic GWs in the upcoming
detectors, such as LISA, BBO, DECIGO, CE, and ET; while at the same time, (ii)
satisfying the constraints imposed by the already existing GW, CMB, and BBN
observations. To be specific, we provide forecasts for GWs corresponding to
the unknown primordial history of the universe comprising of (up to) three
different piece-wise constant EoS phases $\\{w_{1},\,w_{2},\,w_{3}\\}$. As a
concrete example, we apply our method to determine the parameter space of the
duration of different piece-wise constant $w_{i}$ epochs before the
commencement of BBN in the framework of the recently proposed String Theory
inspired model of the post-inflationary universe in Ref. [137], which leads to
a detectable GW background in the aforementioned detectors.

Our paper is organised as follows: we begin with a brief review of the single
field slow-roll inflationary dynamics in Sec. 2, including the spectra of
scalar and tensor quantum fluctuations generated during inflation. We also
provide a succinct discussion about the latest CMB constraints on the
inflationary observables, which serve as the initial conditions for the post-
inflationary evolution. In Sec. 3, we first derive analytical expressions for
the tensor mode functions in the post-inflationary epochs using junction
matching conditions, and then move on to compute the corresponding spectral
energy density of GWs. Sec. 4 is dedicated to determining the parameter space
of the multiple equations of states (and their corresponding duration) which
leads to a potentially detectable signal in the upcoming GW detectors. In Sec.
5, we focus on the post-inflationary evolution of the universe in a specific
(String-inspired) model of the primordial universe, as proposed in Ref. [137].
Finally, we spell out the primary conclusions from our results, and provide
further discussions on our work in Sec. 6. Various appendices provide
supplementary material to the analysis carried out in the main sections.

We work with natural units $\hbar=c=1$ throughout this paper, and define the
reduced Planck mass to be $m_{p}=1/\sqrt{8\pi G}=2.44\times 10^{18}$ GeV. The
background universe is considered to be described by the flat Friedmann-
Lemaitre-Robertson-Walker (FLRW) line element, with metric signature
$(-,\,+,\,+,\,+)$. Derivative with respect to cosmic time ‘$t$’ is denoted
with an overdot $\left(\;\dot{}\;\right)$, while derivative with respect to
conformal time ‘$\tau$’ is denoted with an overprime
$\left(\;{}^{\prime}\;\right)$.

##  2 Inflationary dynamics

We work in the simplest scenario where inflation is sourced by a single
canonical scalar field (inflaton) $\varphi(t,\,\vec{x})$, whose homogeneous
part is denoted as $\phi(t)$, i.e.,
$\varphi(t,\,\vec{x})=\phi(t)+\delta\varphi(t,\,\vec{x})$. The action
governing the dynamics of the scalar field $\varphi$ is given by [6, 9, 14]

$S[g_{\mu\nu},\,\varphi]=\int\,\mathrm{d}^{4}x\,\sqrt{-g}\,\left(\frac{1}{2}\,m_{p}^{2}\,R-\frac{1}{2}\,g^{\mu\nu}\,\partial_{\mu}\varphi\,\partial_{\nu}\varphi\,-V(\varphi)\right)\,,$
(2.1)

where $g_{\mu\nu}$ is the metric tensor, $R$ is the Ricci scalar curvature and
$V(\varphi)$ is the inflaton potential. The background spacetime is described
by the flat FLRW line element

$\mathrm{d}s^{2}=-\mathrm{d}t^{2}+a^{2}(t)\left[\mathrm{d}x^{2}+\mathrm{d}y^{2}+\mathrm{d}z^{2}\right]\,,$
(2.2)

where $a(t)$ is the scale factor whose time evolution is governed by the
Friedmann equation

$H^{2}=\frac{1}{3\,m_{p}^{2}}\,\rho_{\phi}\equiv\frac{1}{3\,m_{p}^{2}}\,\left[\frac{1}{2}\,\dot{\phi}^{2}+V(\phi)\right]\,,$
(2.3)

where $H=\dot{a}/a$ is the Hubble parameter. The homogeneous inflaton
condensate satisfies

$\ddot{\phi}+3\,H\,\dot{\phi}+\frac{\mathrm{d}V}{\mathrm{d}\phi}=0\,.$ (2.4)

Departure from a pure de Sitter (exponential) expansion during inflation is
characterised by the first slow-roll parameter

$\epsilon_{{}_{H}}=-\frac{\mathrm{d}\ln
H}{\mathrm{d}N}=-\frac{\dot{H}}{H^{2}}=\frac{1}{2\,m_{p}^{2}}\,\frac{\dot{\phi}^{2}}{H^{2}}\,,$
(2.5)

while a small enough value of the second slow-roll parameter, defined by

$\eta_{{}_{H}}=-\frac{\ddot{\phi}}{H\,\dot{\phi}}=\epsilon_{{}_{H}}-\frac{1}{2}\,\frac{\mathrm{d}\ln\epsilon_{{}_{H}}}{\mathrm{d}N}\,,$
(2.6)

ensures that inflation lasts for long enough. Note that $N(a)=\ln(a/a_{i})$
marks the number of $e$-folds of expansion during inflation, with $a_{i}$
being the scale factor at an arbitrary initial time during inflation [9].
Slow-roll inflation is defined by

$\epsilon_{{}_{H}}\ll 1\,,\quad{\rm and}\quad|\eta_{{}_{H}}|\ll 1\,.$ (2.7)

###  2.1 Quantum fluctuations during Inflation

The rapid accelerated expansion during inflation stretches the small-scale
(sub-Hubble) vacuum fluctuations to large super-Hubble scales, which then
remain frozen until their subsequent Hubble entry in the post-inflationary
epochs, as illustrated in Fig. 1. In fact, the simplest scenario of single
field slow-roll inflation generates both scalar- and tensor-type fluctuations,
see Refs. [9, 14]; the former induce temperature and density fluctuations in
the primordial plasma, while the latter lead to a stochastic GW background, as
discussed in Sec. 3. Since CMB observations suggest that the primordial
fluctuations are small, their dynamics is well described by the framework of
linear (cosmological) perturbation theory, where the scalar and the tensor
perturbations are decoupled (SVT decomposition theorem). Thus, they can be
treated independently of each other [9].

####  2.1.1 Scalar fluctuations during inflation

Scalar fluctuations during inflation are described by the gauge-invariant
comoving curvature perturbation $\zeta(\tau,\,\vec{x})$, whose action, in the
comoving gauge ($\delta\varphi=0$), takes the form [138, 9]

$S[\zeta(\tau,\vec{x})]=\frac{1}{2}\int\mathrm{d}\tau\,\mathrm{d}^{3}\vec{x}\,\left(\frac{a\,\phi^{\prime}}{\cal
H}\right)^{2}\left[(\zeta^{\prime})^{2}-(\gradient\zeta)^{2}\right]\,,$ (2.8)

where ${\cal H}=a^{\prime}/a$, is the conformal Hubble rate. To transform the
action into its Minkowski equivalent, we define the Mukhanov-Sasaki (MS)
variable $v=(a\,\phi^{\prime}/{\cal H})\,\zeta$. For simplification, we denote
$z=a\,\phi^{\prime}/{\cal H}=\sqrt{2\,\epsilon_{{}_{H}}}\,a\,m_{p}$, where Eq.
(2.5) has been used to get the second equality. The equation of motion in
terms of the MS variable in Fourier space is given by

$v^{\prime\prime}_{k}\,+\,\left(\,k^{2}-\frac{z^{\prime\prime}}{z}\,\right)v_{k}=0\,,$
(2.9)

where

$\frac{z^{\prime\prime}}{z}={\cal
H}^{2}\left[2+2\epsilon_{{}_{H}}-3\eta_{{}_{H}}+2\epsilon_{{}_{H}}^{2}+\eta_{{}_{H}}^{2}-3\epsilon_{{}_{H}}\eta_{{}_{H}}-\frac{\eta_{{}_{H}}^{\prime}}{\cal
H}\,\right]\,.$ (2.10)

To solve the MS equation analytically, we define a new parameter $\nu$ such
that $z^{\prime\prime}/z={\cal H}^{2}(\nu^{2}-1/4)$. For slow-roll inflation
($\epsilon_{{}_{H}},\,|\eta_{{}_{H}}|,\,|\eta^{\prime}_{H}|\ll 1)$, we can
treat $\nu$ to be (approximately) a constant. The MS Eq. (2.9) can be written
in terms of $\nu$ as

$v^{\prime\prime}_{k}\,+\,\left(\,k^{2}-\frac{\nu^{2}-1/4}{\tau^{2}}\,\right)v_{k}=0\,,$
(2.11)

where we have substituted ${\cal H}=-\tau^{-1}$ due to the quasi-de Sitter
expansion during inflation. Eq. (2.11) is similar to the dynamic equation of a
simple harmonic oscillator with a time-dependent frequency. Imposing the
Bunch-Davies initial condition [139, 140], $v_{k}\xrightarrow[]{-k\tau\gg
1}e^{-ik\tau}/\sqrt{2k}$, results in a unique solution [136]

$v_{k}(\tau)=\frac{1}{\sqrt{2k}}\sqrt{\frac{\pi}{2}}\,e^{i(\nu+\frac{1}{2})\frac{\pi}{2}}\sqrt{-k\tau}\,H_{\nu}^{(1)}(-k\tau)\,,$
(2.12)

where $H_{\nu}^{(1)}$ is the Hankel function of first kind [141]. The
derivation of Eq. (2.12) is given in App. A. For slow-roll inflation, $\nu$
can be approximated as $\nu\approx 3/2$. Substituting this into Eq. (2.12)
leads to the following expression for the mode function

$\displaystyle
v_{k}(\tau)=\frac{1}{\sqrt{2k}}\left(1-\frac{i}{k\tau}\right)e^{-ik\tau}\,.$
(2.13)

The power spectrum for scalar fluctuations is defined as

${\cal
P}_{\zeta}(k)=\frac{k^{3}}{2\pi^{2}}\left|\zeta_{k}\right|^{2}=\frac{k^{3}}{2\pi^{2}}\left|\frac{v_{k}}{z}\right|^{2}\,.$
(2.14)

Substituting $z=am_{p}\sqrt{2\epsilon_{{}_{H}}}$, $a=-1/(H\tau)$ and taking
the super-Hubble limit $-k\tau\to 0$, we obtain the scalar power spectrum
during inflation to be

${\cal
P}_{\zeta}(k)=\frac{1}{8\pi^{2}\epsilon_{{}_{H}}}\left(\frac{H}{m_{p}}\right)^{2}\,.$
(2.15)

Figure 1: The figure illustrates the variation of the comoving Hubble radius
$(aH)^{-1}$ with the scale factor $a(t)$ (both in the logarithmic scale).
During inflation, the comoving Hubble radius decreases, causing the
fluctuations to exit the comoving Hubble radius. On super-Hubble scales, the
amplitude of the fluctuations remains frozen. After the end of inflation,
these fluctuations enter the comoving Hubble radius, which begins to increase
due to the decelerating expansion of space. Modes with larger wavelengths exit
the comoving Hubble radius earlier during inflation, as compared to the ones
with smaller wavelengths. The CMB pivot scale (corresponding to a wavenumber
$k_{*}=0.05\,{\rm Mpc}^{-1}$), shown in blue dashed line, makes its Hubble-
exit about $N_{k}$ $e$-folds before the end of inflation, while it enters the
comoving Hubble radius prior to the matter-radiation equality. The black dot
represents the end of reheating. The figure has been plotted assuming a
matter-like reheating equation of state, i.e. $w_{{}_{\rm re}}=0$, for the
purpose of illustration.

####  2.1.2 Tensor fluctuations during inflation

The action for the gauge-invariant, transverse and traceless tensor
fluctuations $h_{ij}$ is given by [138, 10]

$S[h_{ij}]=\frac{m_{p}^{2}}{8}\int\mathrm{d}\tau\leavevmode\nobreak\
\mathrm{d}^{3}\vec{x}\leavevmode\nobreak\
a^{2}\left[(h_{ij}^{\prime})^{2}-(\gradient h_{ij})^{2}\right]\,.$ (2.16)

Using the rotational symmetry of the background metric, we align the z-axis
along the direction of the momentum of the tensor mode, i.e.,
$\vec{k}=(0,\,0,\,k)$, and upon decomposing the tensor perturbation into its
two polarization components, we obtain

$\displaystyle h_{ij}$
$\displaystyle=\frac{1}{\sqrt{2}}\begin{pmatrix}h^{+}&\leavevmode\nobreak\
\leavevmode\nobreak\ h^{\times}&\leavevmode\nobreak\ \leavevmode\nobreak\ 0\\\
h^{\times}&-h^{+}&\leavevmode\nobreak\ \leavevmode\nobreak\ 0\\\
0\leavevmode\nobreak\ &\leavevmode\nobreak\ 0&\leavevmode\nobreak\
\leavevmode\nobreak\
0\end{pmatrix}=\epsilon^{+}_{ij}\,h^{+}+\,\epsilon^{\times}_{ij}\,h^{\times}\,,$
(2.17)

where $+,\,\times$ stands for the two polarisation states of GWs and

$\epsilon^{+}_{ij}=\frac{1}{\sqrt{2}}\begin{pmatrix}1&\leavevmode\nobreak\
\leavevmode\nobreak\ 0&\leavevmode\nobreak\ \leavevmode\nobreak\ 0\\\
0&-1&\leavevmode\nobreak\ \leavevmode\nobreak\ 0\\\ 0&\leavevmode\nobreak\
\leavevmode\nobreak\ 0&\leavevmode\nobreak\ \leavevmode\nobreak\
0\end{pmatrix};\quad\epsilon^{\times}_{ij}=\frac{1}{\sqrt{2}}\begin{pmatrix}0&\leavevmode\nobreak\
\leavevmode\nobreak\ 1&\leavevmode\nobreak\ \leavevmode\nobreak\ 0\\\
1&\leavevmode\nobreak\ \leavevmode\nobreak\ 0&\leavevmode\nobreak\
\leavevmode\nobreak\ 0\\\ 0&\leavevmode\nobreak\ \leavevmode\nobreak\
0&\leavevmode\nobreak\ \leavevmode\nobreak\ 0\end{pmatrix}\,.$ (2.18)

Hence, the action can be rewritten as [14]

$S[h^{+},\,h^{\times}]=\frac{1}{2}\int\mathrm{d}\tau\leavevmode\nobreak\
\mathrm{d}^{3}\vec{x}\,\left(\frac{a\,m_{p}}{2}\right)^{2}\,\sum_{\lambda=+,\times}\,\left[\left(h^{{}^{\prime}\lambda}\right)^{2}-\left(\gradient
h^{\lambda}\right)^{2}\right]\,,$ (2.19)

which is analogous to the action of two massless scalar fields [14] as given
in Eq. (2.8). For each polarisation mode of the tensor fluctuations, we define
a corresponding Mukhanov-Sasaki variable
$v^{\lambda}=\left(a\,m_{p}/2\right)h^{\lambda}$. The full tensor power
spectrum, consisting of both the aforementioned polarisation modes, is given
by

$\displaystyle{\cal P}_{T}(k)$
$\displaystyle=\frac{k^{3}}{2\pi^{2}}\left(\absolutevalue{h_{k}^{+}}^{2}+\absolutevalue{h_{k}^{\times}}^{2}\right)=\frac{k^{3}}{2\pi^{2}}\left(\frac{2}{a\,m_{p}}\right)^{2}\left(\absolutevalue{v_{k}^{+}}^{2}+\absolutevalue{v_{k}^{\times}}^{2}\right)\,.$
(2.20)

Substituting the solution for each of the polarisation mode in Eq. (2.20), as
given in Eq. (2.13), the inflationary tensor power spectrum in the super-
Hubble limit $k\ll aH$ becomes

${\cal P}_{T}(k)=\frac{2}{\pi^{2}}\left(\frac{H}{m_{p}}\right)^{2}\,.$ (2.21)

To facilitate the comparison between the inflationary predictions and the
observational data, it is instructive to express the scalar and tensor power
spectra as power laws around a pivot scale $k_{*}$, such as

$\displaystyle{\cal P}_{T}(k)$
$\displaystyle=A_{{}_{T}}\left(\frac{k}{k_{*}}\right)^{n_{{}_{T}}}\,;$
$\displaystyle{\cal P}_{\zeta}(k)$
$\displaystyle=A_{{}_{S}}\left(\frac{k}{k_{*}}\right)^{n_{{}_{S}}-1}\,,$
(2.22)

where $A_{{}_{S}}$ and $A_{{}_{T}}$ are the amplitudes of the scalar and
tensor power spectra, respectively, at the CMB pivot scale $k_{*}=0.05\,\rm
Mpc^{-1}$; while $n_{{}_{S}}-1$ and $n_{{}_{T}}$ are the scalar and tensor
spectral indices, given by

$n_{{}_{S}}-1=2\,\eta_{{}_{H}}-4\,\epsilon_{{}_{H}}\,;\quad
n_{{}_{T}}=-2\,\epsilon_{{}_{H}}\,.$ (2.23)

The tensor-to-scalar ratio is defined as,

$r=\frac{A_{{}_{T}}}{A_{{}_{S}}}\,.$ (2.24)

From Eqs. (2.15) and (2.21), the tensor-to-scalar ratio can be related to the
Hubble slow-roll parameter as

$r=16\,\epsilon_{{}_{H}}\,,$ (2.25)

which, using Eq. (2.23), leads to the single field consistency relation

$r=-8\,n_{{}_{T}}\,.$ (2.26)

###  2.2 Implications of the observational constraints

Latest observations of the anisotropies of CMB temperature and polarisation
fluctuations, combined with an absence of detection of the primordial tensor
fluctuations on large cosmological scales, allow us to put stringent
constraints on the inflationary power spectra in the single field slow-roll
paradigm. In particular, the primordial scalar fluctuations are observed to be
predominantly adiabatic, highly Gaussian, and nearly scale-invariant on CMB
scales, which is in accordance with the predictions of the standard single
field slow-roll models [142]. The amplitude of the scalar power spectrum at
CMB scales is obtained from the Planck mission to be $A_{{}_{S}}\simeq
2.1\times 10^{-9}$ [25, 35]. The scalar spectral index is constrained to be
$n_{{}_{S}}=0.9649\pm 0.0012$ at $68\%$ confidence [27]. The improved
constraints from the BICEP/Keck Array [28], combined with the Planck 2018
observations, lead to an upper bound on the tensor-to-scalar ratio $r<0.036$
at the CMB pivot scale $k_{*}=0.05\,\rm Mpc^{-1}$. This further constrains the
amplitude of the tensor fluctuations to be $A_{{}_{T}}<7.56\times 10^{-11}$ on
CMB scales.

The corresponding upper bounds on the expansion rate and energy scale during
inflation can be found using Eq. (2.21). The Hubble expansion rate during
(slow-roll) inflation is bounded by $H_{\rm inf}<4.64\times 10^{13}\
\rm{GeV}$. Using the Friedmann equation $H^{2}_{\rm inf}\simeq\rho_{\rm
inf}/{3m_{p}^{2}}$, the bound on $E_{\rm inf}$ is given by

$E_{\rm inf}=(\rho_{\rm inf})^{1/4}\implies E_{\rm inf}<1.39\times 10^{16}\
\rm{GeV}\,,$ (2.27)

which is orders of magnitude below the Planck scale ($\sim 10^{18}\,{\rm
GeV}$). For slow-roll inflation, the constraint on the tensor-to-scalar ratio,
along with Eq. (2.25) and Eq. (2.26) leads to $\epsilon_{{}_{H}}<2.25\times
10^{-3},\,|\eta_{{}_{H}}|\simeq 0.02$ and $|n_{{}_{T}}|<4.5\times 10^{-3}$.
These constraints further reinforce the near-scale invariance of the
inflationary power spectra; in particular, the tensor power spectrum has a
negligible scale dependence [35].

##  3 Primordial gravitational waves in the post-inflationary universe

The rapid accelerated expansion during inflation leads to the generation of
almost scale-invariant primordial tensor fluctuations on super-Hubble scales,
as discussed in Sec. 2.1.2. These tensor fluctuations, upon making their
Hubble-entry in the post-inflationary (decelerating) epochs, then propagate as
waves and thus, they constitute a stochastic gravitational wave (GW)
background. In this section, we carry out a detailed study of the evolution of
these primordial tensor fluctuations in the post-inflationary universe. There
has been a number of important papers in the literature in this direction,
e.g. Refs. [119, 34, 120, 121, 59, 131], to quote a few. In particular, we
would like to emphasize the work carried out in Refs. [119] where the authors
determined the parameter space of the post-inflationary epoch which leads to a
detectable signal in the aLIGO and LISA GW observatories. To be specific,
their analysis primarily focused on the estimation of the post-inflationary
parameter space, such as the expansion rate during (at the end of) inflation
($H_{\rm inf}$), the duration of reheating (characterised by the temperature
achieved at the end of reheating, $T_{\rm r*}$), and the EoS during reheating
($w_{{}_{\rm re}})$.

It is important to note that the authors of Ref. [119] modelled the post-
inflationary reheating epoch by assuming a single EoS. However, as mentioned
earlier, since the dynamics of reheating is quite complex, the universe is
likely to have undergone multiple transitions through different equations of
state during reheating [130, 104]. Furthermore, the early universe might have
undergone a number of phase-transitions after the completion of reheating in
the radiation dominated epoch, resulting in a significant deviation of the EoS
from $w=1/3$. Therefore, it is important to extend the analysis of Ref. [119]
by allowing for the possibility of multiple phases of piece-wise constant EoS.

Before proceeding to determine the post-inflationary evolution of the tensor
fluctuations, let us lay out the primary assumptions under which we base our
analysis of multiple primordial EoS. Our technique is similar to the one
carried out in Ref. [136] for multiple transitions during inflation in the
context of primordial black hole formation.

  1. 1.

We assume the universe to have undergone multiple sharp (instantaneous)
transitions in the EoS, $w_{{}_{\rm re}}$, i.e. $w_{1}\rightarrow w_{2}$ at
$\tau=\tau_{1}$ ; $w_{2}\rightarrow w_{3}$ at $\tau=\tau_{2}$ ;
$w_{3}\rightarrow w_{4}$ at $\tau=\tau_{3}$ ; and so on, in the post-
inflationary epoch before the commencement of BBN. Hence, the EoS can be
written as

$w_{{}_{\rm
re}}(\tau)=w_{1}\,+\,\left(w_{2}-w_{1}\right)\Theta\left(\tau-\tau_{1}\right)\,+\,\left(w_{3}-w_{2}\right)\Theta\left(\tau-\tau_{2}\right)\,+\,\left(w_{4}-w_{3}\right)\Theta\left(\tau-\tau_{3}\right)\,+\,...\,,$
(3.1)

where $\Theta(\tau-\tau_{n})$ is the Heaviside step function.

  2. 2.

The inflationary tensor modes are frozen on the super-Hubble scales, with
amplitude $h_{k,\,{\rm inf}}^{\lambda}$, before their Hubble-entry. Hence they
serve as the initial conditions for the post-inflationary evolution as,

$h_{k}^{\lambda}(\tau_{k})=h_{k,\,{\rm inf}}^{\lambda}\quad;\quad
h_{k}^{{}^{\prime}\lambda}(\tau_{k})=0\,,$ (3.2)

where $\tau_{k}$ is the conformal time when the tensor mode of interest with
the (comoving) wavenumber $k$ enters the comoving Hubble radius.

  3. 3.

We match the solutions of the equation of motion of tensor fluctuations across
two adjacent epochs with different (constant) EoS at the time of the
instantaneous transition, $\tau=\tau_{n}$, by using the Israel Junction
conditions [135], i.e.

$\displaystyle h_{k,\,b}^{\lambda}(\tau_{n})$
$\displaystyle=h_{k,\,a}^{\lambda}(\tau_{n})\,,\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ (\text{Continuity}),$ (3.3)
$\displaystyle h_{k,\,b}^{{}^{\prime}\lambda}(\tau)\Big{|}_{\tau=\tau_{n}}$
$\displaystyle=h_{k,\,a}^{{}^{\prime}\lambda}(\tau)\Big{|}_{\tau=\tau_{n}}\,,\quad(\text{Differentiability}),$
(3.4)

where $h_{k,\,b}^{\lambda}$ and $h_{k,\,a}^{\lambda}$ correspond to the tensor
amplitude before and after the transition respectively.

###  3.1 Evolution of gravitational waves

First, we solve the equation of motion for the tensor mode functions in the
post-inflationary epoch. The action in Eq. (2.19) leads to the field equation

$h_{ij}^{{}^{\prime\prime}\lambda}+2\,\frac{a^{\prime}}{a}\,h_{ij}^{{}^{\prime}\lambda}-\vec{\nabla}^{2}\,h_{ij}^{\lambda}=0\,,$
(3.5)

and the corresponding equation for the evolution of Fourier mode functions
(for each polarisation) becomes

$h_{k}^{{}^{\prime\prime}\lambda}+2\,\frac{a^{\prime}}{a}\,h_{k}^{{}^{\prime}\lambda}+k^{2}\,h_{k}^{\lambda}=0\,.$
(3.6)

For an epoch with a constant EoS, say $w$, the conformal Hubble parameter can
be written as (see App. B)

${\cal
H}=\frac{a^{\prime}}{a}=aH=a_{i}H_{i}\left[1+\frac{a_{i}H_{i}(\tau-\tau_{i})}{\alpha}\right]^{-1}\,,$
(3.7)

where ‘$i$’ denotes the beginning of the epoch and the parameter
$\alpha=2/(1+3w)$. We introduce the variable $y(\tau)=k/[a(\tau)H(\tau)]$ as
the new independent variable in order to rewrite Eq. (3.6) in the desired form

$\frac{\mathrm{d}^{2}h_{k}^{\lambda}}{\mathrm{d}y^{2}}+2\,\frac{\alpha}{y}\,\frac{\mathrm{d}h_{k}^{\lambda}}{\mathrm{d}y}+\alpha^{2}\,h_{k}^{\lambda}=0\,.$
(3.8)

Note that $y<1$ for super-Hubble modes and $y>1$ for sub-Hubble modes.
Defining $f_{k}^{\lambda}(x)=x^{(\alpha-1/2)}\,h_{k}^{\lambda}$ with
$x=\alpha\,y$, we can convert Eq. (3.8) to the following expression

$\frac{\mathrm{d}^{2}f_{k}^{\lambda}}{\mathrm{d}x^{2}}+\frac{1}{x}\,\frac{\mathrm{d}f_{k}^{\lambda}}{\mathrm{d}x}+\left[1-\frac{\left(\alpha-\frac{1}{2}\right)^{2}}{x^{2}}\right]\,f_{k}^{\lambda}=0\,,$
(3.9)

which is of the form of a Bessel equation. Hence, the general solution of this
equation can be written in terms of a linear combination of two independent
Bessel functions of the first-kind, $J_{\left(\alpha-\frac{1}{2}\right)}(x)$
and $J_{-\left(\alpha-\frac{1}{2}\right)}(x)$ [141],

$\displaystyle h_{k}^{\lambda}(y)$
$\displaystyle=\frac{1}{(\alpha\,y)^{\alpha-\frac{1}{2}}}\left[A_{k}\,J_{\left(\alpha-\frac{1}{2}\right)}(\alpha\,y)+B_{k}\,J_{-\left(\alpha-\frac{1}{2}\right)}(\alpha\,y)\right]\;,$
(3.10) $\displaystyle h_{k}^{{}^{\prime}\lambda}(y)$
$\displaystyle=\frac{-k}{(\alpha\,y)^{\alpha-\frac{1}{2}}}\left[A_{k}\,J_{\left(\alpha+\frac{1}{2}\right)}(\alpha\,y)-B_{k}\,J_{-\left(\alpha+\frac{1}{2}\right)}(\alpha\,y)\right]\,.$
(3.11)

The coefficients $A_{k}$ and $B_{k}$ in the solution are obtained by junction
matching conditions at the transitions of two adjacent epochs as specified in
Eqs. (3.3) and (3.4). Specifically, consider a sharp transition from the
$(n-1)^{\rm th}$ epoch to the $n^{\rm th}$ epoch with the equations of state
$w_{(n-1)}$ and $w_{n}$ respectively, at time $\tau=\tau_{(n-1)}$. The
expression for the corresponding tensor amplitude can be written as

$h_{k}^{\lambda}(y)=\begin{cases}\cfrac{\left[A_{k,\,m}\,J_{\left(\alpha_{m}-\frac{1}{2}\right)}(\alpha_{m}\,y)+B_{k,\,m}\,J_{-\left(\alpha_{m}-\frac{1}{2}\right)}(\alpha_{m}\,y)\right]}{(\alpha_{m}\,y)^{\alpha_{m}-\frac{1}{2}}},&y_{(m-1)}\leq
y\leq y_{m}\,,\\\ \\\
\cfrac{\left[A_{k,\,n}\,J_{\left(\alpha_{n}-\frac{1}{2}\right)}(\alpha_{n}\,y)+B_{k,\,n}\,J_{-\left(\alpha_{n}-\frac{1}{2}\right)}(\alpha_{n}\,y)\right]}{(\alpha_{n}\,y)^{\alpha_{n}-\frac{1}{2}}},&y_{m}\leq
y\leq y_{n}\,,\end{cases}$ (3.12)

where $m=n-1$ and $y_{n}=k/[a(\tau_{n})H(\tau_{n})]$. Similarly the
derivatives of the tensor amplitude is given by

$h_{k}^{{}^{\prime}\lambda}(y)=\begin{cases}-k\,\cfrac{\left[A_{k,\,m}\,J_{\left(\alpha_{m}+\frac{1}{2}\right)}(\alpha_{m}\,y)-B_{k,\,m}\,J_{-\left(\alpha_{m}+\frac{1}{2}\right)}(\alpha_{m}\,y)\right]}{(\alpha_{m}\,y)^{\alpha_{m}-\frac{1}{2}}},&y_{(m-1)}\leq
y\leq y_{m}\,,\\\ \\\
-k\,\cfrac{\left[A_{k,\,n}\,J_{\left(\alpha_{n}+\frac{1}{2}\right)}(\alpha_{n}\,y)-B_{k,\,n}\,J_{-\left(\alpha_{n}+\frac{1}{2}\right)}(\alpha_{n}\,y)\right]}{(\alpha_{n}\,y)^{\alpha_{n}-\frac{1}{2}}},&y_{m}\leq
y\leq y_{n}\,.\end{cases}$ (3.13)

Imposing junction matching conditions (3.3) and (3.4) on Eqs. (3.12) and
(3.13), we obtain

$\displaystyle\frac{A_{k,\,n}\,g_{1}+B_{k,\,n}\,f_{1}}{\left(\alpha_{n}\,y_{m}\right)^{\left(\alpha_{n}-\frac{1}{2}\right)}}=\frac{A_{k,\,m}\,g_{2}+B_{k,\,m}\,f_{2}}{\left(\alpha_{m}\,y_{m}\right)^{\left(\alpha_{m}-\frac{1}{2}\right)}},$
(3.14)
$\displaystyle\frac{A_{k,\,n}\,g_{3}-B_{k,\,n}\,f_{3}}{\left(\alpha_{n}\,y_{m}\right)^{\left(\alpha_{n}-\frac{1}{2}\right)}}=\frac{A_{k,\,m}\,g_{4}-B_{k,\,m}\,f_{4}}{\left(\alpha_{m}\,y_{m}\right)^{\left(\alpha_{m}-\frac{1}{2}\right)}},$
(3.15)

where

$\displaystyle g_{1}$
$\displaystyle=J_{\left(\alpha_{n}-\frac{1}{2}\right)}(\alpha_{n}\,y_{m}),$
$\displaystyle f_{1}$
$\displaystyle=J_{-\left(\alpha_{n}-\frac{1}{2}\right)}(\alpha_{n}\,y_{m})\,,$
(3.16) $\displaystyle g_{2}$
$\displaystyle=J_{\left(\alpha_{m}-\frac{1}{2}\right)}(\alpha_{m}\,y_{m}),$
$\displaystyle f_{2}$
$\displaystyle=J_{-\left(\alpha_{m}-\frac{1}{2}\right)}(\alpha_{m}\,y_{m})\,,$
(3.17) $\displaystyle g_{3}$
$\displaystyle=J_{\left(\alpha_{n}+\frac{1}{2}\right)}(\alpha_{n}\,y_{m}),$
$\displaystyle f_{3}$
$\displaystyle=J_{-\left(\alpha_{n}+\frac{1}{2}\right)}(\alpha_{n}\,y_{m})\,,$
(3.18) $\displaystyle g_{4}$
$\displaystyle=J_{\left(\alpha_{m}+\frac{1}{2}\right)}(\alpha_{m}\,y_{m}),$
$\displaystyle f_{4}$
$\displaystyle=J_{-\left(\alpha_{m}+\frac{1}{2}\right)}(\alpha_{m}\,y_{m})\,\,.$
(3.19)

A slight rearrangement of the terms yields the following relations between the
coefficients before and after transition

$\displaystyle
A_{k,\,n}=\frac{\left(\alpha_{n}\,y_{m}\right)^{\left(\alpha_{n}-\frac{1}{2}\right)}}{\left(\alpha_{m}\,y_{m}\right)^{\left(\alpha_{m}-\frac{1}{2}\right)}}\,\frac{\left[A_{k,\,m}\left(g_{2}\,f_{3}+g_{4}\,f_{1}\right)+B_{k,\,m}\left(f_{2}\,f_{3}-f_{4}\,f_{1}\right)\right]}{f_{1}\,g_{3}+g_{1}\,f_{3}}\,,$
(3.20) $\displaystyle
B_{k,\,n}=\frac{\left(\alpha_{n}\,y_{m}\right)^{\left(\alpha_{n}-\frac{1}{2}\right)}}{\left(\alpha_{m}\,y_{m}\right)^{\left(\alpha_{m}-\frac{1}{2}\right)}}\,\frac{\left[A_{k,\,m}\left(g_{2}\,g_{3}-g_{4}\,g_{1}\right)+B_{k,\,m}\left(f_{2}\,g_{3}+f_{4}\,g_{1}\right)\right]}{f_{1}\,g_{3}+g_{1}\,f_{3}}\,.$
(3.21)

Utilizing Eqs. (3.20) and (3.21), we can find the coefficients $A_{k}$ and
$B_{k}$ for any epoch in the post-inflationary universe, except for the first
epoch right after the end of inflation as they have to be calculated using Eq.
(3.2). We stress that the expression for mode functions given in Eq. (3.10) is
applicable at all scales in the post-inflationary universe, irrespective of
whether the mode functions are sub-Hubble222Note that we do not consider the
deep-UV modes that never exited the Hubble-radius during inflation. In this
sense, we have a UV cutoff around $k_{\rm UV}=k_{\rm end}$. For a thorough
discussion on regularization and renormalization of tensor modes, see Refs.
[143, 144]. or super-Hubble, with the junction matching conditions having been
imposed at the time of sharp transition.

####  3.1.1 Universe with a single equation of state during reheating

We begin by providing expressions for the coefficients of the tensor mode
functions in Eq. (3.10) for the simplest case, where the post-inflationary
reheating epoch is characterized by a single constant EoS.

  * •

Reheating epoch

Assume that the reheating epoch is characterised by a single effective EoS,
$w_{{}_{\rm re}}$. Then, the resultant tensor amplitude Eq. (3.10) for
reheating epoch can be written as

$h_{k,\,{\rm re}}^{\lambda}(y)=\frac{1}{(\alpha_{\rm re}\,y)^{\alpha_{\rm
re}-\frac{1}{2}}}\left[A_{k,\,{\rm re}}\,J_{\left(\alpha_{\rm
re}-\frac{1}{2}\right)}(\alpha_{\rm re}\,y)+B_{k,\,{\rm
re}}\,J_{-\left(\alpha_{\rm re}-\frac{1}{2}\right)}(\alpha_{\rm
re}\,y)\right]\,,$ (3.22)

where $\alpha_{\rm re}=2/(1+3\,w_{\rm re})$. In order to use the initial
conditions, we need the super-Hubble limits of the tensor mode and its
derivative, for which we make use of the small argument limit of the Bessel
function, $J_{\nu}(y\ll 1)\approx(y/2)^{\nu}/\Gamma(\nu+1)$, in the super-
Hubble regime. The tensor mode function and its derivative on super-Hubble
scales take the form

$\displaystyle\lim_{y\rightarrow 0}\,h_{k,\,{\rm re}}^{\lambda}(y)$
$\displaystyle=A_{k,\,{\rm re}}\,\frac{1}{2^{\left(\alpha_{\rm
re}-\frac{1}{2}\right)}}\frac{1}{\Gamma\left(\alpha_{\rm
re}+\frac{1}{2}\right)}+B_{k,\,{\rm re}}\,\frac{2^{\left(\alpha_{\rm
re}-\frac{1}{2}\right)}}{\left(\alpha_{\rm re}\,y\right)^{2\left(\alpha_{\rm
re}-\frac{1}{2}\right)}}\frac{1}{\Gamma\left(-\alpha_{\rm
re}+\frac{3}{2}\right)}\,,$ (3.23) $\displaystyle\lim_{y\rightarrow
0}\,h_{k,\,{\rm re}}^{{}^{\prime}\lambda}(y)$
$\displaystyle=-k\,\left[A_{k,\,{\rm re}}\,\frac{\alpha_{\rm
re}\,y}{2^{\left(\alpha_{\rm
re}+\frac{1}{2}\right)}}\frac{1}{\Gamma\left(\alpha_{\rm
re}+\frac{3}{2}\right)}-B_{k,\,{\rm re}}\,\frac{2^{\left(\alpha_{\rm
re}+\frac{1}{2}\right)}}{\left(\alpha_{\rm re}\,y\right)^{2\alpha_{\rm
re}}}\frac{1}{\Gamma\left(-\alpha_{\rm re}+\frac{1}{2}\right)}\right]\,.$
(3.24)

Incorporating this into Eq. (3.2), we obtain

$A_{k,\,{\rm re}}=2^{\left(\alpha_{\rm
re}-\frac{1}{2}\right)}\,\Gamma\left(\alpha_{\rm re}+1/2\right)\,h_{k,\,{\rm
inf}}^{\lambda},\quad B_{k,\,{\rm re}}=0\,.$ (3.25)

  * •

Radiation-dominated epoch

The tensor mode functions in the RD epoch ($w=1/3,\,\alpha=1$) take the form

$h_{k,\,{\rm RD}}^{\lambda}(y)=\frac{1}{y^{\frac{1}{2}}}\left[A_{k,\,{\rm
RD}}\,J_{\frac{1}{2}}(y)+B_{k,\,{\rm RD}}\,J_{-\frac{1}{2}}(y)\right]\,.$
(3.26)

The coefficients for this epoch are determined using Eqs. (3.20) and (3.21) at
the end of reheating (denoted by ‘${\rm r*}$’)

$\displaystyle A_{k,\,{\rm RD}}$ $\displaystyle=\frac{y_{\rm
r*}^{(1-\alpha_{\rm re})}}{\alpha_{\rm re}^{\left(\alpha_{\rm
re}-\frac{1}{2}\right)}}\times A_{k,\,{\rm
re}}\left[\frac{J_{\left(\alpha_{\rm re}-\frac{1}{2}\right)}(\alpha_{\rm
re}\,y_{\rm r*})\,J_{-\frac{3}{2}}(y_{\rm r*})+J_{\left(\alpha_{\rm
re}+\frac{1}{2}\right)}(\alpha_{\rm re}\,y_{\rm r*})\,J_{-\frac{1}{2}}(y_{\rm
r*})}{J_{-\frac{1}{2}}(y_{\rm r*})\,J_{\frac{3}{2}}(y_{\rm
r*})+J_{\frac{1}{2}}(y_{\rm r*})\,J_{-\frac{3}{2}}(y_{\rm r*})}\right]\,,$
(3.27) $\displaystyle B_{k,\,{\rm RD}}$ $\displaystyle=\frac{y_{\rm
r*}^{(1-\alpha_{\rm re})}}{\alpha_{\rm re}^{\left(\alpha_{\rm
re}-\frac{1}{2}\right)}}\times A_{k,{\rm
re}}\left[\;\frac{J_{\left(\alpha_{\rm re}-\frac{1}{2}\right)}(\alpha_{\rm
re}\,y_{\rm r*})\,J_{\frac{3}{2}}(y_{\rm r*})\;-\;J_{\left(\alpha_{\rm
re}+\frac{1}{2}\right)}(\alpha_{\rm re}\,y_{\rm r*})\,J_{\frac{1}{2}}(y_{\rm
r*})}{J_{-\frac{1}{2}}(y_{\rm r*})\,J_{\frac{3}{2}}(y_{\rm
r*})+J_{\frac{1}{2}}(y_{\rm r*})\,J_{-\frac{3}{2}}(y_{\rm r*})}\;\right]\,,$
(3.28)

where $y_{\rm r*}=k/k_{\rm r*}$ with $k_{\rm r*}$ being the comoving
wavenumber that made its Hubble-entry at the end of reheating (beginning of
radiation domination).

  * •

Matter-dominated epoch

The tensor mode functions in the MD epoch ($w=0,\,\alpha=2$) take the form

$h_{k,\,{\rm MD}}^{\lambda}(y)=\frac{1}{(2\,y)^{\frac{3}{2}}}\left[A_{k,\,{\rm
MD}}\,J_{\frac{3}{2}}(2\,y)+B_{k,\,{\rm
MD}}\,J_{-\frac{3}{2}}(2\,y)\right]\,.$ (3.29)

Again, Eqs. (3.20) and (3.21) at matter-radiation equality (denoted by ‘${\rm
eq}$’) lead to

$\displaystyle A_{k,\,{\rm MD}}$ $\displaystyle=2^{\frac{3}{2}}\,y_{\rm
eq}\times\Bigg{\\{}A_{k,\,{\rm RD}}\left[\frac{J_{\frac{1}{2}}(y_{\rm
eq})\,J_{-\frac{5}{2}}(2\,y_{\rm eq})+J_{\frac{3}{2}}(y_{\rm
eq})\,J_{-\frac{3}{2}}(2\,y_{\rm eq})}{J_{-\frac{3}{2}}(2\,y_{\rm
eq})\,J_{\frac{5}{2}}(2\,y_{\rm eq})+J_{\frac{3}{2}}(2\,y_{\rm
eq})\,J_{-\frac{5}{2}}(2\,y_{\rm eq})}\right]$ $\displaystyle\hskip
85.35826pt+B_{k,\,{\rm RD}}\left[\frac{J_{-\frac{1}{2}}(y_{\rm
eq})\,J_{-\frac{5}{2}}(2\,y_{\rm eq})-J_{-\frac{3}{2}}(y_{\rm
eq})\,J_{-\frac{3}{2}}(2\,y_{\rm eq})}{J_{-\frac{3}{2}}(2\,y_{\rm
eq})\,J_{\frac{5}{2}}(2\,y_{\rm eq})+J_{\frac{3}{2}}(2\,y_{\rm
eq})\,J_{-\frac{5}{2}}(2\,y_{\rm eq})}\right]\Bigg{\\}}\,,$ (3.30)
$\displaystyle B_{k,\,{\rm MD}}$ $\displaystyle=2^{\frac{3}{2}}\,y_{\rm
eq}\times\Bigg{\\{}A_{k,\,{\rm RD}}\left[\frac{J_{\frac{1}{2}}(y_{\rm
eq})\,J_{\frac{5}{2}}(2\,y_{\rm eq})-J_{\frac{3}{2}}(y_{\rm
eq})\,J_{\frac{3}{2}}(2\,y_{\rm eq})}{J_{-\frac{3}{2}}(2\,y_{\rm
eq})\,J_{\frac{5}{2}}(2\,y_{\rm eq})+J_{\frac{3}{2}}(2\,y_{\rm
eq})\,J_{-\frac{5}{2}}(2\,y_{\rm eq})}\right]$ $\displaystyle\hskip
85.35826pt+B_{k,\,{\rm RD}}\left[\frac{J_{-\frac{1}{2}}(y_{\rm
eq})\,J_{\frac{5}{2}}(2\,y_{\rm eq})+J_{-\frac{3}{2}}(y_{\rm
eq})\,J_{\frac{3}{2}}(2\,y_{\rm eq})}{J_{-\frac{3}{2}}(2\,y_{\rm
eq})\,J_{\frac{5}{2}}(2\,y_{\rm eq})+J_{\frac{3}{2}}(2\,y_{\rm
eq})\,J_{-\frac{5}{2}}(2\,y_{\rm eq})}\right]\Bigg{\\}}\,,$ (3.31)

where $y_{\rm eq}=k/k_{\rm eq}$ with $k_{\rm eq}$ being the comoving
wavenumber that made its Hubble-entry at the time of matter-radiation
equality.

####  3.1.2 Universe with multiple equations of state during reheating

As stressed before, since reheating is a complex phase involving non-linear
processes, it is natural to represent it by multiple equation of state
parameters, as illustrated in Fig. 2. We model reheating by multiple epochs of
piece-wise constant EoS parameters with sharp (instantaneous) transitions in
between. The procedure to obtain the coefficients appearing in the solution
for the tensor mode functions is therefore the same as that given in Sec.
3.1.1, with some additional (minor) changes that one need to take into account
when considering multiple equations of state parameters, as detailed below.

Figure 2: A schematic depiction of the timeline of the universe with a post-
inflationary reheating phase consisting of multiple equations of state. Note,
the horizontal length is not a representation of the actual duration.

The coefficients for the first epoch following the end of inflation are the
same as that given in Eq. (3.25), for reheating with single EoS. The main
difference comes while labelling them. For example, if the equation of state
of the first epoch is $w_{1}$, then the coefficients of the mode functions are

$A_{k,\,1}=2^{\left(\alpha_{1}-\frac{1}{2}\right)}\,\Gamma\left(\alpha_{1}+1/2\right)\,h_{k,\,{\rm
inf}}^{\lambda},\quad B_{k,\,1}=0\,,$ (3.32)

where $\alpha_{1}=2/(1+3\,w_{1})$. Suppose, there are $n$ epochs during
reheating, each with a constant EoS ($w_{i}$ being the EoS of the $i^{\rm th}$
epoch, with $i\in\mathbb{Z^{+}}$). Then, the coefficients $A_{k,\,i}$ and
$B_{k,\,i}$ corresponding to all subsequent epochs ($i>1$) during reheating
need to be computed using Eqs. (3.20) and (3.21). For the sake of
completeness, we provide closed-form expressions for the coefficients at the
onset of radiation domination following the final transition ($i=n$), which
are different from the ones corresponding to a single transition as given in
Eqs. (3.27) and (3.28), to be

$\displaystyle A_{k,\,{\rm RD}}$ $\displaystyle=\frac{y_{\rm
r*}^{(1-\alpha_{n})}}{\alpha_{n}^{(\alpha_{n}-\frac{1}{2})}}\Bigg{\\{}A_{k,\,i}\left[\frac{J_{\left(\alpha_{n}-\frac{1}{2}\right)}(\alpha_{n}\,y_{\rm
r*})\,J_{-\frac{3}{2}}(y_{\rm
r*})+J_{\left(\alpha_{n}+\frac{1}{2}\right)}(\alpha_{n}\,y_{\rm
r*})\,J_{-\frac{1}{2}}(y_{\rm r*})}{J_{-\frac{1}{2}}(y_{\rm
r*})\,J_{\frac{3}{2}}(y_{\rm r*})+J_{\frac{1}{2}}(y_{\rm
r*})\,J_{-\frac{3}{2}}(y_{\rm r*})}\right]$ $\displaystyle\hskip
56.9055pt+B_{k,\,i}\left[\frac{J_{-\left(\alpha_{n}-\frac{1}{2}\right)}(\alpha_{n}\,y_{\rm
r*})\,J_{-\frac{3}{2}}(y_{\rm
r*})-J_{-\left(\alpha_{n}+\frac{1}{2}\right)}(\alpha_{n}\,y_{\rm
r*})\,J_{-\frac{1}{2}}(y_{\rm r*})}{J_{-\frac{1}{2}}(y_{\rm
r*})\,J_{\frac{3}{2}}(y_{\rm r*})+J_{\frac{1}{2}}(y_{\rm
r*})\,J_{-\frac{3}{2}}(y_{\rm r*})}\right]\Bigg{\\}}\,,$ (3.33) $\displaystyle
B_{k,\,{\rm RD}}$ $\displaystyle=\frac{y_{\rm
r*}^{(1-\alpha_{n})}}{\alpha_{n}^{(\alpha_{n}-\frac{1}{2})}}\Bigg{\\{}A_{k,n}\left[\;\frac{J_{\left(\alpha_{n}-\frac{1}{2}\right)}(\alpha_{n}\,y_{\rm
r*})\,J_{\frac{3}{2}}(y_{\rm
r*})\;-\;J_{\left(\alpha_{n}+\frac{1}{2}\right)}(\alpha_{n}\,y_{\rm
r*})\,J_{\frac{1}{2}}(y_{\rm r*})}{J_{-\frac{1}{2}}(y_{\rm
r*})\,J_{\frac{3}{2}}(y_{\rm r*})+J_{\frac{1}{2}}(y_{\rm
r*})\,J_{-\frac{3}{2}}(y_{\rm r*})}\;\right]$ $\displaystyle\hskip
56.9055pt+B_{k,\,n}\left[\;\frac{J_{-\left(\alpha_{n}-\frac{1}{2}\right)}(\alpha_{n}\,y_{\rm
r*})\,J_{\frac{3}{2}}(y_{\rm
r*})\;+\;J_{-\left(\alpha_{n}+\frac{1}{2}\right)}(\alpha_{n}\,y_{\rm
r*})\,J_{\frac{1}{2}}(y_{\rm r*})}{J_{-\frac{1}{2}}(y_{\rm
r*})\,J_{\frac{3}{2}}(y_{\rm r*})+J_{\frac{1}{2}}(y_{\rm
r*})\,J_{-\frac{3}{2}}(y_{\rm r*})}\;\right]\Bigg{\\}}\,,$ (3.34)

where $\alpha_{n}=2/(1+3\,w_{n})$. While, the analytical expression for
$A_{k,\,{\rm MD}}$ and $B_{k,\,{\rm MD}}$ remains the same as in Eqs. (3.30)
and (3.31).

###  3.2 Spectral energy density of gravitational waves

A cosmological background of stochastic gravitational waves at a given time is
characterized in terms of the spectral energy density of GWs, defined as the
energy density of GWs per unit logarithmic wavenumber interval relative to the
critical energy density of the universe,
$\rho_{c}(\tau)=3\,m_{p}^{2}\,H^{2}(\tau)$ at that epoch [145, 146], i.e.,

$\Omega_{{}_{\text{GW}}}(\tau,k)=\frac{1}{\rho_{c}}\,\frac{\mathrm{d}\rho_{{}_{\rm
GW}}}{\mathrm{d}\ln k}\,.$ (3.35)

The energy density of GWs obtained from the time-time component of the energy-
momentum tensor $\tensor{T}{{}^{\mu}_{\nu}}$ of tensor fluctuations in the
FLRW background (see App. C) is

$\hat{\rho}_{{}_{\rm
GW}}(\tau,\vec{x})=-\tensor{\hat{T}}{{}^{0}_{0}}(\tau,\vec{x})=\frac{m_{p}^{2}}{8a^{2}(\tau)}\left[(\hat{h}^{\prime}_{ij}(\tau,\vec{x}))^{2}+(\gradient\hat{h}_{ij}(\tau,\vec{x}))^{2}\right]\,.$
(3.36)

The vacuum expectation value of $\rho_{{}_{\rm GW}}$ (with summation over the
two polarisation states) is given by

$\rho_{{}_{\rm GW}}(\tau)=\left<0\left|\hat{\rho}_{{}_{\rm
GW}}(\tau,\vec{x})\right|0\right>=\frac{m_{p}^{2}}{8a^{2}(\tau)}\int\mathrm{d}\ln
k\;\frac{k^{3}}{\pi^{2}}\Big{[}\left|h_{k}^{{}^{\prime}\lambda}(\tau)\right|^{2}+k^{2}\left|h_{k}^{\lambda}(\tau)\right|^{2}\Big{]}\,,$
(3.37)

where we used the following Fourier representation of the tensor operator
$\hat{h}_{ij}$

$\hat{h}_{ij}(\tau,\vec{x})=\sum_{\lambda=+,\times}\int\frac{\mathrm{d}^{3}\vec{k}}{(2\pi)^{3}}\,\epsilon^{\lambda}_{ij}\,\left[h_{k}^{\lambda}(\tau)\,\hat{a}_{\vec{k}}\,e^{i\vec{k}\cdot\vec{x}}+\left(h_{k}^{\lambda}(\tau)\right)^{*}\,\hat{a}_{\vec{k}}^{{\dagger}}\,e^{-i\vec{k}\cdot\vec{x}}\right]\,,$
(3.38)

and the vacuum state is defined by $\hat{a}_{\vec{k}}\left|0\right>=0$. Using
Eq. (3.37), expression for the GW spectral energy density becomes

$\Omega_{{}_{\text{GW}}}(\tau,k)=\frac{k^{2}}{24\,a^{2}(\tau)H^{2}(\tau)}\;{\cal
P}_{T}(\tau,k)\,\left[1+\frac{1}{k^{2}}\,\frac{\left|h_{k}^{{}^{\prime}\lambda}(\tau)\right|^{2}}{\left|h_{k}^{\lambda}(\tau)\right|^{2}}\right]\,.$
(3.39)

Considering modes that are deep in the Hubble radius so that the frequency of
GWs is substantially larger than the expansion rate of the universe, i.e.,
$k\gg aH$, the second term inside the square bracket becomes [147, 148]

$\frac{1}{k^{2}}\,\frac{\left|h_{k}^{{}^{\prime}\lambda}(\tau)\right|^{2}}{\left|h_{k}^{\lambda}(\tau)\right|^{2}}=1+{\cal
O}\left(\frac{a^{2}H^{2}}{k^{2}}\right)\simeq 1\,.$ (3.40)

Therefore, in the sub-Hubble regime, expression for the GW spectral energy
density becomes

$\Omega_{{}_{\text{GW}}}(\tau,k)=\frac{k^{2}}{12\,a^{2}(\tau)H^{2}(\tau)}\;{\cal
P}_{T}(\tau,k)\,,$ (3.41)

which coincides with the standard Isaacson approximation [149, 150]. We again
emphasize that this approximation is valid as long as the wavelength of GWs
are much shorter than the Hubble radius of the universe333In the cosmological
context, there are long wavelength tensor modes comparable to the Hubble scale
that might be of observational relevance. Hence, in general, one needs to
derive the energy-momentum tensor without invoking the frequently used
approximation of high frequency GWs propagating in flat background [146, 143].
Nevertheless, in this work we will continue using Eq. (3.41) since we are
primarily interested in high frequency GWs that are deep inside the comoving
Hubble radius in the MD epoch..

In order to calculate the spectral energy density at the present epoch, we
compute tensor mode functions deep inside the Hubble radius, in the limit
$y\gg 1$, during the MD epoch by making use of the large argument expansion of
Bessel function, namely, $J_{\nu}(y)\approx\sqrt{2/(\pi
y)}\,\cos(y-\nu\frac{\pi}{2}-\frac{\pi}{4})$. Accordingly, Eq. (3.29) reduces
to

$h_{k,\,{\rm MD}}^{\lambda}(y\gg
1)=-\frac{1}{(2\,y)^{\frac{3}{2}}}\,\frac{1}{\sqrt{\pi y}}\left[A_{k,\,{\rm
MD}}\,\cos(2\,y)+B_{k,\,{\rm MD}}\,\sin(2\,y)\right]\,.$ (3.42)

Consequently, the average of the modulus squared value of tensor amplitude
becomes

$\overline{|h_{k,\,{\rm MD}}^{\lambda}(y\gg
1)|^{2}}=\frac{1}{(2\,y)^{4}}\,\frac{1}{\pi}\left[|A_{k,\,{\rm
MD}}|^{2}+|B_{k,\,{\rm MD}}|^{2}\right]\,.$ (3.43)

Upon substituting $y(\tau)=y_{\rm eq}\,\frac{a_{\rm eq}H_{\rm
eq}}{a(\tau)H(\tau)}=y_{\rm eq}\,\left(\frac{a(\tau)}{a_{\rm
eq}}\right)^{1/2}$, for $\tau>\tau_{\rm eq}$, into the above expression, we
can observe that the term $\overline{|h_{k,\,{\rm MD}}^{\lambda}(y\gg
1)|^{2}}\propto a^{-2}$. In fact, this is true in general for any epoch, i.e.,
$\overline{|h_{k}^{\lambda}(y\gg 1)|^{2}}\propto a^{-2}$. This implies that
the tensor amplitude dampens as $h_{k}^{\lambda}\propto a^{-1}$ due to the
expansion of the universe, irrespective of the EoS [119]. At late times, as
the universe begins to accelerate again, no new additional tensor modes become
sub-Hubble in this epoch. Instead, a narrow range of long-wavelength tensor
modes, that had freshly become sub-Hubble, begin to exit the Hubble radius
during this accelerated epoch. We therefore consider these tensor modes
inaccessible, and only focus on sub-Hubble modes deep inside the Hubble radius
at the onset of the MD epoch. Substituting Eq. (3.43) and $k=y_{\rm
eq}\,a_{\rm eq}\,H_{\rm eq}$ in Eq. (3.41) we obtain

$\displaystyle\overline{\Omega_{{}_{\rm
GW}}(\tau_{0},k)}=\frac{1}{96\pi^{3}}\,\frac{g_{*,\,{\rm
r*}}}{g_{*,\,0}}\left(\frac{g_{s,\,0}}{g_{s,\,{\rm
r*}}}\right)^{4/3}\,\Omega_{\rm rad,\,0}\,\frac{1}{y_{\rm
eq}^{2}}\left[\tilde{A}^{2}_{k,\,{\rm MD}}+\tilde{B}^{2}_{k,\,{\rm
MD}}\right]\,\frac{H^{2}_{\rm
inf}}{m_{p}^{2}}\left(\frac{k}{k_{*}}\right)^{n_{{}_{T}}}\,,$ (3.44)

where $g_{*},\,g_{s}$ are the effective number of relativistic degrees of
freedom of the universe in energy density, entropy density respectively, with
the second subscript in $g_{*}$ and $g_{s}$ denoting the time of estimation of
its value (with ‘0’ denoting the present epoch). $H_{{\rm inf}}$ corresponds
to the Hubble parameter at the Hubble-exit of the CMB pivot scale during
inflation, i.e., $H_{{\rm inf}}=H_{*}$. Note that in Eq. (3.44), the
coefficients $\tilde{A}_{k,\,{\rm MD}}$ and $\tilde{B}_{k,\,{\rm MD}}$ differ
from $A_{k,\,{\rm MD}}$ and $B_{k,\,{\rm MD}}$ in Eq. (3.29) by the factor
$h_{k,\,{\rm inf}}^{\lambda}$ as discussed in App. D.

From detection perspective, it is important to characterize the GW spectral
energy density in terms of the present-day (physical, as opposed to comoving)
frequency of GWs [77, 119, 64, 117, 147]. The present-day frequency $f_{k}$ of
a mode of comoving wavenumber $k$ is given by

$f_{k}=\frac{k}{2\pi a_{0}}=\frac{1}{2\pi}\frac{a_{k}}{a_{0}}H_{k}\,,$ (3.45)

where $a_{k}$ and $H_{k}$ are the scale factor and Hubble parameter at the
time when the mode with comoving wavenumber $k$ made its Hubble-entry.

In order for us to mark the beginning and the end of each post-inflationary
epoch in our computation, we relate the present-day frequency of a mode to the
temperature at the time of its Hubble-entry. In a radiation-dominated
universe, the Hubble parameter can be expressed in terms of its temperature as

$H^{2}(T)=\frac{\rho_{\rm
rad}(T)}{3\,m_{p}^{2}}=\frac{1}{3\,m_{p}^{2}}\,\frac{\pi^{2}}{30}\,g_{*,\,T}\,T^{4}\,.$
(3.46)

From the conservation of entropy, we have

$\left(\frac{g_{s,\,0}}{g_{s,T}}\right)^{1/3}\,\frac{T_{0}}{T}=\frac{a(T)}{a_{0}}\,.$
(3.47)

By incorporating Eqs. (3.46) and (3.47) to Eq. (3.45), one can relate the
temperature $T_{k}$ of the universe (in GeV) at the Hubble-entry of a mode to
its present-day frequency $f_{k}$ (in Hz) as [34]

$\displaystyle\frac{f_{k}}{{\rm Hz}}=7.43\times
10^{-8}\left(\frac{g_{s,\,0}}{g_{s,T_{k}}}\right)^{1/3}\,\left(\frac{g_{*,T_{k}}}{90}\right)^{1/2}\,\left(\frac{T_{k}}{{\rm
GeV}}\right)\,.$ (3.48)

Note that the above expression is valid for a radiation-dominated universe in
thermal equilibrium. With reheating being primarily a non-thermal process,
achieving thermal equilibrium only at the end of it, Eq. (3.48) cannot be used
during reheating. However, we shall be using Eq. (3.48) to relate the energy
scale of the universe during reheating to the physical frequency of GWs; where
$T_{k}$ is treated as an effective temperature (related to $f_{k}$ via Eq.
(3.48), which serves as a proxy for the energy-scale at the Hubble-entry time
of the mode. In fact, to be precise, we can define the corresponding energy
scale of the universe to be

$E_{k}=\rho_{c}^{1/4}(T_{k})=\left(\frac{\pi^{2}}{30}\,g_{*,\,T_{k}}\right)^{1/4}\,T_{k}\,.$
(3.49)

Therefore, the present-day frequency of a mode making its Hubble-entry at an
energy scale $E_{k}$ is given by

$\frac{f_{k}}{{\rm Hz}}=1.03\times
10^{-8}\left(\frac{g_{s,\,0}}{g_{s,T_{k}}}\right)^{1/3}\,g_{*,T_{k}}^{1/4}\,\left(\frac{E_{k}}{{\rm
GeV}}\right)\,.$ (3.50)

Table 1 lists the present-day frequency of GWs generated at the Hubble-entry
of some of the important cosmological events (and their corresponding energy
scales).

While computing the spectral energy density of GWs, in order to be concrete we
have fixed the energy scale during inflation to be $E_{\rm inf}=5.76\times
10^{15}$ GeV corresponding to the tensor-to-scalar ratio $r=0.001$, which is a
target of next generation CMB experiments. Similarly, in order to be
consistent with the BBN predictions for the abundance of light nuclei,
reheating must end before the onset of BBN, i.e., $T_{\rm r*}\gg T_{\rm
BBN}\sim 1$ MeV. Therefore, we fix the energy scale of the universe at the end
of reheating to be $E_{\rm r*}=1\,{\rm GeV}\;(T_{\rm r*}\simeq 0.45\,{\rm
GeV})$, which is safely above the BBN abound. Note that for a given reheating
history, by decreasing the tensor-to-scalar ratio, the probability of
detection of GWs will become lower. Therefore, our analysis in this work must
be treated as a best case scenario for the detection. By fixing the reheating
temperature, we primarily focus on illustrating how to probe the duration and
EoS of successive epochs during reheating.

Event | Energy scale ($E_{k}$) | Present-day frequency $f_{k}$ (Hz)
---|---|---
Matter-radiation equality |  $\sim 1$ eV | $1.4\times 10^{-17}$
Hubble-entry time of CMB pivot scale |  $\sim 5$ eV | $7.2\times 10^{-17}$
Onset of BBN |  $\sim 1.4$ MeV | $1.8\times 10^{-11}$
QCD Phase Transition |  $\sim 320$ MeV | $3.7\times 10^{-9}$
Electro-weak symmetry breaking |  $\sim 240$ GeV | $2.7\times 10^{-6}$
Table 1:  Present-day frequencies of primordial GWs generated at different
(Hubble-entry) energy scale corresponding to some of the important events in
the thermal history of our universe.

By replacing the comoving wavenumber $k$ with the present-day frequency
$f_{k}$, and denoting the oscillation-averaged spectral energy density as
$\Omega_{{}_{\rm GW}}^{(0)}(f)$ instead of $\overline{\Omega_{{}_{\rm
GW}}^{(0)}(f)}$, we obtain

$\displaystyle h^{2}\,\Omega_{{}_{\rm
GW}}^{(0)}(f_{k})=\frac{1}{96\pi^{3}}\,\frac{g_{*,\,{\rm
r*}}}{g_{*,\,0}}\left(\frac{g_{s,\,0}}{g_{s,\,{\rm
r*}}}\right)^{4/3}\,h^{2}\Omega_{\rm rad,\,0}\,\left(\frac{f_{\rm
eq}}{f_{k}}\right)^{2}\left[\tilde{A}^{2}_{k,\,{\rm
MD}}+\tilde{B}^{2}_{k,\,{\rm MD}}\right]\,\frac{H^{2}_{\rm
inf}}{m_{p}^{2}}\,,$ (3.51)

where $h$ is a dimensionless parameter, related to the Hubble parameter at the
present epoch by $H_{0}=100\,h\,{\rm km}\,{\rm s}^{-1}\,{\rm Mpc}^{-1}$.
Values of the relativistic degrees of freedom at our fixed reheating
temperature are $g_{*,\,{\rm r*}}=g_{s,\,{\rm
r*}}=73.0,\,g_{s,\,0}=3.94,\,{\rm and}\,g_{*,\,0}=3.38$. The present-day
density parameter of radiation is $h^{2}\Omega_{\rm rad,\,0}=4.16\times
10^{-5}$ (see Eq. (E.7)). The spectral-tilt of relic GWs is defined as [58]

$n_{{}_{\rm GW}}=\frac{\mathrm{d}\ln\Omega_{{}_{\rm
GW}}^{(0)}(f_{k})}{\mathrm{d}\ln f_{k}}\,,$ (3.52)

which, for a range of frequencies corresponding to modes making their Hubble-
entry at the beginning and end of a post-inflationary cosmic epoch with
constant EoS $w$, reduces to444 Note that we have ignored contribution from
the tiny red-tilt $n_{{}_{T}}\simeq-r/8$ of the inflationary tensor modes,
since they are very small for $r\leq 0.001$, see Ref. [14]. [34]

$n_{{}_{\rm GW}}=2\left(\frac{w-1/3}{w+1/3}\right)\,.$ (3.53)

Figure 3: The spectral energy density of first-order inflationary GWs is shown
as a function of the present-day GW frequencies for a single constant post-
inflationary EoS during reheating; along with the sensitivity curves of
various GW detectors [151], as well as the BBN constraint. The spectrum is
shown for three different possibilities of reheating EoS above $f\gtrsim
10^{-8}\,{\rm Hz}$, corresponding to a blue-tilted stiff matter phase with
$w_{{}_{\rm re}}=0.8$ (solid black), radiative phase with $w_{{}_{\rm
re}}=1/3$ (dotted black), and softer matter-like phase with $w_{{}_{\rm
re}}=0$ (dashed black). The tensor-to-scalar ratio has been set to $r=0.001$
and the reheating energy scale is fixed at $E_{\rm r*}=1$ GeV.

Fig. 3 shows the GW spectral energy density as a function of the present-day
frequency of GWs for the case of a single constant EoS during reheating.
During the MD epoch, the energy density of relic GWs redshifts faster than the
background density, implying that $\Omega_{{}_{\rm GW}}(f)$ exhibits a red-
tilt, with slope $n_{{}_{\rm GW}}\simeq-2$ as determined from Eq. (3.53).
Similarly, $\Omega_{{}_{\rm GW}}$ is (nearly) flat for the range of
frequencies corresponding to modes that made their Hubble-entry in the RD
epoch, since the GWs redshift at an equal rate in comparison to the background
radiation. Finally, the spectrum of those modes that made their Hubble-entry
during reheating can be either blue-tilted, or flat or red-tilted depending
upon the EoS during reheating. For example, we see a red-tilt for a softer EoS
$w_{{}_{\rm re}}<1/3$ (black dashed curve), and a blue-tilt for a stiff-matter
phase $w_{{}_{\rm re}}>1/3$ (solid black curve); while for radiation like EoS
$w_{{}_{\rm re}}=1/3$ (black dotted curve), the spectrum is flat, as expected.
Therefore, we observe that if the post-inflationary reheating EoS is stiff-
matter like, then the relic GWs have a greater likelihood of getting detected
by the upcoming GW detectors [152, 119, 34, 153].

Keeping in mind that GW energy density contributes to the total energy density
of the universe, the observations on the light nuclei abundance impose an
upper bound on the GW energy density during BBN to be [61]

$\rho_{{}_{\rm GW}}(T_{\rm BBN})\leq\Delta\rho_{{}_{\rm rad}}(T_{\rm
BBN})=\frac{\pi^{2}}{30}\times\frac{7}{8}\,(2\,\Delta N_{\nu})\,T_{\rm
BBN}^{4}=\frac{7}{8}\,\rho_{\gamma}(T_{\rm BBN})\,\Delta N_{\nu}\,,$ (3.54)

where $\Delta N_{\nu}$ is the extra neutrino (relativistic) degrees of freedom
outside the Standard Model from the Standard Model value $N_{\nu}=3$, and
$\rho_{\gamma}$ is the energy density of photons. Using entropy conversation
from Eq. (3.47), we obtain the expression for the total energy density of
primordial GWs that are sub-Hubble before the BBN to be

$\frac{\rho_{{}_{\rm GW}}(T_{\rm BBN})}{\rho_{{}_{\rm
GW},\,0}}=\left(\frac{a_{0}}{a_{\rm BBN}}\right)^{4}=\left(\frac{g_{s,\,T_{\rm
BBN}}}{g_{s,\,0}}\right)^{4/3}\left(\frac{T_{\rm BBN}}{T_{0}}\right)^{4}\,.$
(3.55)

Figure 4: The spectral energy density ($\Omega_{{}_{\rm GW}}$) of inflationary
GWs associated with two (top row), three (middle row) and four (bottom row)
epochs of piece-wise constant EoS parameters during reheating. Plots in the
left column correspond to scenarios where $\Omega_{{}_{\rm GW}}$ has a minimal
overlap with the LISA sensitivity curve, while those in the right column
belong to scenarios where $\Omega_{{}_{\rm GW}}$ has a substantial overlap
with LISA sensitivity region, without violating the existing BBN and aLIGO
constraints. The energy scales $E_{i}$ (with $i\in\mathbb{Z}^{+}$)
corresponding to the end of each epoch, are given in the top labels of each
plot. The energy scales at the end of reheating and inflation are the same in
all plots, namely, $E_{\rm r*}=1$ GeV and $E_{\rm inf}=5.76\times 10^{15}$ GeV
respectively, so as the tensor-to-scalar ratio, $r=0.001$.

We can further modify Eq. (3.55) using the relation involving energy density
of photons, $\rho_{\gamma}(T)=(\pi^{2}/15)\,T^{4}$,

$\frac{\rho_{{}_{\rm GW},\,0}}{\rho_{c,\,0}}=\frac{\rho_{{}_{\rm GW}}(T_{\rm
BBN})}{\rho_{\gamma}(T_{\rm BBN})}\,\left(\frac{g_{s,\,0}}{g_{s,\,T_{\rm
BBN}}}\right)^{4/3}\,\Omega_{\gamma,\,0}\,.$ (3.56)

Using Eqs. (3.35) and (3.54), we can rewrite Eq. (3.56) as

$\displaystyle h^{2}\int_{f_{\rm BBN}}^{f_{e}}\mathrm{d}\ln f\;\Omega_{{}_{\rm
GW}}(\tau_{0},\,f)$ $\displaystyle\leq\frac{7}{8}\,\Delta
N_{\nu}\,\left(\frac{g_{s,\,0}}{g_{s,\,T_{\rm
BBN}}}\right)^{4/3}\,h^{2}\,\Omega_{\gamma,\,0}$ $\displaystyle=5.64\times
10^{-6}\,\Delta N_{\nu}<1.13\times 10^{-6}\,,$ (3.57)

where we have used the values $h^{2}\,\Omega_{\gamma,\,0}=2.46\times 10^{-5}$
(see Eq. (E.3))$,\,g_{s,\,0}=3.94,\,g_{s,\,T_{\rm BBN}}=10.75$. The upper
bound, $N_{\nu}<3.2$, is obtained with a $95\%$ C.L. from combining the
deuterium-to-hydrogen ratio along with CMB baryon density [154]. Eq. (3.2) is
integrated from the frequencies of those modes that have entered the comoving
Hubble radius at the onset of BBN $(f_{\rm BBN})$ to those at the end of
inflation $(f_{e})$. For single-epoch reheating dominated by stiff-matter,
with EoS $w_{\rm re}$, the aforementioned integral is dominated by the upper
bound on $\Omega_{{}_{\rm GW}}$ at $f_{e}$ corresponding to the end of
inflation. Hence, Ref. [119] used the following weaker constraint instead of
Eq. (3.2) for numerical computation

$h^{2}\Omega_{{}_{\rm GW}}(\tau_{0},\,f_{e})<(1-\alpha_{\rm re})\,2.26\times
10^{-6}\,,$ (3.58)

where $\alpha_{\rm re}=2/(1+3w_{\rm re})$. In the case of reheating with
multiple epochs, we also use a similar constraint

$h^{2}\Omega_{{}_{\rm
GW}}(\tau_{0},\,f_{e})<\begin{cases}(1-\alpha_{1})\,2.26\times
10^{-6},&w_{1}>\frac{1}{3}\,,\\\ 2.24\times
10^{-8},&w_{1}<\frac{1}{3}\,,\end{cases}$ (3.59)

where $\alpha_{1}=2/(1+3w_{1})$ and $w_{1}$ is the EoS parameter of the first
epoch during reheating. Throughout this paper, we use Eq. (3.59) as the BBN
constraint instead of Eq. (3.2) for numerical computation as the latter is
computationally expensive. However, we emphasise that this particular choice
leading to a weaker constraint may not be reliable for other values of $r$ and
$E_{\rm r*}$, which results in an increased value of the maximum allowed area
under the GW spectrum from $f_{\rm BBN}$ to $f_{e}$.

Fig. 4, displays the GW spectral energy density for the cases of reheating
with two, three, and four epochs of piece-wise constant EoS parameters; with
each epoch having a different EoS and duration. The figure helps to highlight
post-inflationary scenarios where the predicted signals for the spectral
energy density of GWs have, approximately, either a minimal overlap (left
column) or maximal overlap (right column) with the LISA sensitivity curve,
without violating the existing aLIGO and BBN constraints. Hence, it is
important to analyse the possible combinations of EoS parameters and duration
of these epochs resulting in signals that are detectable in the near future,
not only by LISA, but also by other upcoming GW detectors.

The case of instantaneously sharp transitions, apart from being
computationally simple, gives a reasonable insight into the detectable post-
inflationary properties of the universe. However, we must stress that
realistic sharp cosmological transitions are not expected to be instantaneous,
but rather smooth [59]. Hence, our $\Omega_{{}_{\rm GW}}$ curves are not fully
accurate for frequencies right around the transitions.

##  4 Inflationary gravitational waves as a probe of the primordial equation
of state

In this section, we determine the subspace of equation of state parameters
(and their corresponding duration) during reheating, comprising of two and
three epochs of piece-wise constant EoS, which leads to a detectable signal
for the inflationary GWs in the upcoming GW detectors555Important additional
constraints come from the bound on CMB spectral distortions [155, 156] which
are more relevant for relatively lower frequency GW signals in between the CMB
and PTA sensitivity curves. such as LISA, CE, ET, BBO, and DECIGO. As
mentioned earlier, all of the following results are obtained assuming a
tensor-to-scalar ratio $r=0.001$, and for a fixed reheating energy scale
$E_{\rm r*}=1$ GeV. Note that even though the post-inflationary EoS can vary
between $-1/3<w\leq 1$, we restrict ourselves to $-0.28\leq w\leq 1$ due to
computational limitations associated with the values of EoS closer to
$w=-1/3$, which lead to divergences in our numerical computation. The
aforementioned GW detectors can be classified into three groups based on their
optimal frequency range of operation (see Fig. 5): $10^{-4}-10^{-1}\,{\rm Hz}$
for LISA, $10^{-3}-10\,{\rm Hz}$ for BBO and DECIGO, and $1-10^{3}\,{\rm Hz}$
for aLIGO, CE and ET [157, 158, 88, 159, 160, 161].

Figure 5: GW spectral energy density for single-epoch reheating with
$w_{{}_{\rm re}}=0.73$. In this scenario, $\Omega_{{}_{\rm GW}}$ has a (near)
minimal overlap with the LISA sensitivity curve. The spectrum is plotted for
$r=0.001$ and $E_{\rm r*}=1$ GeV.

The EoS parameter space, $\\{w_{1}\,,\leavevmode\nobreak\ w_{2}\\}$, for
reheating with two epochs is shown in Fig. 6, with $E_{1}$ being the energy
scale of the universe at the transition between the two epochs. In order for
the signal to possess a greater likelihood of detection with the upcoming GW
detectors (without violating BBN and aLIGO bounds), a higher value of $w_{2}$
is preferred, as evident from Fig. 6. However, having higher values of both
$w_{1}$ and $w_{2}$ would violate the BBN constraint given in Eq. (3.59) as
indicated by the grey shaded region. Similarly, the dark maroon-shaded region
represents the subspace outside the BBN constraint that are ruled out by the
null detection of stochastic GWs by aLIGO. From Fig. 6, we also notice that
with an increase in the value of $E_{1}$, the range of values of $w_{2}$ that
leads to a detection in the future detectors actually decreases; while that of
$w_{1}$ becomes relatively unconstrained. This is because the physical
frequency of GWs corresponding to the Hubble-entry at the end of the first
epoch during reheating becomes higher than the frequency ranges of the
aforementioned GW detectors. As a result, these GW detectors fails to put any
constraint on $w_{1}$; as evident from Figs. 6, 6 and 6. Similar things happen
for each detector, albeit at different values of $E_{1}$. Specifically, the
detectable values of $w_{2}$ for LISA reaches its minimum around $E_{1}\geq
10^{6}$ GeV. For BBO and DECIGO, this happens around $E_{1}\geq 10^{8}$ GeV,
while for CE and ET, operating at higher frequencies, the same takes place
around $E_{1}\geq 10^{10}$ GeV.

Furthermore, we can also constrain the duration of each epoch, in terms of the
number of $e$-folds $N$, corresponding to an EoS ($w$), given the energy
scales at the beginning ($E_{i}$) and end ($E_{f}$) of that epoch. Since the
energy density scales as $\rho\simeq E^{4}\propto a^{-3(1+w)}$, the relation
between $N$ and $w$ is given by

$N=\ln(\frac{a_{f}}{a_{i}})=\frac{4}{3(1+w)}\,\ln(\frac{E_{i}}{E_{f}})\,.$
(4.1)

Accordingly, we have also provided the corresponding number of $e$-folds in
the plots appearing in Fig. 6 for each value of $w$.

Following the above procedure, one can obtain the parameter space for a three-
epoch reheating phase with the EoS parameters $\\{w_{1},\,w_{2},\,w_{3}\\}$.
In order to clearly describe the constraints on such a 3-dimensional parameter
space, we have considered illustrating the results by using the 2-dimensional
projections of the full plots. The projected two-dimensional plots for
different values of $E_{1}$ and $E_{2}$ (energy scales of the universe at the
end of the first and the second epochs, respectively) are shown in Figs. 7, 8
and 9.

A careful reader must take caution in interpreting these plots since they are
projections. In fact, the subspace of $\\{w_{1},\,w_{2},\,w_{3}\\}$ that leads
to a potentially detectable GW signal in the future detectors, is actually
smaller than that shown in these projected plots. Specifically, the
interpretation of whether a point in the projected space will lead to a
potential detection is quite murky towards the edges of the shaded regions.
Keeping this in mind, we shall briefly discuss a general prescription that the
reader may follow while interpreting these plots, where we use a tuple
$(w_{1},\,w_{2},\,w_{3})$ to represent a point in the EoS parameter space
$\\{w_{1},\,w_{2},\,w_{3}\\}$. The prescription may seem complicated at first,
but with practice, it becomes much easier to interpret the 2-dimensional plots
correctly.

Step-1:

Start by selecting a GW detector that is of interest to the reader, amongst
LISA, BBO, CE, ET and DECIGO.

Step-2:

Pick a value for $w_{1}$ (denoting it as $w_{1}^{*}$) within the shaded region
(of interest) in the left column of Figs. 7, or 8 or 9 for that detector. For
this choice ($w_{1}^{*}$), find the range of $w_{2}$ and $w_{3}$ from the left
and the middle columns, respectively, which leads to a detectable signal.

Step-3:

Succeeding this step, choose a value for $w_{2}$ (denoting it as $w_{2}^{*}$)
from the aforementioned range of $w_{2}$. For that particular choice
($w_{2}^{*}$), determine the range of detectable values of $w_{3}$ from the
right columns.

Step-4:

Next, obtain the intersection from aforementioned ranges of values of $w_{3}$,
as obtained from the second and the third steps. We denote this resultant
range as $W_{3}$, which may lead to a detectable GW signal. However, one has
to be careful as we are dealing with projected 2-dimensional constraints of
the fully three-dimensional parameter space. Thus, not all values in $W_{3}$
will guarantee a detection. To confirm whether any chosen
$(w_{1}^{*},\,w_{2}^{*},\,w_{3}^{*})$, with $w_{3}^{*}\in W_{3}$, results in a
detection, we have to check if $w_{3}^{*}>w_{3}^{\rm min}$, where $w_{3}^{\rm
min}\in W_{3}$ is the minimum value that leads to a detectable signal.

In general, it is difficult to find $w_{3}^{\rm min}$ for a given $w_{1}^{*}$
and $w_{2}^{*}$ from the projected plots. However, if the lower boundary of
the sensitivity curve of the detector near $w_{1}^{*}$ is (nearly) flat in the
left column plots, then $w_{3}^{\rm min}$ will approximately be equal to the
lower bound of $W_{3}$.

Step-5:

If the projection of a given point ($w_{1}^{*},\,w_{2}^{*},\,w_{3}^{*}$) is
outside BBN and aLIGO constraint curves in any of the three columns, then we
can conclude that the combination ($w_{1}^{*},\,w_{2}^{*},\,w_{3}^{*}$)
satisfies the BBN and aLIGO constraints. On the other hand, it is difficult,
in general, to conclusively infer whether a point violates the BBN or aLIGO
constraints or not.

Figure 6: The parameter space for two-epoch reheating phase with EoS
parameters $w_{1}$ and $w_{2}$ that leads to a potentially detectable GW
signal in the LISA, BBO, DECIGO, CE, ET, and aLIGO detectors. The grey-shaded
region in each plot corresponds to the combination of $w_{1}$ and $w_{2}$ that
violate the BBN constraint (as provided in Eq. (3.59)), while the maroon-
shaded region represents those that are ruled out by the aLIGO null detection.
Energy scale of the universe at the end of the first post-inflationary epoch
is labelled above in each plot. Note that, the duration of each epoch (number
of $e$-folds) has been provided opposite to the corresponding EoS axis.

Figure 7: The parameter space for three-epoch reheating phase with EoS
parameters $w_{1},\,w_{2}$ and $w_{3}$ that leads to a potentially detectable
GW signal in the LISA, BBO, DECIGO, CE, ET, and aLIGO detectors. The plots are
obtained by taking 2-dimensional projections of the full 3-dimensional
subspace of EoS parameters. Energy scales at which the first and second post-
inflationary epochs end are given as $E_{1}$ and $E_{2}$. The duration of each
epoch (number of $e$-folds) has been provided opposite to the corresponding
EoS axis. The projections of the BBN-constrained and aLIGO-constrained regions
are enclosed by the dashed and dotted black lines respectively.

Figure 8: The parameter space for three-epoch reheating phase with EoS
parameters $w_{1},\,w_{2}$ and $w_{3}$ that leads to a potentially detectable
GW signal in the LISA, BBO, DECIGO, CE, ET, and aLIGO detectors. The plots are
obtained by taking 2-dimensional projections of the full 3-dimensional
subspace of EoS parameters. Energy scales at which the first and second post-
inflationary epochs end are given as $E_{1}$ and $E_{2}$. The duration of each
epoch (number of $e$-folds) has been provided opposite to the corresponding
EoS axis. The projections of the BBN-constrained and aLIGO-constrained regions
are enclosed by the dashed and dotted black lines respectively.

Figure 9: The parameter space for three-epoch reheating phase with EoS
parameters $w_{1},\,w_{2}$ and $w_{3}$ that leads to a potentially detectable
GW signal in the LISA, BBO, DECIGO, CE, ET, and aLIGO detectors. The plots are
obtained by taking 2-dimensional projections of the full 3-dimensional
subspace of EoS parameters. Energy scales at which the first and second post-
inflationary epochs end are given as $E_{1}$ and $E_{2}$. The duration of each
epoch (number of $e$-folds) has been provided opposite to the corresponding
EoS axis. The projections of the BBN-constrained and aLIGO-constrained regions
are enclosed by the dashed and dotted black lines respectively.

In the following, we will demonstrate how to utilize the aforementioned
prescription with a number of key examples.

Example 1:

Consider the case with $E_{1}=10^{6}$ GeV and $E_{2}=10^{3}$ GeV which is
shown in Fig. 7. Suppose we are interested in determining whether a given
combination of the three EoS parameters results in a detectable signal in the
LISA detector (purple-shaded region) (following step-1). The range of $w_{1}$
that leads to a potential detection is $[-0.28,\,0.74]$. Fixing
$w_{1}^{*}=0.4$ in the left and middle columns, we obtain the range of $w_{2}$
which likely to give detection to be $[0.54,\,1]$, and that of $w_{3}$ to be
$[0.56,\,1]$ (following step-2). Note that this does not imply any combination
$(0.4,\,w_{2}^{*},\,w_{3}^{*})$ with $w_{2}^{*}$ and $w_{3}^{*}$ within the
aforementioned ranges will result in a potential detection, as stressed
before. Choosing $w_{2}^{*}=0.6\in[0.54,\,1]$, we observe that the allowed
range of $w_{3}$ to be $[0.84,\,1]$ from the right column in Fig. 7 (following
step-3), with this range being independent of $w_{1}^{*}$. After taking the
intersection of the two aforementioned ranges of $w_{3}$, we obtain
$W_{3}=[0.84,\,1]$. In the left column, we observe that the lower boundary
near $w_{1}^{*}=0.4$ is quite flat for LISA. Hence, the value $w_{3}^{\rm
min}$ will be approximately equal to the lower bound of $W_{3}$, i.e.,
$w_{3}^{\rm min}\simeq 0.84$ (following step-4). We can, therefore, be almost
certain that the combination, (0.4, 0.6, $w_{3}^{*}$) for
$w_{3}^{*}\in[0.84,\,1]$ will result in a GW signal that will be detectable by
LISA.

Finally, we carry out the last check regarding whether any point (0.4, 0.6,
$w_{3}^{*}$) for $w_{3}^{*}\in[0.84,\,1]$ violates the BBN and aLIGO
constraints. From the left column, we find that $w_{1}=0.4$ lies outside the
BBN-constrained region. Therefore, for $w_{1}=0.4$, all values of
$w_{2}\in[0.54,\,1]$ and $w_{3}\in[0.56,\,1]$ satisfy the BBN constraint
(following step-5). Similarly, we find that there is no region ruled out by
the aLIGO constraints in Fig. 7. Hence, the aforementioned collection of EoS
parameters leads to a potentially detectable GW signal in the LISA detector,
without violating the existing constraints.

Example 2:

Consider the case of $E_{1}=10^{8}$ GeV and $E_{2}=10^{5}$ GeV shown in Fig.
8. This time, let us look for a detectable signal in the CE detector (red-
shaded region) (following step-1). The detectable range for $w_{1}$ from the
left column is $[-0.05,1]$. Fixing $w_{1}^{*}=0$, we obtain $w_{2}\in[0.66,1]$
and $w_{3}\in[0.76,1]$ from the left and the middle columns, respectively
(following step-2).

Choosing $w_{2}^{*}=0.7\in[0.66,1]$ leads to $w_{3}\in[0.4,1]$ from the right
column plot (following step-3). Taking the intersection of the aforementioned
two ranges of $w_{3}$ results in $W_{3}=[0.76,1]$. In the left column, the
lower boundary of the parameter space for CE around $w_{1}^{*}=0$ is curved in
this case. Hence, we cannot be fully certain about $w_{3}^{\rm min}$
(following step-4). For example, we find that the point (0, 0.7, 0.8) which
appears to lie inside the projected CE region will actually not be detectable
by CE.

Example 3:

Let us now consider $E_{1}=10^{7}$ GeV and $E_{2}=10^{4}$ GeV, as shown in
Fig. 8. Focusing on detection with the ET detector (orange-shaded region)
(following step-1), the range for $w_{1}$ is $[0.12,0.97]$. Note that there is
a short interval $w_{1}\in[0.97,1]$, which does not result in a detection with
ET for any values of $w_{2}$ and $w_{3}$. This is due to the fact that for
$w_{1}\in[0.97,1]$ (independent of the values of $w_{2}$ and $w_{3}$), either
the GW signal is not large enough to be detected by ET or it violates the BBN
constraint.

Let us fix $w_{1}^{*}=0.8\in[0.12,0.97]$, leading to $w_{2}\in[0.28,1]$ (from
the left column plot) and $w_{3}\in[0.37,1]$ (from the middle column plot)
(following step-2). Choosing $w_{2}^{*}=0.6\in[0.28,1]$, the range of $w_{3}$
(from the right column plot) is found to be $[0.48,1]$ (following step-3). The
intersection of these two ranges of $w_{3}$ results in $W_{3}=[0.48,1]$. From
the left column plot in Fig. 8, we observe that lower boundary of ET region
near $w_{1}^{*}=0.8$ is again not flat. Hence, we can only conclude that for
$w_{1}^{*}=0.8$ and $w_{2}^{*}=0.6$, and not-all-values of $w_{3}\in[0.48,1]$
will result in a detection in ET (following step-4).

For example, as per step-5, the point (0.8, 0.6, 0.8) appears to fall inside
the region which is ruled out by the BBN constraint in all of three subplots
shown in Fig. 8. But we restrict ourselves from conclusively concluding
whether this combination of EoS parameters actually satisfies the BBN
constraint or not, since these are just 2-dimensional projections of the full
three-dimensional parameter space.

Since interpreting conclusively from the projected plots for three-epoch
reheating phase is quite difficult, and at times ambiguous, we provide a link
to the GitHub repository  which contains a Python code to directly plot the
spectral energy density of first-order GWs. The reader can use this as a tool
to conclusively determine whether a particular combination of EoS parameters
would result in a potentially detectable GW signal (which is consistent with
BBN and aLIGO constraints) for multiple transitions during reheating.

As mentioned in the case of two-epoch reheating, higher values of $E_{1}$
implies a shorter duration of the first epoch during reheating following the
end of inflation. As a result, the aforementioned GW detectors will be able to
impose weaker constraints on $w_{1}$ for large values of $E_{1}$, as evident
from Fig. 9. This is because if the first post-inflationary epoch is short,
then the detectability essentially relies on the subsequent epochs,
specifically on the final epoch before BBN. Similarly, higher values of
$E_{2}$ widen the range of detectable values of $w_{2}$ further, with the
range of $w_{3}$ being shortened (see Figs. 9 and 9). To conclude, our results
in this section indicate that higher values of $w_{2}$ and $w_{3}$ would
result in detectable GW signal, as expected.

##  5 Post-inflationary evolution in a scenario inspired from String Theory

In our analysis of the spectral energy density of GWs in Sec. 4, we remained
agnostic about the exact values of the post-inflationary EoS parameters and
focused on determining the sub-space of the post-inflationary parameter space
that results in a potentially detectable signal in the upcoming GW detectors,
without violating the existing constraints. In this section, we focus on a
particular post-inflationary scenario considered in a recent study [137],
which is phenomenologically inspired from String Theory [162, 163, 164, 165,
166].

In type 2B String Theory within the large volume scenario (LVS) [165], a
resultant moduli field $\chi$ (which is an artifact of string
compactification) may persist throughout the universe, potentially long enough
until the onset of BBN. During inflation, the moduli field is sub-dominant and
remains nearly frozen in its potential, while after the end of inflation, the
volume modulus comes to dominate the universe and commences rolling swiftly
down its steep exponential potential. The moduli field eventually reaches the
minimum of its potential around which it oscillates and later decays, thereby
reheating the universe to the hot Big Bang phase. Thus, the dynamics of the
universe from the end of inflation until the onset of BBN is primarily
dictated by the moduli field in this scenario [137].

The succeeding epoch after the end of inflation is therefore a kination
epoch666See Refs. [131, 132] for discussions on kination cosmology. during
which the moduli exhibits a fast-roll phase, with $\dot{\chi}^{2}\gg V(\chi)$,
down a steep region of its potential of the form

$V(\chi)\simeq V_{0}\,e^{-\lambda\,\frac{\chi}{m_{p}}}\,,$ (5.1)

where $\lambda=\sqrt{27/2}$ and $V_{0}=m_{p}^{4}$ for LVS. The EoS during this
kination epoch is
$w=\frac{\dot{\chi}^{2}-2V(\chi)}{\dot{\chi}^{2}+2V(\chi)}\simeq 1$. While the
universe, at this epoch, is also filled with a small fraction of radiative
degrees of freedom that was generated from the decay of the inflaton at the
end of inflation. Due to a sharp decline of the moduli energy density during
kination, the universe enters into a radiation-dominated epoch, following the
kination phase.

The moduli field, rolling down its steep potential in the radiation dominated
epoch, then enters into a tracker regime [167, 168, 169, 170, 171], as
illustrated in Fig. 10, where its energy density scales proportionally to the
energy density of the background universe. Eventually, the field leaves the
steep region of its potential and rolls down to its stable (quadratic) local
minimum about which it begins to oscillate. The oscillating moduli field then
dominates the energy budget of the universe again, leading to a matter-
dominated universe. The oscillating field then decays, potentially non-
perturbatively (hence efficiently), thereby reheating the universe to give
rise to the hot Big Bang phase.

Figure 10: A schematic depiction of the chronology of post-inflationary epochs
in the String-inspired model considered in Ref. [137], where the final epoch
of reheating is assumed to be instantaneous. Note that the horizontal length
is not an accurate representation of the duration of each epoch. Figure 11:
The spectral energy density ($\Omega_{{}_{\rm GW}}$) of inflationary GWs for
the String-inspired scenario as illustrated in Fig. 10, where we have an early
kination epoch ($w\simeq 1$) following the end of inflation, which was later
followed by a radiation-dominated epoch ($w=1/3$) and subsequently by a
matter-dominated epoch ($w=0$). The final epoch of reheating is assumed to be
instantaneous at an energy scale $E_{\rm r*}=1$ GeV. The energy scales $E_{i}$
(with $i\in\mathbb{Z}^{+}$) corresponding to the end of each epoch, are given
in the top label of the plot. The energy scale at the end of inflation is the
same as in our previous plots, namely, $E_{\rm inf}=5.76\times 10^{15}$ GeV,
so as to have the tensor-to-scalar ratio $r=0.001$.

In short, the specific chronology of post-inflationary epochs in this scenario
involves an early kination epoch $(w\approx 1)$ following the end of
inflation, followed by a radiation-dominated epoch $(w=1/3)$, and
subsequently, a matter-dominated epoch $(w=0)$, before reheating the universe
prior to the commencement of the BBN. Note that this final epoch of reheating
is assumed to be instantaneous. Fig. 11 depicts a typical plot of spectral
energy density against the current frequency of GWs for this particular model.

Figure 12: The parameter-space plot for the energy scale at the end of the
kination epoch ($E_{1}$), and that at the beginning of the early matter
domination ($E_{2}$), which leads to a detectable GW signal in the future
detectors (assuming the final epoch of reheating to be instantaneous) for the
case of a String-inspired scenario considered in Ref. [137] and illustrated in
Fig. 10.

Following the methodology as spelled out systematically in Sec. 4, we
determine the parameter space of the energy scales at the end of each of the
aforementioned epochs. To be consistent with our previous analyses, let us
assume the energy scale at the end of inflation to be $E_{\rm inf}=5.76\times
10^{15}$ GeV, and that at the end of (the instantaneous period of) reheating
to be $E_{\rm r*}=1$ GeV. Specifically, we determine the range of energy scale
at the end of the kination epoch ($E_{1}$) and that at the end of (the early)
radiation domination epoch ($E_{2}$), which results in detectable GW signals.
The corresponding parameter-space plot is shown in Fig. 12.

From Fig. 12, we observe that this particular scenario results in a GW signal
which will be detectable by BBO, DECIGO, and CE (for our choice of $E_{\rm
inf}$ and $E_{\rm r*}$) for a large range of $E_{1}$ and $E_{2}$. However, the
combination of $E_{1}$ and $E_{2}$, that leads to a detectable signal in LISA
and ET, actually violates the BBN constraint given in Eq. (3.59). Furthermore,
the spectrum of only those tensor modes that entered the Hubble radius during
the kination epoch are high enough to intersect the sensitivity curves of the
detectors, since the kination epoch results in a highly blue-tilted GW
spectrum. The plot also demonstrates that a detection of primordial GWs is
possible in this scenario only if $E_{2}\lesssim{\cal O}(10^{2})$ GeV, or
equivalently, if the duration of the early matter-dominated epoch is roughly
less than 7 $e$-folds, i.e., $N_{\rm matter}\lesssim 7$.

##  6 Discussion and conclusions

One of the key predictions of the inflationary paradigm is the generation of
tensor fluctuations which later constitute a stochastic background of
gravitational waves (GWs) upon their Hubble-entry after the end of inflation
[30]. These primordial GWs encode crucial information about the state of the
universe, both at the time of their production, and the subsequent epochs
through which they propagate before reaching the terrestrial and celestial
detectors. Below the Planck scale, the GWs interact minimally with matter and
energy, making them the most ideal candidate for probing the dynamics of the
very early universe [30, 58, 60, 62, 63, 59]. Numerous studies have recently
examined the feasibility of detecting these GWs using both the current and
upcoming GW detectors [119, 172, 147, 173, 174, 175, 63, 64, 176].

In this paper, we focused on the possibility of probing the unknown post-
inflationary history of the universe via primordial GWs. Since the dynamics of
reheating involves a number of complex phases during which the universe
transits from a highly non-thermal state at the end of inflation to the
thermal hot Big Bang phase [177, 93, 110], it is important to describe
reheating by a series of epochs with distinct EoS parameters [130, 104].
Hence, it is phenomenologically important to determine the possible post-
inflationary scenarios that would result in a detectable stochastic GW signal
in the upcoming GW detectors, without violating the existing constraints.
Therefore, we have dedicated this paper in computing the spectral energy
density of (first-order) inflationary GWs at the present epoch by considering
the post-inflationary universe to feature multiple successive phases, each
with a (nearly) constant EoS parameter. We determined the subspace of the EoS
parameters of different epochs, as well as the corresponding duration of each
epoch, during reheating that leads to a potentially detectable GW signal.

After a brief review of the inflationary dynamics, both at background and
perturbation levels, in Sec. 2, we focused on studying the evolution of
inflationary (first-order) tensor fluctuations in the post-inflationary
universe in Sec. 3. We began by solving the evolution equation for the tensor
mode functions, where we obtained the general solution for a given post-
inflationary epoch described by a constant EoS parameter in terms of the
Bessel functions. By assuming the transitions between successive epochs of
constant EoS to be instantaneous, we obtained expressions for the unknown
coefficients appearing in the general solutions for the tensor mode functions;
be it for the epoch immediately following the end of inflation, as given by
Eq. (3.32), or any other epoch succeeding that, as given by Eqs. (3.20) and
(3.21). We computed the spectral energy density $\Omega_{{}_{\rm GW}}(f)$ of
GWs using Eq. (3.51), and checked whether their resulting spectrum will be
detectable by the upcoming GW observatories, such as LISA, ET, CE, DECIGO and
BBO; for a variety of reheating scenarios as demonstrated in Fig. 4. In the
case of a reheating scenario consisting of two successive phases, the EoS
parameter space which results in a detectable GW signal is plotted in Fig. 6.
Similarly for a three-epoch reheating phase, the projections of the three-
dimensional EoS parameter space are shown in Figs. 7-9.

Our results primarily indicate that if inflationary GWs are detected by the
aforementioned observatories in the near future, then it is highly likely that
the temperature at the end of reheating was low enough, and that the universe
evolved through one or more epochs of stiff-matter ($w>1/3$) dominated phase
towards the end of reheating. More importantly, our analysis demonstrates that
a predicted GW signal, that could have violated the BBN and/or aLIGO
constraints when modelled by a single stiff EoS, $w_{\rm stiff}>1/3$, may
actually be consistent with the aforementioned constraints if the post-
inflationary history is comprised of additional phases with softer EoS
parameters, $w\leq 1/3$, along with the stiff-matter phase. Therefore, our
work highlights the drawback of using a single average EoS to describe the
post-inflationary expansion, and opens up a broad range of parameter space
which generates robust forecasts for the upcoming GW observatories.

Finally, in Sec. 5 we applied our techniques to a particular post-inflationary
scenario described in Ref. [137], which is phenomenologically inspired from
String Theory, in which a scalar moduli field plays a cosmologically important
role, and dictates the evolution of the universe in between the end of
inflation and the onset of BBN. In this scenario, the post-inflationary
universe transits through a series of chronological epochs, namely, a
kination-dominated epoch, an intermediate radiation-dominated epoch, and an
early matter-dominated epoch (see Fig. 10), prior to the standard radiative
hot Big Bang phase. We explored the possibility of observing the imprints of
this scenario in the inflationary GW spectrum via upcoming GW observatories.
Representing the duration of each epoch in terms of the energy scale at the
end of that epoch, we obtained the corresponding parameter space that results
in a detectable signal, as shown in Fig. 12; with the assumption that the
final stage reheating is instantaneous. We found that GWs, corresponding to
tensor modes that made their Hubble-entry during the kination phase, will be
detected by BBO, DECIGO and CE, and not by LISA and ET (due to violation of
BBN constraints). Furthermore, we obtain a necessary condition for this
scenario to result in a detectable GW signal to be that the early matter-
dominated epoch prior to the end of reheating must be shorter in duration,
lasting no more than $7\,e$-folds of expansion.

Before concluding, we emphasize that the specific analysis carried out in this
work comes with a number of caveats, which we list in the following.

  * •

Throughout this paper, we have fixed the value of the tensor-to-scalar ratio
at $r=0.001$ which is a standard target for the upcoming CMB missions [178,
179]. In this sense, our analysis was focused on the question: ‘if low-
frequency primordial GWs get detected by their imprints on the CMB B-mode
polarization in the upcoming decade, then what will the implications be for
the upcoming GW observatories functioning at higher frequencies?’. However,
one can easily consider lower values of tensor-to-scalar ratio and incorporate
our computational scheme to generate updated forecast for GWs. Nevertheless,
we provide a detailed discussion on the effect of changing reheating
temperature on the GW spectral energy density in App. F.

  * •

Similarly, the energy scale of the universe at the end of reheating was fixed
to be $E_{\rm r*}=1$ GeV in this paper in order to demonstrate our approach
and to focus on the EoS parameter and duration of each post-inflationary (pre-
BBN) epoch. However, the hot Big bang phase might have commenced at a much
lower temperature (as low as ${\cal O}(1)\,{\rm MeV}$) or even at higher
temperatures (as high as ${\cal O}(10^{15})\,{\rm GeV}$). The lower reheating
temperature is of greater phenomenological interest, and one can utilize our
scheme to generate the corresponding GW forecast. However, much higher
reheating temperature usually results in a low amplitude of GWs which may be
of phenomenological interest for the high-frequency resonant cavity GW
detectors [180, 181, 182, 183].

  * •

More importantly, while realistic transitions (no matter how sharp they are)
in the early universe are expected to be smooth [102, 119], in this work we
modelled the transitions between different post-inflationary epochs to be
instantaneous, primarily for two reasons: (i) for the convenience of solving
the tensor mode equations analytically in terms of Bessel functions, and (ii)
to remain agnostic about the exact functional form of the EoS $w(\tau)$.
However, it is possible to carry out a fully numerical approach by either
using a (phenomenologically) smoothed EoS parameter, or by considering a
specific model where the exact evolution of $w(\tau)$ is known. It has been
shown that a smoothed reheating EoS in the perturbative regime leads to
transient oscillatory features in the GW spectrum [59]. Nevertheless, since we
primarily focused on determining the detectable parameter space of post-
inflationary dynamics, the implementation of instantaneous transitions
suffices for our purpose.

  * •

In our analysis, we did not take into account the finer effects coming from
the variation in the effective number of relativistic degrees of freedom,
$g_{*}(T)$, in the early universe [116] as well as due to the anisotropic
stress induced by freely streaming relativistic neutrinos [114] on the
spectrum of GWs.

  * •

We assumed a (nearly-) flat primordial tensor power spectrum (at the end of
inflation) to be the initial condition for evolving the tensor modes post
inflation and for computing the present-day GW spectrum. However, towards the
end of inflation, the tensor power-spectrum acquires a scale dependence which
affects the GW spectrum [59] closer to the highest frequency values $f_{e}$.
Additionally, we did not consider the UV tail of GWs [184] corresponding to
tensor modes that never became super-Hubble during inflation. However, the
aforementioned effects are more relevant for detecting ultra-high frequency
GWs, potentially via resonant cavity detectors, which was not the primary
focus of our work.

  * •

Our treatment of tensor modes in this paper was based on first-order
inflationary GWs. However, at second-order in cosmological perturbation
theory, GWs can be induced by the presence of potentially large scalar
perturbations in the primordial universe [185, 43, 145, 186]. They get
amplified by the enhancement of scalar curvature perturbations at small scales
and/or by the presence of an early matter-dominated epoch in the post-
inflationary universe [187, 145], which might lead to a larger background of
second-order GWs in comparison to their first-order counterparts. The CMB
observations can probe scalar fluctuations with comoving wavenumber only in
the range around $k=0.05\,{\rm Mpc}^{-1}$, hence they are insensitive to small
scale primordial fluctuations. Large enough fluctuations at smaller
cosmological scales may also lead to the formation of copious amounts of
primordial black holes (PBHs), a plausible candidate for dark matter. Hence, a
potential detection of these induced GWs would allow us to put strong
constraints on the small-scale primordial density perturbations [188] and on
the abundance of PBHs [189, 190, 191]. We will revert to second-order GWs in a
future project.

  * •

Finally, it is plausible that non-linear solitonic and non-solitonic
configurations can contribute considerably towards the energy budget during
the early stages of the post-inflationary evolution [112, 100, 192]. They
further induce second and higher-order GWs that is expected to be detected by
the planned GW detectors. For example, high-frequency (GHz-scale) stochastic
GWs can result from oscillon formation (and inflaton fragmentation in general)
during reheating [112, 193]. However, such signals are beyond the
sensitivities of the current and proposed detectors, and are relevant for the
resonant cavity detectors [194]. Another example would be the stochastic GW
associated with cosmic strings which can be picked up by LISA, DECIGO and BBO
[133, 134, 195, 196]. Thus, it is important to classify the possible GW signal
resulting from such differing scenarios, which is outside the scope of this
paper.

##  7 Acknowledgements

SSM is supported by the STFC Consolidated Grant [ST/T000732/1] at the
University of Nottingham. AKS was supported by the INSPIRE scholarship of the
Department of Science and Technology (DST), Govt. of India, during his
Master’s thesis work. SSM is grateful to IUCAA for the hospitality during the
earliest stage of this work. SSM thanks Varun Sahni for insightful discussions
over the years which led to the inception of this project. SSM is also
thankful to Ed Copeland, Oliver Gould, Sanjit Mitra, Ranjeev Misra and P. P.
Divakaran for helpful discussions at various stages of this project. Numerical
computations were performed on the Padmanabha HPC cluster at IISER TVM. For
the purpose of open access, the authors have applied a CC BY public copyright
license to any Author Accepted Manuscript version arising.

Data Availability Statement: This work is entirely theoretical and has no
associated data. A Python code to generate the plots of spectral energy
density of first-order inflationary GWs can be found at the GitHub repository
.

## Appendices

##  Appendix A Analytical solution to the Mukhanov-Sasaki equation

The Mukhanov-Sasaki (MS) Eq. (2.9) can be written as

$v_{k}+\left(k^{2}-\frac{\nu^{2}-\frac{1}{4}}{\tau^{2}}\right)v_{k}=0\,,$
(A.1)

where $\frac{z^{\prime\prime}}{z}=(aH)^{2}(\nu^{2}-1/4)$ with
$z=am_{p}\sqrt{2\epsilon_{{}_{H}}}$. For slow-roll inflation, at early times
all relevant scales were sub-Hubble,
$(\tau\xrightarrow{}-\infty\;\;\text{or}\;\;|k\tau|\gg 1)$. In this limit, Eq.
(A.1) reduces to the form

$v_{k}+k^{2}v_{k}=0\,,$ (A.2)

which is similar to the equation for a simple harmonic oscillator with a time-
independent frequency. A unique solution to Eq. (A.2) is obtained by fixing
the vacuum in the sub-Hubble regime to be the Bunch-Davies vacuum [139, 140],
namely,

$\lim_{k\tau\rightarrow-\infty}v_{k}\simeq\frac{1}{\sqrt{2k}}\,e^{-ik\tau}\,.$
(A.3)

In terms of the variable to be $T=-k\tau$, the MS equation becomes

$\frac{\mathrm{d}^{2}v_{k}}{\mathrm{d}T^{2}}+\left(1-\frac{\nu^{2}-\frac{1}{4}}{T^{2}}\right)v_{k}=0\,.$
(A.4)

Furthermore, we also introduce the variable $v_{k}=F\sqrt{T}$, in order to
transform Eq. (A.4) to the Bessel equation

$\frac{\mathrm{d}^{2}F}{dT^{2}}+\frac{1}{T}\frac{\mathrm{d}F}{\mathrm{d}T}+\left(1-\frac{\nu^{2}}{T^{2}}\right)F=0\,,$
(A.5)

whose general solution can be written in terms of Hankel functions as

$F(T)=C_{1}\,H_{\nu}^{(1)}(T)+C_{2}\,H_{\nu}^{(2)}(T)\,.$ (A.6)

Hence, the general solution of the MS equation can be written as

$v_{k}(T)=\sqrt{T}\left[C_{1}\,H_{\nu}^{(1)}(T)+C_{2}\,H_{\nu}^{(2)}(T)\right]\,.$
(A.7)

In the sub-Hubble limit, $T\gg 1$, the Hankel functions become

$\displaystyle H_{\nu}^{(1)}(T)\Big{|}_{T\rightarrow\infty}$
$\displaystyle\simeq\sqrt{\frac{2}{\pi
T}}\,e^{iT}\,e^{-i(\nu+\frac{1}{2})\frac{\pi}{2}}\,,$ (A.8) $\displaystyle
H_{\nu}^{(2)}(T)\Big{|}_{T\rightarrow\infty}$
$\displaystyle\simeq\sqrt{\frac{2}{\pi
T}}\,e^{-iT}\,e^{i(\nu+\frac{1}{2})\frac{\pi}{2}}\,.$ (A.9)

From the Bunch-Davies initial condition in Eq. (A.3), we have

$v_{k}(T)\Big{|}_{T\rightarrow\infty}=\frac{1}{\sqrt{2k}}e^{iT}=\sqrt{T}\,C_{1}\,H_{\nu}^{(1)}(T)\Big{|}_{T\rightarrow\infty}\,,$
(A.10)

which implies that the coefficients in Eq. (A.7) are

$C_{1}=\frac{1}{\sqrt{2k}}\,\sqrt{\frac{\pi}{2}}\,e^{i(\nu+\frac{1}{2})\frac{\pi}{2}}\quad,\quad
C_{2}=0\,.$ (A.11)

Ultimately, the expression for the mode function becomes

$v_{k}(T)=\frac{1}{\sqrt{2k}}\,\sqrt{\frac{\pi}{2}}\,e^{i(\nu+\frac{1}{2})\frac{\pi}{2}}\,\sqrt{T}\,H_{\nu}^{(1)}(T)\,.$
(A.12)

##  Appendix B Expression for the conformal Hubble parameter

Consider an epoch with EoS, $w$. Let ‘$i$’ denote the beginning of that epoch,
then the energy density evolves as

$\rho(\tau)=\rho_{i}\left(\frac{a(\tau)}{a_{i}}\right)^{-3(1+w)}\,.$ (B.1)

Hence, we can write the Hubble rate during the epoch to be,
$H^{2}=H_{i}^{2}\,\left(a/a_{i}\right)^{-3(1+w)}$. Integrating the above
expression from $a_{i}$ to an arbitrary $a$ during this epoch, we obtain

$\displaystyle\int_{a_{i}}^{a}\mathrm{d}a\left(\frac{a^{(3(1+w)/2)-1}}{a_{i}^{3(1+w)/2}}\right)$
$\displaystyle=\int_{\tau_{i}}^{\tau}a_{i}H_{i}\mathrm{d}\tau\,,$ (B.2)

which leads to the expression for the scale factor $a(\tau)$ during that epoch
as

$a(\tau)=a_{i}\left[1+\frac{a_{i}H_{i}(\tau-\tau_{i})}{\alpha}\right]^{\alpha}\,,$
(B.3)

where $\alpha=2/(1+3w)$. Hence, the conformal Hubble parameter ${\cal H}$
becomes

$\displaystyle{\cal H}$ $\displaystyle=aH=\frac{\mathrm{d}\ln
a}{\mathrm{d}\tau}\,,$ $\displaystyle\ln a$ $\displaystyle=\ln
a_{i}+\alpha\,\ln{\left[1+\frac{a_{i}H_{i}(\tau-\tau_{i})}{\alpha}\right]}\,,$
$\displaystyle\Rightarrow aH$
$\displaystyle=\frac{a_{i}H_{i}}{\left[1+\frac{a_{i}H_{i}(\tau-\tau_{i})}{\alpha}\right]}\,.$
(B.4)

##  Appendix C Energy-momentum tensor of gravitational waves

The quadratic action for the tensor fluctuations $h_{ij}$ is given by

$S_{2}[h_{ij}]=\frac{m_{p}^{2}}{8}\int\mathrm{d}^{4}x\;\sqrt{-\bar{g}}\left[-\bar{g}^{\mu\nu}(\partial_{\mu}h_{ij})(\partial_{\nu}h^{ij})\right]\,,$
(C.1)

where $\bar{g}_{\mu\nu}$ is the background metric and $\bar{g}\equiv{\rm
det}\,\bar{g}_{\mu\nu}$. The energy-momentum tensor is obtained from the
variation of this action with respect to the background (inverse) metric [148]

$\delta
S=-\frac{1}{2}\int\mathrm{d}^{4}x\sqrt{-\bar{g}}\,T_{\mu\nu}\,\delta\bar{g}^{\mu\nu}\,,$
(C.2)

leading to

$\displaystyle T_{\alpha\beta}$
$\displaystyle=\frac{m_{p}^{2}}{4}\left[-\frac{1}{2}\,\bar{g}_{\alpha\beta}\,\bar{g}^{\mu\nu}\,(\partial_{\mu}h_{ij})\,(\partial_{\nu}h^{ij})+(\partial_{\alpha}h_{ij})\,(\partial_{\beta}h^{ij})\right]\,,$
$\displaystyle\tensor{T}{{}^{\mu}_{\nu}}$
$\displaystyle=\frac{m_{p}^{2}}{4}\,\left[-\frac{1}{2}\,\tensor{\delta}{{}^{\mu}_{\nu}}\,\bar{g}^{\alpha\beta}(\partial_{\alpha}h_{ij})\,(\partial_{\beta}h^{ij})+(\partial^{\mu}h_{ij})\,(\partial_{\nu}h^{ij})\right]\,.$
(C.3)

The time-time component of the energy-momentum tensor is the energy density of
GWs, given by

$\rho_{{}_{\rm
GW}}(\tau,\vec{x})=-\tensor{T}{{}^{0}_{0}}=\frac{m_{p}^{2}}{8a^{2}(\tau)}\left[(h^{\prime}_{ij}(\tau,\vec{x}))^{2}+(\gradient
h_{ij}(\tau,\vec{x}))^{2}\right]\,.$ (C.4)

##  Appendix D Derivation of the present-day GW spectral energy density

The derivation is done for the general case of reheating consisting of
multiple equation of state parameters. For $\tau=\tau_{0}$ in Eq. (3.41)

$\Omega_{{}_{\text{GW}}}(\tau_{0},k)=\frac{k^{2}}{12\,a^{2}(\tau_{0})H^{2}(\tau_{0})}\;{\cal
P}_{T}(\tau_{0},k)\,.$ (D.1)

Substituting $k=y_{\rm eq}\,a_{\rm eq}\,H_{\rm eq}$ and ${\cal
P}_{T}(\tau_{0},k)=(k^{3}/\pi^{2})\,\overline{|h_{k,\,{\rm MD}}^{\lambda}(y\gg
1)|^{2}}$ in Eq. (D.1), we obtain

$\displaystyle\Omega_{{}_{\rm GW}}(\tau_{0},k)$
$\displaystyle=\frac{1}{12\,a_{0}^{2}H_{0}^{2}}\,\left(y_{\rm eq}\,a_{\rm
eq}\,H_{\rm
eq}\right)^{2}\,\frac{k^{3}}{\pi^{2}}\,\frac{1}{(2\,y_{0})^{4}}\,\frac{1}{\pi}\left[|A_{k,\,{\rm
MD}}|^{2}+|B_{k,\,{\rm MD}}|^{2}\right]\,,$ (D.2)

where $\overline{|h_{k,\,{\rm MD}}^{\lambda}(y\gg 1)|^{2}}$ was substituted
from Eq. (3.43). We can rewrite $y_{0}$ as

$\displaystyle y_{0}$ $\displaystyle=\frac{k}{a_{0}H_{0}}=y_{\rm
eq}\,\frac{a_{\rm eq}H_{\rm eq}}{a(\tau_{0})H(\tau_{0})}=y_{\rm
eq}\,\left[1+\frac{a_{\rm eq}H_{\rm eq}(\tau_{0}-\tau_{\rm eq})}{2}\right]\,,$
$\Rightarrow y_{0}=y_{\rm eq}\,\left(\frac{a(\tau_{0})}{a_{\rm
eq}}\right)^{1/2}\,.$ (D.3)

Incorporating Eq. (D.3) into Eq. (D.2), we obtain

$\displaystyle\Omega_{{}_{\rm GW}}(\tau_{0},k)$
$\displaystyle=\frac{1}{192}\frac{a^{4}_{\rm eq}H^{2}_{\rm
eq}}{a_{0}^{4}H_{0}^{2}}\,\frac{k^{3}}{\pi^{3}}\frac{1}{y_{\rm
eq}^{2}}\left[\tilde{A}^{2}_{k,\,{\rm MD}}+\tilde{B}^{2}_{k,\,{\rm
MD}}\right]\,\left|h_{k,\,{\rm inf}}^{\lambda}\right|^{2}\,,$
$\displaystyle=\frac{1}{96\pi^{3}}\frac{a^{4}_{\rm r*}H^{2}_{\rm
r*}}{a_{0}^{4}H_{0}^{2}}\,\frac{1}{y_{\rm
eq}^{2}}\left[\tilde{A}^{2}_{k,\,{\rm MD}}+\tilde{B}^{2}_{k,\,{\rm
MD}}\right]\,\frac{H^{2}_{\rm
inf}}{m_{p}^{2}}\left(\frac{k}{k_{*}}\right)^{n_{{}_{T}}}\,,$ (D.4)

where $\tilde{A}_{k,\,{\rm MD}}=\frac{A_{k,\,{\rm MD}}}{h_{k,{\rm
inf}}^{\lambda}}$ and similarly, $\tilde{B}_{k,\,{\rm MD}}=\frac{B_{k,\,{\rm
MD}}}{h_{k,{\rm inf}}^{\lambda}}$. Furthermore, we have made use of the
following relations:

${\cal P}_{T,\,{\rm
inf}}(\tau_{0},\,k)=\frac{k^{3}}{\pi^{2}}\,\left|h_{k,\,{\rm
inf}}^{\lambda}\right|^{2}=\frac{2}{\pi^{2}}\left(\frac{H_{\rm
inf}}{m_{p}}\right)^{2}\,\left(\frac{k}{k_{*}}\right)^{n_{{}_{T}}}\,,$ (D.5)
$a^{4}_{\rm eq}H^{2}_{\rm eq}=(a_{\rm eq}H_{\rm eq})^{2}a^{2}_{\rm
eq}=\left[\frac{a_{\rm r*}H_{\rm r*}}{1+a_{\rm r*}H_{\rm r*}(\tau_{\rm
eq}-\tau_{\rm r*})}\right]^{2}a_{\rm r*}^{2}[1+a_{\rm r*}H_{\rm r*}(\tau_{\rm
eq}-\tau_{\rm r*})]^{2}=a_{\rm r*}^{4}H_{\rm r*}^{2}\,.$ (D.6)

where $H_{{\rm inf}}$ corresponds to the Hubble scale during inflation at the
Hubble-exit time of the CMB pivot scale, i.e., $H_{{\rm inf}}=H_{*}$. The
Hubble parameter can be written as

$H^{2}(\tau)=\frac{\rho_{\rm rad}}{3\,m_{p}^{2}\,\Omega_{\rm
rad}}=\frac{\rho_{\rm rad}\,\Omega_{\rm rad,\,0}\,H_{0}^{2}}{\rho_{\rm
rad,\,0}\,\Omega_{\rm rad}}\,.$ (D.7)

Radiation energy density is given by

$\rho_{\rm rad}(T)=\frac{\pi^{2}}{30}\,g_{*}(T)\,T^{4}\,,$ (D.8)

while the conservation of entropy implies

$T\propto g_{s}^{-1/3}a^{-1}\,.$ (D.9)

Substituting Eq. (D.8) and Eq. (D.9) in Eq. (D.7), we obtain

$H^{2}(T)=\frac{g_{*}(T)\;g_{s}^{-4/3}(T)\;a^{-4}(T)\;\Omega_{\rm
rad,\,0}\;H_{0}^{2}}{g_{*,\,0}\;g_{s,\,0}^{-4/3}a_{0}^{-4}\Omega_{\rm
rad}(T)}\,,$ (D.10)

which can be rewritten as

$\frac{a_{k}^{4}H_{k}^{2}}{a_{0}^{4}H_{0}^{2}}=\frac{g_{*,k}}{g_{*,\,0}}\left(\frac{g_{s,\,0}}{g_{s,k}}\right)^{4/3}\frac{\Omega_{\rm
rad,\,0}}{\Omega_{\rm rad}(\tau_{k})}\,,$ (D.11)

where $g_{*,\,k},\,g_{s,\,k}$ are the number of relativistic degrees of
freedom in energy density and entropy density, respectively, at the time when
the tensor mode with comoving wavenumber $k$ makes its Hubble-entry, while
$g_{*,\,0},\,g_{s,\,0}$ correspond to the same at the present epoch. When the
tensor mode of interest enters during the RD epoch, $\Omega_{\rm rad}\approx
1$. So, (D.11) becomes

$\frac{a_{k}^{4}H_{k}^{2}}{a_{0}^{4}H_{0}^{2}}\approx\frac{g_{*,k}}{g_{*,\,0}}\left(\frac{g_{s,\,0}}{g_{s,k}}\right)^{4/3}\,\Omega_{\rm
rad,\,0}\,,$ (D.12)

where the label ‘$k$’ refers to any mode that entered the Hubble radius during
the RD epoch. Using Eq. (D.12) in Eq. (D.4) at $\tau=\tau_{\rm r*}$ (beginning
of RD epoch), we arrive at the final expression for the spectral energy
density of GWs

$\displaystyle\Omega_{{}_{\rm
GW}}(\tau_{0},k)=\frac{1}{96\pi^{3}}\,\frac{g_{*,\,{\rm
r*}}}{g_{*,\,0}}\left(\frac{g_{s,\,0}}{g_{s,\,{\rm
r*}}}\right)^{4/3}\,\Omega_{\rm rad,\,0}\,\frac{1}{y_{\rm
eq}^{2}}\left[\tilde{A}^{2}_{k,\,{\rm MD}}+\tilde{B}^{2}_{k,\,{\rm
MD}}\right]\,\frac{H^{2}_{\rm
inf}}{m_{p}^{2}}\left(\frac{k}{k_{*}}\right)^{n_{{}_{T}}}\,.$ (D.13)

##  Appendix E Energy density of photons and neutrinos

From the present-day CMB temperature, $T_{0}=2.73\,{\rm K}$ [197], we obtain
the energy density of CMB photons

$\rho_{\gamma,\,0}=\frac{\pi^{2}}{30}\times 2\times T_{0}^{4}=2.00\times
10^{-51}\,{\rm GeV}^{4}\,,$ (E.1)

where the factor ‘2’ accounts for the internal degrees of freedom of photons.
The present-day critical density of the universe is

$\rho_{c,\,0}=3\,H_{0}^{2}\,m_{p}^{2}=8.12\times 10^{-47}\,h^{2}\,{\rm
GeV}^{4}\,,$ (E.2)

where $H_{0}=100\,h\,{\rm km}\,{\rm s}^{-1}\,{\rm Mpc}^{-1}=2.132\times
10^{-42}\,h\,{\rm GeV}$ [25]. Hence, the density parameter of CMB photons is
given by

$h^{2}\,\Omega_{\gamma,\,0}=h^{2}\,\frac{\rho_{\gamma,\,0}}{\rho_{c,\,0}}=2.46\times
10^{-5}\,.$ (E.3)

Post electron-positron annihilation, the neutrino temperature is given by [23]

$T_{\nu}=\left(\frac{4}{11}\right)^{1/3}\,T_{\gamma}\,,$ (E.4)

from which, we can derive a relation between the energy density of SM
neutrinos and CMB photons to be

$\rho_{\nu}=\frac{7}{8}\,N_{\rm
eff}\,\left(\frac{4}{11}\right)^{4/3}\,\rho_{\gamma}\,,$ (E.5)

where $N_{\rm eff}=3.046$ [198] is the effective number of neutrino species in
the SM. Using Eq. (E.5), the present-day density parameter for neutrinos
becomes

$h^{2}\,\Omega_{\nu,\,0}=\frac{7}{8}\,N_{\rm
eff}\,\left(\frac{4}{11}\right)^{4/3}\,h^{2}\,\Omega_{\gamma,\,0}=1.70\times
10^{-5}\,.$ (E.6)

Hence, the present-day density parameter of radiation becomes

$h^{2}\,\Omega_{{\rm
rad},\,0}=h^{2}\,(\Omega_{\gamma,\,0}+\,\Omega_{\nu,\,0})=4.16\times
10^{-5}\,.$ (E.7)

##  Appendix F Effect of varying reheating temperature on GW spectral energy
density

In the main text of this paper, we had fixed the reheating temperature to be
$T_{\rm r*}=0.45\,{\rm GeV}$, or equivalently, the energy scale at the end of
reheating to be $E_{\rm r*}=1\,{\rm GeV}$, for the purpose of illustration.
However, in general, the reheating energy scale can range from the end of
inflation all the way down to the onset of BBN. Hence, in this appendix, we
discuss the impact of varying the reheating temperature on the spectral energy
density of inflationary GWs.

Additionally, in the main text of the paper, we had fixed the tensor-to-scalar
ratio during inflation to be $r=0.001$, or equivalently, the energy scale at
the end of inflation to be $E_{\rm inf}=5.76\times 10^{15}\,{\rm GeV}$,
keeping in mind that the upcoming CMB detectors will be able to probe $r$ up
to ${\cal O}\left(10^{-3}\right)$. Since inflation does not make a prediction
for $r$, it is possible that $r\ll 10^{-3}$. Therefore, we proceed to first
discuss the impact of small $r$ on the detectability of inflationary GWs via
the upcoming GW observatories, before moving on to determine the detectable
parameter space of equation of states for fixed $r=0.001$, but different
values of the reheating energy scales $E_{\rm r*}$.

###  F.1 Effect of varying tensor-to-scalar ratio on GW signal

A quantity of central importance is the minimum value of reheating temperature
corresponding to a given value of tensor-to-scalar ratio which leads to a
detectable GW signal in various GW detectors. In fact, for a given GW
detector, and for a fixed value of $r$, we find that upon increasing the
reheating temperature, there is a critical value of $T_{\rm r*}$ (or $E_{\rm
r*}$) above which the inflationary GW signal will lie outside the sensitivity
region of that detector, even if the EoS during reheating is that of a stiff-
matter, i.e. $w\simeq 1$.

In Fig. 13, we show the (shaded) region in the parameter space of $r$ and
$T_{\rm r*}$ (or $E_{\rm r*}$) for different GW detectors that will never lead
to a detectable inflationary GW signal. This does not necessarily mean that
the unshaded region corresponding to each detector will guarantee a detection
in that detector, which of course depends upon the EoS and duration of
reheating, as well as on whether the GW spectrum satisfies the current
observational constraints.

Figure 13: The figure illustrates the effect of varying tensor-to-scalar ratio
($r$), and reheating temperature ($T_{\rm r*}$) on the detectability of
inflationary GWs in various GW detectors. The coloured shaded regions indicate
the parameter space of the $r$, and $T_{\rm r*}$, that can never yield a
potential detection of inflationary GW signal in CE, ET, BBO, LISA and DECIGO
detectors. The grey-shaded region is ruled out by the CMB observations,
namely, $r\leq 0.036$. The Hubble parameter during (at the end of) inflation
($H_{\rm inf}$) corresponding to each value of $r$ and the energy scale
corresponding to each value of $T_{\rm r*}$ are provided in the respective
opposite axes.

For large enough values of $r\gtrsim{\cal O}\left(10^{-2}\right)$, height of
flat (nearly scale-invariant) region of $\Omega_{{}_{\rm GW}}(f)$,
corresponding to modes making their Hubble-entry during the thermal RD epoch,
increases enough to reach the sensitivity curves of BBO and DECIGO, as long as
the reheating temperature is also large enough. In such cases, at large
reheating temperatures, these detectors will detect a GW signal with
certainty. That is why we see nearly vertical curves for BBO and DECIGO at
high values of $r$ in Fig. 13 (blue and green curves).

However, the same is not true for other GW detectors, because their
sensitivity curves are above the highest allowed height of $\Omega_{{}_{\rm
GW}}$, corresponding to $r=0.036$, for modes that entered in the RD epoch.
Therefore, if $T_{\rm r*}$ is high enough, then the GW signal will be absent
in these detectors, independent of the EoS during reheating.

Note that the small kinks seen in the solid curves are due to the sudden drop
in the relativistic degrees of freedom in the early Universe, within the
Standard Model of particle physics. In particular, the prominent horizontal
shift in all curves around $T_{\rm r*}\simeq 150$ MeV is due to a large drop
in relativistic degrees of freedom closer to QCD phase transition.

###  F.2 Two-epoch reheating with different reheating temperature

In Sec. 4, we determined the parameter space of EoS $w_{1}$ and $w_{2}$, that
results in a potentially detectable inflationary GW signal in future GW
detectors, for fixed values $E_{\rm r*}=1\,{\rm GeV}$ and $r=0.001\,(E_{\rm
inf}=5.76\times 10^{15}\,{\rm GeV})$. As stressed before, the aforementioned
values of $r$ and reheating temperature were used purely for the purpose of
illustrating how to to use our calculations of $\Omega_{{}_{\rm GW}}$ to check
whether a particular set of parameters leads to a detectable signal.

However, for a given value of $r$, reheating temperature can be very
different, ranging in between $1\,{\rm MeV}\leq T_{\rm r*}\lesssim E_{\rm
inf}$. Therefore, in this appendix, we include the parameter space plots of
$w_{1}$ and $w_{2}$ corresponding to a two-epoch reheating scenario (similar
to Fig. 6), specifically for $E_{\rm r*}=100\,{\rm GeV}$ and $10\,{\rm MeV}$
in Figs. 14 and 15, respectively. We fix the value of $E_{\rm inf}$ to be the
same as before. Note that the BBN constraint that we have imposed in this
appendix is slightly weaker than the one quoted in Eq. (3.59). Here, we only
make sure whether the GW spectrum intersects the BBN constraint curve at
$h^{2}\,\Omega_{{}_{\rm GW}}(\tau_{0})=1.13\times 10^{-6}$.

Figure 14: The parameter space for two-epoch reheating phase with EoS
parameters $w_{1}$ and $w_{2}$ that leads to a potentially detectable GW
signal in the LISA, BBO, DECIGO, CE, and ET detectors for $E_{\rm r*}=100$ GeV
and $r=0.001$. The grey-shaded region in each plot corresponds to the
combination of $w_{1}$ and $w_{2}$ that crosses the BBN constraint curve at
$h^{2}\,\Omega_{{}_{\rm GW}}(\tau_{0})=1.13\times 10^{-6}$, while the maroon-
shaded region represents those that are ruled out by the aLIGO null detection.
Energy scale of the universe at the end of the first post-inflationary epoch
is labelled above in each plot. The duration of each epoch (number of
$e$-folds) has been provided opposite to the corresponding EoS axis.

Figure 15: The parameter space for two-epoch reheating phase with EoS
parameters $w_{1}$ and $w_{2}$ that leads to a potentially detectable GW
signal in the LISA, BBO, DECIGO, CE, and ET detectors for $E_{\rm r*}=10$ MeV
and $r=0.001$. The grey-shaded region in each plot corresponds to the
combination of $w_{1}$ and $w_{2}$ that crosses the BBN constraint curve at
$h^{2}\,\Omega_{{}_{\rm GW}}(\tau_{0})=1.13\times 10^{-6}$, while the maroon-
shaded region represents those that are ruled out by the aLIGO null detection.
Energy scale of the universe at the end of the first post-inflationary epoch
is labelled above in each plot. The duration of each epoch (number of
$e$-folds) has been provided opposite to the corresponding EoS axis.

## References

  * [1] Alexei A. Starobinsky  “A New Type of Isotropic Cosmological Models Without Singularity”  In _Phys. Lett. B_ 91, 1980, pp. 99–102  DOI: 10.1016/0370-2693(80)90670-X
  * [2] Alan H. Guth  “The Inflationary Universe: A Possible Solution to the Horizon and Flatness Problems”  In _Phys. Rev. D_ 23, 1981, pp. 347–356  DOI: 10.1103/PhysRevD.23.347
  * [3] Andrei D. Linde  “A New Inflationary Universe Scenario: A Possible Solution of the Horizon, Flatness, Homogeneity, Isotropy and Primordial Monopole Problems”  In _Phys. Lett. B_ 108, 1982, pp. 389–393  DOI: 10.1016/0370-2693(82)91219-9
  * [4] Andreas Albrecht and Paul J. Steinhardt  “Cosmology for Grand Unified Theories with Radiatively Induced Symmetry Breaking”  In _Phys. Rev. Lett._ 48, 1982, pp. 1220–1223  DOI: 10.1103/PhysRevLett.48.1220
  * [5] Andrei D. Linde  “Chaotic Inflation”  In _Phys. Lett. B_ 129, 1983, pp. 177–181  DOI: 10.1016/0370-2693(83)90837-7
  * [6] Andrei D. Linde  “Particle physics and inflationary cosmology”, 1990  arXiv:hep-th/0503203
  * [7] William H. Kinney  “TASI Lectures on Inflation”, 2009  arXiv:0902.1529 [astro-ph.CO]
  * [8] Jerome Martin, Christophe Ringeval and Vincent Vennin  “Encyclopædia Inflationaris”  In _Phys. Dark Univ._ 5-6, 2014, pp. 75–235  DOI: 10.1016/j.dark.2014.01.003
  * [9] Daniel Baumann  “Inflation”  In _Theoretical Advanced Study Institute in Elementary Particle Physics: Physics of the Large and the Small_ , 2011, pp. 523–686  DOI: 10.1142/9789814327183_0010
  * [10] Daniel Baumann  “Primordial Cosmology”  In _PoS_ TASI2017, 2018, pp. 009  DOI: 10.22323/1.305.0009
  * [11] Hideo Kodama and Misao Sasaki  “Cosmological Perturbation Theory”  In _Progress of Theoretical Physics Supplement_ 78, 1984, pp. 1–166  DOI: 10.1143/PTPS.78.1
  * [12] Antonio Riotto  “Inflation and the theory of cosmological perturbations”  In _ICTP Lect. Notes Ser._ 14, 2003, pp. 317–413  arXiv:hep-ph/0210162
  * [13] John Ellis and David Wands  “Inflation (2023)”, 2023  arXiv:2312.13238 [astro-ph.CO]