DESY-23-075

Reheating after Inflaton Fragmentation

Marcos A. G<EMAIL_ADDRESS>and Mathias
<EMAIL_ADDRESS>
♠ Departamento de Física Teórica, Instituto de Física, Universidad Nacional
Autónoma de México, Ciudad de México C.P. 04510, Mexico
♣ Deutsches Elektronen-Synchrotron DESY, Notkestr. 85, 22607 Hamburg, Germany

ABSTRACT

In the presence of self-interactions, the post-inflationary evolution of the
inflaton field is driven into the non-linear regime by the resonant growth of
its fluctuations. The once spatially homogeneous coherent inflaton is
converted into a collection of inflaton particles with non-vanishing momentum.
Fragmentation significantly alters the energy transfer rate to the inflaton’s
offspring during the reheating epoch. In this work we introduce a formalism to
quantify the effect of fragmentation on particle production rates, and
determine the evolution of the inflaton and radiation energy densities,
including the corresponding reheating temperatures. For an inflaton potential
with a quartic minimum, we find that the efficiency of reheating is
drastically diminished after backreaction, yet it can lead to temperatures
above the big bang nucleosynthesis limit for sufficiently large couplings. In
addition, we use a lattice simulation to estimate the spectrum of induced
gravitational waves, sourced by the scalar inhomogeneities, and discuss
detectability prospects. We find that a Boltzmann approach allows to
accurately predict some of the main features of this spectrum.

June 2023

###### Contents

  1. 1 Introduction
  2. 2 Post-inflationary dynamics in a quartic potential
    1. 2.1 Coherent oscillations
    2. 2.2 Parametric resonance
    3. 2.3 Backreaction and fragmentation
    4. 2.4 The Boltzmann approximation
  3. 3 Reheating in a quartic potential
    1. 3.1 Production rates
    2. 3.2 Reheating temperatures
  4. 4 Gravitational waves
    1. 4.1 Simulated spectrum
    2. 4.2 Signal and detection prospects
  5. 5 Summary and conclusions
  6. Acknowledgments

##  1 Introduction

Ever since the inception of cosmic inflation as a potential solution to the
initial condition problems of standard Big Bang cosmology, the question of how
the universe transitions from a cold, empty, quasi-de Sitter state to a
radiation dominated stage in thermal equilibrium, has been an active research
topic. The original inflationary proposal (old inflation) was in particular
abandoned due to its impossibility to lead to a phenomenologically successful
reheating of the universe [1, 2, 3, 4]. The overwhelming majority of the
subsequent proposals, based on the slow-roll of an elementary scalar field
called the inflaton, incorporate reheating mechanisms tied to the coherent
oscillation of this inflaton field about its minimum.111The assumption of a
dominant classical, coherent component to the inflaton is also present in
oscillationless reheating mechanisms such as instant preheating [5, 6].

At the perturbative (Boltzmann) level, reheating is modeled as the dissipation
of the oscillating inflaton into elementary particles, with a decay rate
determined by the quantum mechanical transition amplitude from the time-
dependent vacuum to the corresponding particle states. This amplitude is
evaluated by averaging the transition rate over the fast oscillation of the
inflaton about its minimum, assuming therefore that particle production
evolves adiabatically following the expansion of the universe. For two-body
decays, the dissipation rate is mostly sensitive to the time-dependence of the
effective mass of the oscillating inflaton, which can lead to time-dependent
kinematic blocking effects, and in general to a decreasing (or increasing)
decay efficiency depending on the shape of the minimum of the inflaton
potential [7, 8, 9, 10, 11, 12].

It is known, however, that the perturbative picture is insufficient to
describe the dynamics of reheating in the presence of strong couplings, or
when the short time-scale oscillation of the inflaton field leads to the
resonant enhancement (bosons) or suppression (fermions) of the growth of the
quantum fields associated with the decay products. When active, these
collective effects, known as preheating, generically lead to a significant
modification of the reheating dynamics [13, 14, 15, 7, 16, 17, 18]. In
scenarios with a quadratic inflaton potential, the parametric resonance of the
momentum modes of the produced fields can lead to the exponential growth of
perturbations in the early stages of reheating. This growth is transient and
disordered, due to a mismatch between the particle production rates and the
expansion rate of the universe. Nevertheless, if the resonance is sufficiently
strong and maintained for a sufficient amount of time, it can bring the growth
of fluctuations into the non-linear regime. The inflaton condensate is
fragmented in favor of a quasi-thermal bath of inflaton particles and its
decay products [19, 20, 21, 22, 23]. In some scenarios, the result is that the
Boltzmann approximation correctly describes the production of daughter fields
only in a very narrow range of the available parameter space [24].

Preheating in non-quadratic minima has also been studied extensively [25, 26,
27, 21, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37]. Notably, the presence of
inflaton self-interactions, even if small, can accumulate over several
oscillations and source the resonant growth of non-zero momenta inflaton
modes, eventually backreacting with the homogeneous component. This self-
fragmentation of the inflaton not only re-distributes the initial energy
density of the field, but it can also lead to the formation of localized,
soliton-like objects such as oscillons [38, 39, 40, 28, 31, 32]. The large
inhomogeneities triggered by the self-fragmentation of the inflaton can also
efficiently source gravitational wave (GW) production. The typical frequency
for the resulting stochastic gravitational wave background is larger than its
counterpart sourced during inflation $f\gtrsim\,\text{MHz}$ [41, 42, 43]. Such
high frequency is beyond the reach of current future earth or space based
interferometers such as LISA or the DECIGO but might be in the reach of future
dedicated experiments [44]. The resulting gravitational wave spectrum would
carry precious information about the dynamics of the inflaton in the first
instants following inflation and shed light on the reheating
epoch.222Consequences of the presence of a scalar spectator field on the GW
spectrum has been explored in several works such as Refs. [45, 33]

The depletion of the coherent condensate eventually shuts down the parametric
resonance. Following fragmentation, one must therefore follow the population
of the relativistic bath of decay products by perturbative means, accounting
for both the inflaton zero and non-zero mode dissipation. The decay rates will
differ from those prior to the backreaction epoch, significantly modifying any
estimates for the reheating temperature of the universe. A formalism to
estimate the evolution of the radiation energy density during reheating in the
pre- and post-fragmentation regimes, including the reheating temperature is
the main goal of this work. To showcase the differences with respect to the
purely perturbative evolution of the inflaton-radiation sector we specialize
to the case of a quartic inflaton potential near its minimum.

This paper is organized as follows. In Section 2 we discuss in detail the
evolution of an inflaton field in a quartic potential during reheating, in the
absence of interactions with other fields. Sec. 2.1 is devoted to the dynamics
of the background evolution of the homogeneous inflaton condensate. In Sec.
2.2 we study the parametric growth of inflaton fluctuations at linear order,
and in Sec. 2.3 we study the non-linear regime of their growth, including
backreaction and fragmentation effects. For Sec. 2.4 we include the
perturbative Boltzmann determination of the spectrum of inflaton fluctuations.
The inflaton-matter/radiation couplings are introduced in Section 3. The
analysis is divided into the study of the dissipation of the coherent inflaton
condensate (Sec. 3.1.1) and the decay of the inflaton particles after
fragmentation (Sec. 3.1.2). Our main results, the corresponding reheating
temperatures, are presented in Sec. 3.2. Finally, in Sec. 4 we determine the
induced gravitational waves from the non-linear scalar dynamics, and discuss
their potential for detectability. Our conclusions are presented in Sec. 5.

##  2 Post-inflationary dynamics in a quartic potential

The current measurements of the primordial curvature power spectrum, and
bounds on the tensor power spectrum, are compatible with the presence of a
single, slowly rolling neutral scalar field which drives inflation. We denote
this inflaton by $\phi$, and for it we will assume the following form for the
action,

$\mathcal{S}\;=\;\int\mathop{}\\!\mathrm{d}^{4}x\,\sqrt{-g}\left[\frac{1}{2}(\partial_{\mu}\phi)^{2}-V(\phi)+\mathcal{L}_{\rm
int}\right]\,.$ (2.1)

Here $g\equiv\det(g_{\mu\nu})$ is the metric determinant of a flat Friedmann-
Robertson-Walker metric with scale factor $a$, $V(\phi)$ denotes the inflaton
potential, and $\mathcal{L}_{\rm int}$ denotes its interactions with the rest
of the (extended) Standard Model. The potential $V(\phi)$ must be chosen in
order to match the measured amplitude and tilt of the scalar power spectrum,
$A_{S*}\simeq 2.1\times 10^{-9}$, $n_{s}\simeq 0.966$, and to avoid the upper
bound on the tensor-to-scalar ratio, $r<0.036$ [46, 47, 48]. Moreover, in this
work we impose the condition that this potential is quartic during the post-
inflationary reheating. Among several suitable candidates, we choose for
definiteness the quartic T-model attractor [49, 12],

$\displaystyle V(\phi)\;$ $\displaystyle=\;\lambda
M_{P}^{4}\left[\sqrt{6}\tanh\left(\frac{\phi}{\sqrt{6}M_{P}}\right)\right]^{4}$
(2.2) $\displaystyle\simeq\;\lambda\phi^{4}\,,\qquad(\phi\ll M_{P})$ (2.3)

where $M_{P}=1/\sqrt{8\pi G_{N}}\simeq 2.435\times 10^{18}\,{\rm GeV}$ denotes
the reduced Planck mass. The constant $\lambda$ is determined by the scalar
spectrum amplitude, evaluated at the horizon exit time of the Planck pivot
scale, $k_{*}=0.05\,{\rm Mpc}^{-1}$. For a generic potential, the number of
$e$-folds before the end of inflation at horizon exit time are determined from
the expression [50, 51]

$\displaystyle N_{*}\;=\;$
$\displaystyle\ln\left[\frac{1}{\sqrt{3}}\left(\frac{\pi^{2}}{30}\right)^{1/4}\left(\frac{43}{11}\right)^{1/3}\frac{T_{0}}{H_{0}}\right]-\ln\left(\frac{k_{*}}{a_{0}H_{0}}\right)-\frac{1}{12}\ln
g_{\rm reh}$
$\displaystyle\qquad+\frac{1}{4}\ln\left(\frac{V_{*}^{2}}{M_{P}^{4}\rho_{\rm
end}}\right)+\frac{1-3w_{\rm int}}{12(1+w_{\rm int})}\ln\left(\frac{\rho_{\rm
rad}}{\rho_{\rm end}}\right)\,.$ (2.4)

Here $H_{0}=67.36\,{\rm km}\,{\rm s}^{-1}{\rm Mpc}^{-1}$ [46],
$T_{0}=2.7255\,{\rm K}$ [52] and $a_{0}=1$ denote the present Hubble
parameter, photon temperature and scale factor, respectively. The energy
density at the end of inflation is denoted by $\rho_{\rm end}$, and the energy
density at the beginning of the radiation dominated era by $\rho_{\rm rad}$.
The effective number of degrees of freedom during reheating is denoted by
$g_{\rm reh}$. The $e$-fold averaged equation of state parameter during
reheating corresponds to

$w_{\rm int}\;\equiv\;\frac{1}{N_{\rm rad}-N_{\rm end}}\int_{N_{\rm
end}}^{N_{\rm rad}}w(n)\,\mathop{}\\!\mathrm{d}n\,.$ (2.5)

As is well know, and we show explicitly in the following section, for quartic
reheating $w\simeq 1/3$. The last term in (2) vanishes, and therefore the
number of $e$-folds is uniquely determined, $N_{*}\simeq 56$ for T-model
inflation, and so is the coupling constant $\lambda\simeq 3.3\times 10^{-12}$.
We now study the post-inflationary evolution of the inflaton sector, from the
coherent oscillation stage, to the backreaction and post-fragmentation phase.

###  2.1 Coherent oscillations

We first consider the dynamics of the inflaton field at the background level
after the end of inflation. Inflation ends when $\ddot{a}=0$, or equivalently
when $\dot{\phi}^{2}=V(\phi)$, which for the T-model (2.2) corresponds to

$\phi_{\rm end}\;\simeq\;1.52\,M_{P}\,,\quad\rho_{\rm end}\;\simeq\;(4.5\times
10^{15}\,{\rm GeV})^{4}\,.$ (2.6)

Our focus here will be the early stage of reheating, with $a/a_{\rm
end}\lesssim\mathcal{O}(10^{2})$. Importantly, we will assume that the main
decay channel of the inflaton is to fermionic fields, which do not manifest
the exponential growth from the parametric resonance of their mode equations.
This allows us to assume that the coherence of the inflaton oscillation is
maintained, and that the third term in (2.1) can be disregarded at early times
[11, 12]. If this is the case, variation of the action with respect to the
homogeneous inflaton field and metric yields the equations of motion

$\displaystyle\ddot{\phi}+3H\dot{\phi}+4\lambda\phi^{3}\;$
$\displaystyle\simeq\;0\,,$ (2.7)
$\displaystyle\frac{1}{2}\dot{\phi}^{2}+\lambda\phi^{4}\;$
$\displaystyle\simeq\;3M_{P}^{2}H^{2}\,,$ (2.8)

where $H=\dot{a}/a$ is the Hubble parameter, and a dot represents
differentiation with respect to cosmic time. These equations describe an
underdamped anharmonic oscillator, which can be parametrized in terms of an
envelope function $\phi_{0}(t)$, encoding the redshift due to expansion, and a
quasi-periodic function $\mathcal{P}(t)$, which encodes the short time-scale
oscillation,

$\phi(t)\;\simeq\;\phi_{0}(t)\,\mathcal{P}(t)\,.$ (2.9)

To determine the time-dependence of the envelope function, we multiply Eq.
(2.7) by $\phi$ and average over one oscillation. This yields,

$\langle\phi\ddot{\phi}+3H\phi\dot{\phi}+4\lambda\phi^{4}\rangle\;\simeq\;-\langle\dot{\phi}^{2}\rangle+4\langle
V(\phi)\rangle\;=\;0\,,$ (2.10)

implying that the oscillation-averaged energy and pressure densities can be
written as [12]

$\displaystyle\rho_{\phi}\;$
$\displaystyle\equiv\;\frac{1}{2}\langle\dot{\phi}^{2}\rangle+\langle
V(\phi)\rangle\;\simeq\;3\langle V(\phi)\rangle\;=\;V(\phi_{0})\,,$ (2.11)
$\displaystyle p_{\phi}\;$
$\displaystyle\equiv\;\frac{1}{2}\langle\dot{\phi}^{2}\rangle-\langle
V(\phi)\rangle\;\simeq\;\langle V(\phi)\rangle\;=\;\frac{1}{3}V(\phi_{0})\,,$
(2.12)

where we have used $\langle\mathcal{P}^{4}\rangle=\frac{1}{3}$, as it can be
checked from the solution for $\mathcal{P}$, see (2.15) below. Hence, the
equation of state parameter is $w_{\phi}\simeq 1/3$, and the equation of
motion (2.7) can be rewritten as

$\dot{\rho}_{\phi}+4H\rho_{\phi}\;\simeq\;0\,.$ (2.13)

Thus, as is well known, the energy density of the inflaton redshifts as
radiation, $\rho_{\phi}\propto a^{-4}$. As a result, we determine that the
decaying envelope of the oscillation redshifts as

$\phi_{0}(t)\;\simeq\;\phi_{\rm end}\left(\frac{a_{\rm end}}{a(t)}\right)\,.$
(2.14)

The time-dependence of the quasi-periodic function $\mathcal{P}(t)$ can be
obtained by exploiting the underdamped nature of the oscillations.
Approximating the envelope as constant during one oscillation, the equation of
motion (2.7) can be written as

$\dot{\mathcal{P}}^{2}\;=\;\frac{1}{6}m_{\phi}^{2}\left(1-\mathcal{P}^{4}\right)\,.$
(2.15)

where the effective mass of the inflaton is defined as

$m_{\phi}^{2}(t)\;\equiv\;V_{\phi\phi}(\phi_{0}(t))\;=\;12\lambda\phi_{0}^{2}(t)\;\simeq\;12\lambda\phi_{\rm
end}^{2}\left(\frac{a_{\rm end}}{a(t)}\right)^{2}\;\equiv\;m_{\rm
end}^{2}\left(\frac{a_{\rm end}}{a(t)}\right)^{2}\,,$ (2.16)

with $V_{\phi\phi}=\partial^{2}V/\partial\phi^{2}$. Equation (2.15) can be
solved in terms of Jacobi elliptic functions, $\mathcal{P}(t)\;=\;{\rm
sn}\,(m_{\phi}t/\sqrt{6},-1)$ which can be expanded in terms of Fourier
coefficients as

$\mathcal{P}(t)=\sum_{n=-\infty}^{\infty}\mathcal{P}_{n}e^{-in\omega_{\phi}t}\,,$
(2.17)

where the sum is taken over the harmonic modes of the $\phi$ oscillation with
the fundamental angular frequency given by

$\omega_{\phi}\;=\;m_{\phi}\sqrt{\frac{2}{3}\pi}\,\frac{\Gamma(3/4)}{\Gamma(1/4)}\,.$
(2.18)

In order to account for expansion, we can indirectly exploit the conformal
invariance of the model. Since $a\;\propto\;t^{1/2}$, conformal time is
$\tau\equiv\int a^{-1}\,\mathop{}\\!\mathrm{d}t\propto a$, and
$m_{\phi}t\;=\;m_{\rm end}(\tau-\tau_{\rm end})$. We can therefore
characterize the stage of coherent oscillations of the inflaton by means of
the approximate solution

$\phi(\tau)\;\simeq\;\phi_{\rm end}\left(\frac{a_{\rm
end}}{a(\tau)}\right){\rm sn}\,\left(\frac{m_{\rm
end}}{\sqrt{6}}(\tau-\tau_{\rm end}),-1\right)\,.$ (2.19)

###  2.2 Parametric resonance

In the Introduction we have made the case for the study of the growth of the
inflaton fluctuations during reheating, driven by the self-interaction. We now
proceed to study the resonant growth of the inflaton momentum modes in the
linear approximation. For these purposes, we denote the non-homogeneous
inflaton perturbation as $\delta\phi(t,\boldsymbol{x})$. Variation of the
action (2.1), disregarding couplings of $\phi$ to other degrees of freedom,
leads to the following dynamical equation,333In this work we ignore the
fluctuations of the metric, which can potentially be important during the
preheating stage [53, 54, 55, 56, 57, 58].

$\ddot{\delta\phi}+3H\dot{\delta\phi}-\frac{\nabla^{2}\delta\phi}{a^{2}}+12\lambda\phi(t)^{2}\,\delta\phi\;=\;0\,.$
(2.20)

Introducing the canonically normalized fluctuation $X\equiv a\,\delta\phi$,
and switching to conformal time, we can write the canonically quantized field
as

$X(\tau,\boldsymbol{x})=\int\frac{\mathop{}\\!\mathrm{d}^{3}\boldsymbol{k}}{(2\pi)^{3/2}}\,e^{-i\boldsymbol{k}\cdot\boldsymbol{x}}\left[X_{k}(\tau)\hat{a}_{\boldsymbol{k}}+X_{k}^{*}(\tau)\hat{a}^{\dagger}_{-\boldsymbol{k}}\right]\,,$
(2.21)

where $\boldsymbol{k}$ denotes the comoving momentum, and
$\hat{a}_{\boldsymbol{k}}$ and $\hat{a}^{\dagger}_{\boldsymbol{k}}$ are the
annihilation and creation operators, respectively, satisfying the canonical
commutation relations
$[\hat{a}_{\boldsymbol{k}},\hat{a}^{\dagger}_{\boldsymbol{k}^{\prime}}]=\delta(\boldsymbol{k}-\boldsymbol{k}^{\prime})$,
$[\hat{a}_{\boldsymbol{k}},\hat{a}_{\boldsymbol{k}^{\prime}}]=[\hat{a}^{\dagger}_{\boldsymbol{k}},\hat{a}^{\dagger}_{\boldsymbol{k}^{\prime}}]=0$.
We ensure that the corresponding canonical commutation relations between the
field, $X_{k}$, and its momentum conjugate, $X_{k}^{\prime}$, are fulfilled by
imposing the Wronskian constraint
$X_{k}X^{*\prime}_{k}-X_{k}^{*}X_{k}^{\prime}\;=\;i$. In this decomposition,
the equation of motion (2.20) reduces to

$X_{\boldsymbol{k}}^{\prime\prime}+\left(k^{2}-\frac{a^{\prime\prime}}{a}+12\lambda\phi^{2}a^{2}\right)X_{\boldsymbol{k}}\;=\;0\,,$
(2.22)

where ′ denotes the derivative with respect to conformal time. Examination of
the factor in parenthesis shows that the second term decreases over time; for
it we have

$\frac{a^{\prime\prime}}{a}\;=\;\frac{a^{2}}{6M_{P}^{2}}(4V-\dot{\phi}^{2})\;\sim\;\frac{\lambda\phi^{4}a^{2}}{M_{P}^{2}}\;\rightarrow\;0\,.$
(2.23)

On the other hand, the third term inside the parenthesis of Eq. (2.22) does
not decrease with time (cf. (2.14)), and its time dependence is only modulated
by $\mathcal{P}$. Introducing the dimensionless variable $z\equiv m_{\rm
end}(\tau-\tau_{\rm end})$, we then find that soon after the onset of
reheating, Eq. (2.22) can be approximated as

Figure 1: Floquet chart for Eq. (2.24). Due to the conformal nature of the
model (2.7), a mode that begins inside a resonance band (shown in gray) stays
indefinitely in the band, regardless of the expansion of the universe.

$\dfrac{\mathop{}\\!\mathrm{d}^{2}X_{\boldsymbol{k}}}{\mathop{}\\!\mathrm{d}z^{2}}+\left[\left(\frac{k}{m_{\rm
end}}\right)^{2}+{\rm
sn}^{2}\left(\frac{z}{\sqrt{6}},-1\right)\right]X_{\boldsymbol{k}}\;=\;0\,.$
(2.24)

Eq. (2.24) has the form of Hill’s equation (known as Mathieu’s equation in the
harmonic case), and presents parametric resonance. Floquet’s theorem
guarantees that the solutions to this equation have the form [59]

$X_{\boldsymbol{k}}(z)\;=\;e^{\mu_{k}z}g_{1}(z)+e^{-\mu_{k}z}g_{2}(z)\,,$
(2.25)

where $g_{1}(z)$ and $g_{2}(z)$ are periodic functions, and $\mu_{k}$ is a
complex number, called the Floquet exponent. Exponentially growing solutions
are therefore found if ${\rm Re}\,\mu_{k}\;>\;0$. In order to determine the
Floquet exponents, we follow the eigenvalue method described in detail in
[22]. The result of this numerical computation is shown in Fig. 1. The
simplicity of the Floquet chart is noteworthy, as only two main resonance
bands appear. One at small momentum values, $k/m_{\rm end}\lesssim 1.8\times
10^{-4}$, and a second one, with a narrow width and Floquet exponent more than
two orders of magnitude larger, at $0.71\lesssim k/m_{\rm end}\lesssim 0.76$.
Due to the conformal nature of the quartic potential, the Floquet exponents
are insensitive to the expansion of the universe, and a mode within a
resonance band will stay in the band while the oscillation of the inflaton
lasts. Therefore, modes with $k/m_{\rm end}\simeq 0.7$ will rapidly grow at
the onset of reheating, and become non-linear after a few $e$-folds, even for
a small coupling constant $\lambda$.

###  2.3 Backreaction and fragmentation

At linear order, the growth of the momentum modes of $X$ can be tracked by
solving Eq. (2.22) with the appropriate initial conditions, which we take as
the positive-frequency Bunch-Davies vacuum

$X_{k}(\tau_{0})\;=\;\frac{1}{\sqrt{2\omega_{k}}}\,,\quad
X_{k}^{\prime}(\tau_{0})\;=\;-\frac{i\omega_{k}}{\sqrt{2\omega_{k}}}\,,$
(2.26)

where

$\omega_{k}^{2}\;=\;k^{2}-\frac{a^{\prime\prime}}{a}+12\lambda\phi^{2}a^{2}\,.$
(2.27)

Our main quantity of interest will be the phase space distribution (PSD) of
the fluctuations, which coincides with the comoving occupation number of
$\phi$. Its UV-finite form is given by [18]

$f_{\delta\phi}(k,t)\;=\;n_{k}\;=\;\frac{1}{2\omega_{k}}\left|\omega_{k}X_{k}-iX^{\prime}_{k}\right|^{2}\,.$
(2.28)

Two other important quantities can be obtained from this PSD, the number
density and the energy density of the fluctuations. Their UV-regular forms are
computed as [18, 23]

$\displaystyle n_{\delta\phi}\;$
$\displaystyle=\;\frac{1}{(2\pi)^{2}a^{3}}\int\mathop{}\\!\mathrm{d}^{3}\boldsymbol{k}\,n_{k}\,,$
(2.29) $\displaystyle\rho_{\delta\phi}\;$
$\displaystyle=\;\frac{1}{(2\pi)^{2}a^{4}}\int\mathop{}\\!\mathrm{d}^{3}\boldsymbol{k}\,\omega_{k}n_{k}\,.$
(2.30)

The resulting form of the PSD in the linear approximation is shown as the
dashed lines in the left panel of Fig. 2. As expected, the PSD is peaked
around $k/m_{\rm end}=0.73$, in full agreement with the results of the Floquet
analysis of the previous section. The height of this peak increases with time,
as the resonant growth accumulates during the oscillation of $\phi$. We only
show three snapshots of the PSD for $a\leq 130$, since the system rapidly
evolves toward the non-linear regime, making the analysis based on the
solution of (2.22) unsuitable.

Figure 2: The particle phase space distribution (PSD) for the inflaton
fluctuations for a selection of scale factors, coded by color (left), and the
corresponding scale factor dependence of the comoving number density (right).
The lattice results are shown in both panels as solid curves. In the left
panel, only for the first three values of $a$ we also show the PSD as computed
from the linear spectral analysis, as dashed lines. The main Floquet band is
shown in light gray.

The impact of the rapid growth of the resonant modes on the energy density of
the inflaton fluctuations can also be appreciated in Fig. 3. We follow the
scale factor dependence of $\rho_{\delta\phi}$ (the orange curve) by means of
(2.30) only up to $a/a_{\rm end}=10^{2}$. Within this range, the growth of the
energy density in fluctuations is clear, ramping up at $a/a_{\rm end}\sim 30$,
and growing exponentially fast afterwards. The transition from the linear to
the non-linear regime is clear.

From this discussion it follows that, despite the smallness of $\lambda$, the
growth of the momentum modes in the main resonance band is strong enough to
eventually be able to draw an $\mathcal{O}(1)$ fraction of the energy density
of the zero-mode of $\phi$. The inflaton “fluctuation” is no longer such, and
the system enters the non-linear regime. Mode-mode couplings become important,
redistributing the energy of the resonant mode into other modes
(rescattering), leading to large configuration space gradients, and
effectively fragmenting the once homogeneous condensate. This phenomenon is
known as backreaction, and it does not allow for a straightforward spectral
analysis. Instead, the full non-linear configuration-space equation of motion
for the full operator $\phi(t,\boldsymbol{x})$ must be solved, which is also a
non-trivial task. To get around some of the complications, it is argued that,
since in the non-linear regime occupation numbers are large, $n_{k}\gg 1$, the
dynamics of the system can be adequately tracked by approximating $\phi$ as a
classical field. The full non-linear PDE

$\ddot{\phi}+3H\dot{\phi}-\frac{\nabla^{2}\phi}{a^{2}}+V_{\phi}\;=\;0\,,$
(2.31)

is then solved over a configuration-space lattice. The energy density of the
inflaton is computed from the spatial average of the energy-momentum tensor of
$\phi$, which we denote with an over-bar,

$\rho_{\phi}\;=\;\mkern
1.5mu\overline{\mkern-1.5mu\frac{1}{2}\dot{\phi}^{2}+\frac{1}{2a^{2}}(\nabla\phi)^{2}+V(\phi)\mkern-1.5mu}\mkern
1.5mu\,.$ (2.32)

Spectral data in turn is obtained upon Fourier transformation of
configuration-space quantities. The quartic inflaton system has been
extensively studied in the past by means of lattice methods, see e.g. [21, 30,
31, 32, 33, 34, 35, 37]. For our analysis we use the publicly available code
$\mathcal{C}\text{osmo}\mathcal{L}\text{attice}$ [60, 61].

Figure 3: Top: Pre- and post-fragmentation evolution of the total inflaton
energy density (blue), the energy density of the inhomogeneous component of
$\phi$ (orange) and the energy density of the homogeneous component (green).
For $\rho_{\delta\phi}$, the linear approximation (2.30) is depicted up to
$a/a_{\rm end}\leq 10^{2}$. At later times the lattice result
$\rho_{\delta\phi}=\rho_{\phi}-\mkern
1.5mu\overline{\mkern-1.5mu\rho_{\phi}\mkern-1.5mu}\mkern 1.5mu$ is shown.
Bottom: Oscillation-averaged inflaton equation of state parameter.

Our lattice results are summarized in Figs. 2 and 3. The left panel of Fig. 2
shows, as continuous curves, the form of the PSD of the inflaton as computed
by the lattice code. At early times it shares the features of the linear
approximation, showing a main growing peak, located at the momentum band
predicted by the Floquet analysis. For the earliest times, the lattice PSD
coincides with the linear prediction. Nevertheless, for $a/a_{\rm end}\gtrsim
10^{2}$, the effect of rescattering becomes evident. The growth of the
resonant mode is traded for the population of the PSD for lower momentum
modes, and larger momentum adjacent modes. As time passes, the redistribution
of energy becomes more efficient, and a UV tail in the distribution appears.
At first this tail is populated more efficiently at some discrete values of
$k$. Notably, these wavenumbers can be estimated by means of the Boltzmann
approximation, as we show in Section 2.4. In any case, these features are
quickly erased, and a smooth limiting distribution emerges, shown in red. The
right panel of Fig. 2 conveys this asymptotic behavior of the PSD, in terms of
the comoving number density of $\phi$. At first $n_{\delta\phi}$ grows
rapidly, driven by the parametric resonance. However, for $a/a_{\rm
end}\gtrsim 300$, this growth stops, the comoving number density freezes, and
the effect of rescatterings is merely the kinetic redistribution of energy
among modes.

Fig. 3 contains the lattice results for the total energy density of $\phi$,
the density of the condensate, and that of free particles, for $a/a_{\rm
end}>10^{2}$. For definiteness, and in order to connect with the results of
the spectral analysis, we define the condensate component of the energy
density as follows [23],

$\overline{\rho_{\phi}}\;=\;\frac{1}{2}\mkern
1.5mu\overline{\mkern-1.5mu\dot{\phi}\mkern-1.5mu}\mkern 1.5mu^{2}+V(\mkern
1.5mu\overline{\mkern-1.5mu\phi\mkern-1.5mu}\mkern 1.5mu)\,,$ (2.33)

that is, the energy density of the spatially averaged inflaton field. In turn,
we take $\rho_{\delta\phi}=\rho_{\phi}-\mkern
1.5mu\overline{\mkern-1.5mu\rho_{\phi}\mkern-1.5mu}\mkern 1.5mu$. This
definition of the fluctuation energy density matches excellently with the
spectral result, and is shown as the orange curve in Fig. 3. Its value rapidly
grows until $a/a_{\rm end}\simeq 180$, point at which the backreaction in the
oscillating condensate is noticeable. The total energy density of $\phi$
redshifts as radiation even after rescattering becomes important, but the
energy density of the zero mode is now efficiently transferred to
fluctuations. Interestingly, despite the initial rapid growth in
$\rho_{\delta\phi}$, the dissipation of $\mkern
1.5mu\overline{\mkern-1.5mu\rho_{\phi}\mkern-1.5mu}\mkern 1.5mu$ after
entering backreaction is gradual, $\mkern
1.5mu\overline{\mkern-1.5mu\rho_{\phi}\mkern-1.5mu}\mkern 1.5mu\propto
a^{-5.3}$, unlike the exponential decrease of preheating scenarios in
quadratic potentials [23]. The survival of this coherent component will be an
important ingredient in our exploration of the decay of $\phi$ of Section 3.

###  2.4 The Boltzmann approximation

In the previous sections we have shown that a combination of spectral and
lattice methods is necessary to correctly determine the distribution of the
inflaton fluctuations. We now take the opportunity to show that, in the linear
regime, it is possible to extract non-trivial information for the PSD by means
of the integration of the Boltzmann equation for the inflaton fluctuations.
Inflaton quanta $\delta\phi$ sourced from the inflaton background $\phi$ are
induced by the interaction Lagrangian

$\mathcal{L}_{I}=6\lambda\phi^{2}\delta\phi^{2}\,.$ (2.34)

Following [9, 12, 24], the Boltzmann equation in the presence of anharmonic
oscillations of the inflaton background takes the form

$\begin{split}\frac{\partial f_{\delta\phi}}{\partial
t}-H|\boldsymbol{P}|\frac{\partial
f_{\delta\phi}}{\partial|\boldsymbol{P}|}\,=\,\frac{1}{P^{0}}&\sum_{n=1}^{\infty}\int\frac{\mathop{}\\!\mathrm{d}^{3}\boldsymbol{K}_{n}}{(2\pi)^{3}n_{\phi}}\frac{\mathop{}\\!\mathrm{d}^{3}\boldsymbol{P}^{\prime}}{(2\pi)^{3}2P^{\prime
0}}(2\pi)^{4}\delta^{(4)}\left(K_{n}-P-P^{\prime}\right)\left|\mkern
1.5mu\overline{\mkern-1.5mu\mathcal{M}_{n}\mkern-1.5mu}\mkern
1.5mu\right|^{2}\\\
&\times\left[f_{\phi}\left(K_{n}\right)\left(1+f_{\delta\phi}(P)\right)\left(1+f_{\delta\phi}\left(P^{\prime}\right)\right)-f_{\delta\phi}(P)f_{\delta\phi}\left(P^{\prime}\right)\left(1+f_{\phi}\left(K_{n}\right)\right)\right].\end{split}$
(2.35)

Here we denote the physical four-momenta of the fluctuations by $P$ and
$P^{\prime}$. $K_{n}=\left(E_{n},\bf{0}\right)$ is the physical four-momentum
of the inflaton condensate, where $E_{n}=n\,\omega_{\phi}$ denotes the energy
of the $n^{\text{th}}$ oscillating mode. $\mathcal{M}_{n}$ represents the
transition amplitude corresponding to the production of a pair of inflaton
quanta $|f\rangle=|\delta\phi\delta\phi\rangle$ from the $n^{\text{th}}$
Fourier mode of the coherently oscillating inflaton background field.

$\big{|}\big{\langle}f\big{|}i\int\mathop{}\\!\mathrm{d}^{4}x\,\mathcal{L}_{I}\big{|}0\big{\rangle}\big{|}^{2}=\operatorname{Vol}_{4}\sum_{n=-\infty}^{\infty}\left|\overline{\mathcal{M}_{n}}\right|^{2}(2\pi)^{4}\delta^{(4)}\left(K_{n}-P-P^{\prime}\right).$
(2.36)

where Vol4 denote the space-time volume. From Eq. (2.34), the matrix elements
$\mathcal{M}_{n}$ are found to be

$\mathcal{M}_{n}\;=\;12\lambda\phi_{0}^{2}\,\hat{\mathcal{P}}_{n}\qquad\Rightarrow\qquad\left|\overline{\mathcal{M}_{n}}\right|^{2}=72\lambda^{2}\phi_{0}^{4}\,|\hat{\mathcal{P}}_{n}|^{2}\,.$
(2.37)

The mean squared amplitude $|\overline{\mathcal{M}_{n}}|^{2}$ accounts for a
factor of 2 for identical particles in the final state. The
$\hat{\mathcal{P}}_{n}$ are the Fourier coefficients of the harmonic
decomposition over one oscillation of the square of the quasi-periodic
function $\mathcal{P}(t)$ defined in (2.9),

$\mathcal{P}^{2}(t)\;=\;\sum_{n=-\infty}^{\infty}\hat{\mathcal{P}}_{n}e^{-in\omega_{\phi}t}\,.$
(2.38)

The spatially homogeneous background inflaton condensate distribution can be
expressed as the zero-mode

$f_{\phi}(\boldsymbol{K},t)=(2\pi)^{3}n_{\phi}(t)\delta^{(3)}(\boldsymbol{K})\,.$
(2.39)

Neglecting the backreaction of the inflaton quanta onto the condensate, the
Boltzmann equation takes the following form

$\displaystyle\frac{\partial f_{\delta\phi}}{\partial t}$
$\displaystyle\;-\;H|\boldsymbol{P}|\frac{\partial
f_{\delta\phi}}{\partial|\boldsymbol{P}|}$
$\displaystyle\simeq\;\sum_{n=1}^{\infty}\frac{1}{P^{0}}\int\frac{\mathop{}\\!\mathrm{d}^{3}\boldsymbol{K}_{n}}{(2\pi)^{3}n_{\phi}}\frac{\mathop{}\\!\mathrm{d}^{3}\boldsymbol{P}^{\prime}}{(2\pi)^{3}2P^{\prime
0}}(2\pi)^{4}\delta^{(4)}\left(K_{n}-P-P^{\prime}\right)\left|\overline{\mathcal{M}_{n}}\right|^{2}f_{\phi}\left(K_{n}\right)\left(1+f_{\delta\phi}(P)+f_{\delta\phi}\left(P^{\prime}\right)\right)$
$\displaystyle\simeq\;72\pi\lambda^{2}\phi_{0}^{4}\sum_{n=1}^{\infty}\int\frac{\mathop{}\\!\mathrm{d}^{3}\boldsymbol{P}^{\prime}}{P^{0}P^{\prime
0}}\delta\left(n\omega_{\phi}-P^{0}-P^{\prime
0}\right)\delta^{(3)}\left(\boldsymbol{P}+\boldsymbol{P}^{\prime}\right)\,|\hat{\mathcal{P}}_{n}|^{2}\left(1+f_{\delta\phi}(P)+f_{\delta\phi}\left(P^{\prime}\right)\right)$
$\displaystyle=\;144\pi\lambda^{2}\phi_{0}^{4}\sum_{n=1}^{\infty}\frac{|\hat{\mathcal{P}}_{n}|^{2}}{n^{2}\omega_{\phi}^{2}\beta_{n}}\delta\left(|\boldsymbol{P}|-\frac{1}{2}n\omega_{\phi}\beta_{n}\right)\Big{(}1+2f_{\delta\phi}(|\boldsymbol{P}|)\Big{)}\,,$
(2.40)

where

$\beta_{n}\;=\;\sqrt{1-\frac{4m_{\phi}^{2}}{n^{2}\omega_{\phi}^{2}}}\;=\;\sqrt{1-\left(\frac{2}{nc}\right)^{2}},\quad\quad
c\;\equiv\;\sqrt{\frac{2\pi}{3}}\frac{\Gamma(3/4)}{\Gamma(1/4)}\,,$ (2.41)

denotes the (time-independent) kinematic suppression factor. The quantum
mechanical Bose enhancement factor, which appears in the right side of Eq.
(2.40), can be scaled out by defining the “classical” distribution
$f_{\chi}^{c}$ [62, 63, 64]:

$f_{\delta\phi}(|\boldsymbol{P}|,t)\;\equiv\;\frac{1}{2}\Big{[}\exp\big{(}2f_{\delta\phi}^{c}(|\boldsymbol{P}|,t)\big{)}-1\Big{]}\,.$
(2.42)

The Boltzmann equation reduces then to

$\displaystyle\frac{\partial f^{c}_{\delta\phi}}{\partial
t}-H|\boldsymbol{P}|\frac{\partial
f^{c}_{\delta\phi}}{\partial|\boldsymbol{P}|}\;=\;144\pi\lambda^{2}\phi_{0}^{4}\sum_{n=1}^{\infty}\frac{|\hat{\mathcal{P}}_{n}|^{2}}{n^{2}\omega_{\phi}^{2}\beta_{n}}\delta\left(|\boldsymbol{P}|-\frac{1}{2}n\omega_{\phi}\beta_{n}\right)\,.$
(2.43)

As the right-hand side of Eq. (2.43) is $f^{c}_{\delta\phi}$-independent, we
can immediately integrate this equation to get [65, 66]

$\displaystyle f^{c}_{\delta\phi}(|\boldsymbol{P}|,t)$
$\displaystyle\simeq\;144\pi\lambda^{2}\sum_{n=1}^{\infty}\frac{|\hat{\mathcal{P}}_{n}|^{2}}{n^{2}c^{2}\beta_{n}}\int_{t_{\mathrm{end}}}^{t}\mathop{}\\!\mathrm{d}t^{\prime}\frac{\phi_{0}^{4}\left(t^{\prime}\right)}{m_{\phi}^{2}\left(t^{\prime}\right)}\delta\left(\frac{a(t)}{a\left(t^{\prime}\right)}|\boldsymbol{P}|-\frac{1}{2}nc\beta_{n}m_{\phi}(t^{\prime})\right)$
(2.44)
$\displaystyle\simeq\;\frac{\pi}{c^{2}}\left(\frac{m_{\mathrm{end}}}{H_{\mathrm{end}}}\right)\left(\frac{a(t)}{a_{\mathrm{end}}}-1\right)\sum_{n=1}^{\infty}\frac{|\hat{\mathcal{P}}_{n}|^{2}}{n^{2}\beta_{n}}\delta\left(q-\frac{1}{2}nc\beta_{n}\right)\,,$

with $q\equiv a(t)|\boldsymbol{P}|/(a_{\text{end}}m_{\text{end}})$ the re-
scaled comoving momentum referred to the end of inflation, denoted as
$k/m_{\rm end}$ in Fig. 2.

Peak structure. The distribution function corresponds to a collection of Dirac
delta functions, located at different momenta $(1/2)nc\beta_{n}$. The
$\beta_{n}$ coefficient is a decreasing function of $n$ whose first (real)
non-vanishing value is achieved for $n=5$. By denoting the location of the
$i^{\text{th}}$ peak as $\hat{q}_{i}$, the distribution can be equivalently
expressed as as sum of contributing modes via

$f^{c}_{\delta\phi}(q,t)\,=\,\sum_{i=1}^{\infty}\mathcal{F}^{c}_{i}(t)\,\delta\left(q-\hat{q}_{i}\right)\,,$
(2.45)

with time-dependent coefficients $\mathcal{F}^{c}_{i}(t)$ and peak locations
determined by

$\hat{q}_{i}\,=\,\sqrt{\frac{2\pi}{3}}\,\left(\dfrac{i+4}{2}\right)\dfrac{\Gamma(3/4)}{\Gamma(1/4)}\sqrt{1-\left(\frac{6}{\pi(i+4)^{2}}\right)\left(\dfrac{\Gamma(1/4)}{\Gamma(3/4)}\right)^{2}}\,,$
(2.46)

where $i=n-4$. This ratio is remarkably independent of any model parameter. It
is straightforward to check that the first peak $\hat{q}_{1}$ is predicted to
be at $\hat{q}_{1}\simeq 0.7$, in excellent agreement with the Floquet
analysis and the linear and lattice results. The relative location of peaks at
larger momenta can be straightforwardly inferred from Eq. (2.46) and
corresponds to the following ratio

$\dfrac{\hat{q}_{i}}{\hat{q}_{j}}\,=\,\dfrac{(i+4)\beta_{i+4}}{(j+4)\beta_{j+4}}\,.$
(2.47)

Limitation of the Boltzmann approach. The previous computation resulted in a
distribution peaked at discrete values of momenta, with an infinitesimally
narrow width and a large amplitude peak characteristic of a Dirac delta
function. In reality however, as we have seen in the previous sections, the
distribution is spread around the peaks. Immediately after the end of
inflation, from our linear approach we found a distribution function at the
first peak $f_{\delta\phi}(\hat{q}_{1},t_{\text{end}})\sim 0.2$. Moreover, the
mode amplitudes are expected to grow continuously
$\mathcal{F}^{c}_{i}(t)\propto a(t)/a_{\text{end}}$ during reheating, becoming
rapidly even larger. For such large amplitudes the Boltzmann approach breaks
down almost immediately [63, 24].

The amplitude of the peaks, therefore, cannot be accurately estimated from the
Boltzmann approach starting from the first instants after the end of
inflation. However, we find that this approach gives a qualitatively good
description of the location of the PSD peaks appearing at the onset of
backreaction, and notably, for the peaks in the induced gravitational wave
spectrum. We postpone the detailed comparison to Section 4.

##  3 Reheating in a quartic potential

We now turn to the study of the decay of the inflaton into light degrees of
freedom, necessary to complete reheating and populate the universe with a
relativistic plasma in thermal equilibrium. As mentioned in the Introduction,
our analysis will be based on the assumption that fields coupled to the
inflaton are not strongly sourced via parametric resonance. To satisfy this
assumption, we will implicitly assume that the main decay products of $\phi$
are spin 1/2 fermions, denoted as $\psi$. Moreover, we will assume for
simplicity that the decay of the inflaton occurs at tree-level, and produces
two particles in the final state. Thus, generically, we consider

$\mathcal{L}_{\rm
int}\;=\;g\phi\bar{\psi}\psi+ig^{\prime}\phi\bar{\psi}\gamma_{5}\psi\,.$ (3.1)

It is worth noting that the Yukawa term in the previous expression can lead to
fermion preheating [67, 68, 69, 70, 71], for which the production of particles
is resonantly suppressed, a consequence of Fermi-Dirac statistics.
Perturbatively, the oscillating $\phi$ induces an effective mass for $\psi$
than can kinematically block the decay into $\psi$ very efficiently [12]. The
exploration of these interesting effects is not the main purpose of the
present work, and will be postponed for a future study. We will therefore
disregard these mechanisms, a possibility that arises for $g^{\prime}\gg g$.

###  3.1 Production rates

Schematically, the Boltzmann equation that must be integrated to track the
phase space distribution of the decay products can be written as

$\frac{\partial f_{\psi}}{\partial t}-H|\boldsymbol{P}|\frac{\partial
f_{\psi}}{\partial|\boldsymbol{P}|}\;=\;\mathcal{C}[f_{\psi}]\,,$ (3.2)

with $\mathcal{C}$ denoting the collision term. For generality (and future
applications) let us for now assume a non-vanishing, and potentially time-
dependent mass $m_{\psi}(t)$ for the decay products. Applying the operator
$\int\mathop{}\\!\mathrm{d}^{3}\boldsymbol{P}\,P^{0}/(2\pi)^{3}\,$ on both
sides of (3.2), we obtain the following general form for the continuity
equation for the energy density $\rho_{\psi}$,

$\dot{\rho}_{\psi}+4H\rho_{\psi}-\int\frac{\mathop{}\\!\mathrm{d}^{3}\boldsymbol{P}}{(2\pi)^{3}P^{0}}\,f_{\psi}(P)\,m_{\psi}\left(\dot{m}_{\psi}+Hm_{\psi}\right)\;=\;\int\frac{\mathop{}\\!\mathrm{d}^{3}\boldsymbol{P}}{(2\pi)^{3}}\,P^{0}\mathcal{C}[f_{\psi}]\;\equiv\;R(t)\,.$
(3.3)

Here we introduce $R(t)$, the radiation-energy production rate per unit of
volume.

In the following we assume that the states $\psi$ eventually thermalize among
themselves and the rest of the Standard Model states. We therefore identify
$\rho_{\psi}=\rho_{R}$, where the later denotes the energy density of the
relativistic radiation plasma. The evaluation of the production rate in Eq.
(3.3) will be split in two. We estimate the contributions to the production
rate from the oscillating inflaton condensate $\phi$ and from the fragmentated
inflaton quanta $\delta\phi$

$R(t)\,\equiv\,R_{\phi}(t)+R_{\delta\phi}(t)\,.$ (3.4)

####  3.1.1 Decay of the coherent oscillations

Let us first briefly discuss the decay of the inflaton at early times, when it
is composed mainly of the coherent, oscillating zero-mode. Under the
assumption that quantum statistics do not play a role in the dissipation
process, the Boltzmann equation which describes the decay of the inflaton can
be written as [10, 8, 12]

$\dot{\rho}_{\phi}+3H(1+w_{\phi})\rho_{\phi}\;=\;-R_{\phi}(t)\,.$ (3.5)

The rate in this equation is evaluated as

$R_{\phi}(t)\,=\,(1+w_{\phi})\Gamma_{\phi}\rho_{\phi}\,,$ (3.6)

where the decay rate $\Gamma_{\phi}$ for a generic process $\phi\rightarrow
A+B$ is

$\Gamma_{\phi}\;=\;\frac{1}{8\pi(1+w_{\phi})\rho_{\phi}}\sum_{n=1}^{\infty}|\mathcal{M}_{n}|^{2}E_{n}\sqrt{\left(1-\frac{(m_{A}+m_{B})^{2}}{E_{n}^{2}}\right)\left(1-\frac{(m_{A}-m_{B})^{2}}{E_{n}^{2}}\right)}\,.$
(3.7)

Here the sum is taken over the harmonic modes of the function $\mathcal{P}(t)$
over one oscillation, with associated energy $E_{n}=n\omega_{\phi}$ (see Eq.
(2.17)). $\mathcal{M}_{n}$ denotes the transition amplitude in one oscillation
for each mode from the coherent state to the two-particle state that can be
defined analogously to Eq. (2.36), with $|f\rangle=|A,B\rangle$ and
$\mathcal{L}_{I}=\mathcal{L}_{\text{int}}$, cf. (2.1). In the case of an
inflaton oscillating about a quartic potential, assuming a fermionic decay for
which the masses can be disregarded, the decay rate $\Gamma_{\phi}$ can be
evaluated explicitly, taking the form [12]

$\Gamma_{\phi\rightarrow\bar{\psi}\psi}(t)\;=\;\alpha^{2}\,\frac{y^{2}}{8\pi}m_{\phi}(t)\,,$
(3.8)

where $y=\\{g,g^{\prime}\\}$, $m_{\phi}$ is the envelope inflaton mass (2.16),
and $\alpha\simeq 0.71$ is an efficiency factor that codifies the
anharmonicity of the $\phi$ oscillation [12]. The continuity-Friedmann system
that completes Eq. (3.5) corresponds to

$\displaystyle\dot{\rho}_{R}+4H\rho_{R}\;$
$\displaystyle=\;(1+w_{\phi})\,\Gamma_{\phi}\,\rho_{\phi}\,,$ (3.9)
$\displaystyle\rho_{\phi}+\rho_{R}\;$ $\displaystyle=\;3M_{P}^{2}H^{2}\,,$
(3.10)

with $w_{\phi}=1/3$.

####  3.1.2 Decay of the fragmented inflaton

In Section 2.3 we discuss how the inflaton is fragmented after it enters the
backreaction regime. When this occurs, the production of particles can no
longer be characterized by (3.9) with decay rate (3.8). Instead, one must now
determine the form of the continuity equation for $\rho_{R}$ from the
microscopic Boltzmann equation, assuming a population of inflaton particles
with PSD $f_{\delta\phi}(K)$. Disregarding inverse decays and Pauli blocking /
Bose enhancement, the collision term in the Boltzmann equation (3.2) has the
form

$\displaystyle\mathcal{C}_{\delta\phi}[f_{\psi}]\;=\;\frac{1}{P^{0}}\int$
$\displaystyle\frac{\mathop{}\\!\mathrm{d}^{3}\boldsymbol{K}}{(2\pi)^{3}2K^{0}}\frac{\mathop{}\\!\mathrm{d}^{3}\boldsymbol{P}^{\prime}}{(2\pi)^{3}2P^{0\prime}}(2\pi)^{4}\delta^{(4)}(K-P-P^{\prime})|\mathcal{M}_{\delta\phi\rightarrow\bar{\psi}\psi}|^{2}f_{\delta\phi}(K)\,.$
(3.11)

The particle production rate can then be evaluated in a straightforward way,
and yields

$R_{\delta\phi}(t)\;=\;\int\frac{\mathop{}\\!\mathrm{d}^{3}\boldsymbol{P}}{(2\pi)^{3}}\,P^{0}\,\mathcal{C}_{\delta\phi}[f_{\psi}]\;=\;\Gamma_{\delta\phi}m_{\phi}n_{\delta\phi}\,,$
(3.12)

where

$\Gamma_{\delta\phi}\;=\;\frac{|\mathcal{M}_{\delta\phi\rightarrow\bar{\psi}\psi}|^{2}}{16\pi
m_{\phi}}\sqrt{1-\frac{4m_{\psi}^{2}}{m_{\phi}^{2}}}\,,$ (3.13)

denotes the decay rate of free inflaton particles $\delta\phi$. In the case
corresponding to (3.1),
$\Gamma_{\delta\phi\rightarrow\bar{\psi}\psi}=y^{2}m_{\phi}/8\pi$ for
negligible masses for the outgoing states.

####  3.1.3 Total rate

Summarizing our previous findings, the total rate (3.4) accounting for
production from the inflaton condensate and particle rates reads

$\displaystyle R(t)\;$
$\displaystyle=\;\frac{4}{3}\Gamma_{\phi}\overline{\rho_{\phi}}+\Gamma_{\delta\phi}m_{\phi}n_{\delta\phi}$
(3.14)
$\displaystyle\;=\;\frac{y^{2}}{8\pi}m_{\phi}\left(\frac{4}{3}\alpha^{2}\overline{\rho_{\phi}}+m_{\phi}n_{\delta\phi}\right)\,.$
(3.15)

The second line corresponds to the two-body fermionic decay channel, in the
limit where $m_{\psi}\ll m_{\phi}$.

Figure 4: Scale factor dependence of the effective inflaton mass (left), and
of the effective Boltzmann rate for the relativistic decay products $R(t)$
(right). The total rate is the sum of the condensate contribution (black) and
the particle contribution (orange).

Fig. 8 shows the scale factor dependence of the effective mass of the
inflaton, inherited from its coherent oscillation (left), and that of the rate
$R(t)$ decomposed into the condensate ($R_{\phi}(t)$) and particle
($R_{\delta\phi}(t)$) contributions (right). In both panels we have multiplied
the corresponding variable by its redshift factor prior to $\phi$
fragmentation. For the effective mass, we have $m_{\phi}(a/a_{\rm
end})\simeq{\rm const.}$ quickly after the beginning of reheating, until
$a/a_{\rm end}\simeq 160$, at the onset of strong rescattering. Due to the
efficient conversion of the inflaton zero-mode into finite momentum particles,
$m_{\phi}$ is further reduced, first in a non-monotonic manner while strong
backreaction takes place, and monotonically for $a\gtrsim 300\,a_{\rm end}$,
with $m_{\phi}\propto(a/a_{\rm end})^{1.32}$, as determined from the lattice
data.

In the right panel of Fig. 8 we show the condensate (black) and particle
(orange) contributions to the total rate $R(t)$, rescaled by its scale factor
dependence and by the amplitude parameter $y^{2}/8\pi$. Similarly to the
effective mass, for $a/a_{\rm end}\lesssim 160$, the rate is solely determined
by the anharmonic oscillation of the inflaton, $R(t)\propto
m_{\phi}\rho_{\phi}\propto(a/a_{\rm end})^{-5}$. However, at later times, the
parametric growth of non-vanishing $k$ modes seeps energy from the classical
$\phi$, reducing accordingly the corresponding rate. For the condensate
component, at late times, we find $R_{\phi}\propto m_{\phi}\mkern
1.5mu\overline{\mkern-1.5mu\rho_{\phi}\mkern-1.5mu}\mkern
1.5mu\propto(a/a_{\text{end}})^{-6.6}$. On the other hand, following the
orange curve, we see that the contribution to $R$ from the fragmented inflaton
becomes quickly important, and in fact dominates the dissipation process for
$a/a_{\rm end}\gtrsim 200$. Since at this stage the comoving number density of
$\delta\phi$ is conserved, we have $R_{\delta\phi}\propto
m_{\phi}^{2}n_{\delta\phi}\propto a^{-5.65}$.

###  3.2 Reheating temperatures

With the Boltzmann rate $R(t)$ at hand, we can now estimate the effect of
fragmentation on the instantaneous temperature of the primordial plasma during
reheating, and its value at the beginning of the domination by the thermal
bath, which we denote by $T_{\rm reh}$. For simplicity we assume the
instantaneous thermalization of the inflaton decay products, and for
definiteness we assume an effective number of degrees of freedom $g_{\rm
reh}=427/4$ above the electroweak scale, which is the Standard Model value.

####  3.2.1 Condensate decay

If we neglect the effect of backreaction, or assume a large coupling constant
so that the decay occurs prior to fragmentation, we can use the condensate
form for $R(t)$. Fixing for simplicity $a_{\rm end}=1$ for now, Eq. (3.3) can
be written as

$\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}a}\left(\rho_{R}a^{4}\right)\;=\;\frac{R(a)\,a^{3}}{H}\,,$
(3.16)

in general, and as

$\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}a}\left(\rho_{R}a^{4}\right)\;=\;\frac{\alpha^{2}y^{2}\lambda^{1/4}(\rho_{\phi}a^{4})^{5/4}M_{P}}{\pi\sqrt{(\rho_{\phi}+\rho_{R})a^{4}}}\;\simeq\;\frac{\alpha^{2}}{\pi}y^{2}\lambda^{1/4}\rho_{\rm
end}^{3/4}M_{P}\,,$ (3.17)

for $R=R_{\phi}$. In the second equality we have used the fact that, to a good
approximation, $\rho_{\phi}a^{4}\simeq\rho_{\rm end}\gg\rho_{R}a^{4}$ up until
the end of reheating. Straightforward integration then yields

$\displaystyle\rho_{R}\;$
$\displaystyle\simeq\;\frac{\alpha^{2}}{\pi}y^{2}\lambda^{1/4}\rho_{\rm
end}^{3/4}M_{P}\left(\frac{a_{\rm end}}{a}\right)^{3}\,,$ (3.18)
$\displaystyle T\;$ $\displaystyle=\;\left(\frac{30\rho_{R}}{\pi^{2}g_{\rm
reh}}\right)^{1/4}\,,$ (3.19)

that is, $T\propto a^{-3/4}$ during reheating. The end of reheating is
determined by the inflaton-radiation equality condition,
$\rho_{R}=\rho_{\phi}$. Solving with the above expression for $\rho_{R}$, one
obtains [11, 12]

$\displaystyle\frac{a_{\rm reh}}{a_{\rm end}}\;$
$\displaystyle\simeq\;\frac{\pi\rho_{\rm
end}^{1/4}}{\alpha^{2}y^{2}\lambda^{1/4}M_{P}}\,,$ (3.20) $\displaystyle
T_{\rm reh}\;$ $\displaystyle\simeq\;\left(\frac{30\lambda}{\pi^{6}g_{\rm
reh}}\right)^{1/4}\alpha^{2}y^{2}M_{P}\,.$ (3.21)

Figure 5: Scale factor dependence of the instantaneous temperature of the
relativistic plasma during and after reheating, for a selection of effective
inflaton-fermion couplings. Shown as dashed curves is the numerical
integration of Eq. (3.3) assuming no fragmentation of the inflaton condensate.
For those curves, the end of reheating is marked with a triangle. Continuous
curves correspond to the solution of (3.3) accounting for the fragmentation of
$\phi$. The end of reheating in this case is marked by a star.

Fig. 5 shows the evolution of the instantaneous temperature during and after
reheating in the condensate approximation, shown as the dashed lines, for
$y=\\{10^{-1},10^{-2},10^{-3}\\}$. This temperature rapidly rises from zero at
the end of inflation to reach a maximum $T_{\rm max}$, which may be
approximated as [12]

$T_{\rm max}\;\simeq\;6\times 10^{14}\,y^{1/2}\,{\rm GeV}\,.$ (3.22)

Soon after reaching this maximum temperature, the relativistic bath of
inflaton decay products continues being populated by these decays, while being
redshifted by expansion, following the relation (3.19). Finally, upon reaching
the end of reheating, given by (3.20), the thermal plasma simply redshifts,
$T\propto a^{-1}$ during the radiation dominated epoch.

Fig. 6 shows the reheating temperature $T_{\rm reh}$, as a function of the
effective coupling $y$, corresponding to the orange dashed curve in the pure
condensate approximation. In the range of couplings selected the reheating
temperature is always well above the lower bound imposed by successful big
bang nucleosynthesis (BBN), $T_{\rm BBN}\sim 1\,{\rm MeV}$ [72].

####  3.2.2 Fragmented inflaton decay

Taking now into account the preheating phase, we can approximate the total
rate $R\simeq R_{\delta\phi}$ at late times, $a/a_{\rm
end}\gg\mathcal{O}(10^{2})$. Parametrizing

$R_{\delta\phi}\;=\;\frac{y^{2}}{8\pi}m_{\phi}^{2}n_{\delta\phi}\;\equiv\;\frac{\gamma_{\phi}y^{2}}{8\sqrt{3}\pi}M_{P}^{5}\left(\frac{a_{\rm
end}}{a}\right)^{5+x}\,,$ (3.23)

with $\gamma_{\phi}\simeq 2.49\times 10^{-15}$ and $x\simeq 0.65$ (as
determined from the lattice, cf. Fig. 8), we can immediately integrate (3.16)
to obtain

$\rho_{R}\;\simeq\;\frac{\gamma_{\phi}y^{2}}{8\pi(1-x)}\left(\frac{M_{P}^{4}}{\rho_{\rm
end}}\right)^{1/2}M_{P}^{4}\left(\frac{a}{a_{\rm end}}\right)^{-3-x}\,,$
(3.24)

during reheating. This implies that the instantaneous temperature redshifts as
$T\propto a^{-3/4-x/4}\simeq a^{-0.91}$, that is, more rapidly than in the
pure-condensate scenario. In this case, reheating ends when

$\frac{a_{\rm reh}}{a_{\rm
end}}\;\simeq\;\left[\frac{8\pi(1-x)}{\gamma_{\phi}y^{2}}\left(\frac{\rho_{\rm
end}}{M_{P}^{4}}\right)^{3/2}\right]^{\frac{1}{1-x}}\,,$ (3.25)

and

$T_{\rm reh}\;\simeq\;\left(\frac{30\rho_{\rm end}}{\pi^{2}g_{\rm
reh}}\right)^{1/4}\left[\frac{\gamma_{\phi}y^{2}}{8\pi(1-x)}\left(\frac{M_{P}^{4}}{\rho_{\rm
end}}\right)^{3/2}\right]^{\frac{1}{1-x}}\,,$ (3.26)

which is the main result of this work. In particular, we note that due to
fragmentation, $T_{\rm reh}\propto y^{5.68}$. Therefore, low reheating
temperatures are expected even for comparatively large values of $y$ with
respect to the pure condensate scenario.

Figure 6: Reheating temperature as a function of the inflaton-fermion
coupling, accounting for $\phi$ fragmentation (continuous), and without
accounting for fragmentation (dashed). Shown in gray is the region forbidden
from big bang nucleosynthesis considerations.

The effect of the self-resonance of the inflaton on the evolution of the
instantaneous temperature of the relativistic plasma can be clearly seen in
Fig. 5. Here the continuous curves correspond to the fragmentation scenario.
Immediately after the beginning of reheating the temperature follows the
condensate result, as the bulk of $\rho_{\phi}$ is still contained in the
coherent component of $\phi$. However, after the onset of backreaction, the
condensate is depleted and the decay process is dominated by the decay of
$k\neq 0$ modes, with a rate that redshifts faster than that of the pure
condensate case (cf. Fig. 8). The result is a faster redshift of $T$, and
consequently a reduced reheating temperature, shown in the figure as stars.

The difference in reheating temperatures can be better appreciated in Fig. 6.
For $y\gtrsim 0.17$ the decay of $\phi$ occurs earlier than the onset of
backreaction, and the full result is indistinguishable from the pure
condensate scenario. On the other hand, for smaller couplings, $T_{\rm reh}$
is smaller relative to the condensate result, with the difference increasing
for decreasing $y$. The small kinks that may be noted in the $T_{\rm reh}$ vs
$y$ curve correspond to the effect of the change in the effective number of
relativistic degrees of freedom. We finally note that only for $y\gtrsim
2.8\times 10^{-4}$ can the BBN constrain be averted.

##  4 Gravitational waves

In cosmological perturbation theory, tensor metric perturbations are sourced
by scalar metric perturbations, starting at second order. This contribution,
sensibly negligible prior to fragmentation, would become significant once non-
linearity induces sizable mode-mode couplings, efficiently sourcing
gravitational waves. In this section we compute the gravitational wave
production associated to the large inhomogeneities generated posterior to the
end of inflation by the self-fragmentation of the inflaton condensate.
Accounting for tensor metric perturbations, the isotropic and homogeneous
Friedmann-Robertson-Walker metric is modified to

$\mathop{}\\!\mathrm{d}s^{2}\,=\,a(\tau)^{2}\Big{[}\mathop{}\\!\mathrm{d}\tau^{2}-\Big{(}\delta_{ij}+h_{ij}\Big{)}\mathop{}\\!\mathrm{d}x^{i}\mathop{}\\!\mathrm{d}x^{j}\Big{]}\,,$
(4.1)

where the transverseness and tracelessness (TT) conditions
$\partial_{i}h_{ij}=0$ and $h^{i}_{i}=0$ are satisfied. The transverse-
traceless component of the anisotropic stress
$\Pi^{\text{TT}}_{ij}=\big{[}\partial_{i}\phi\partial_{j}\phi\big{]}^{\text{TT}}$
appears as a source term in the equation of motion for the tensor metric
perturbations444$h_{ij}$ represent two independent tensor degrees of freedom
satisfying the same equation.

$h^{\prime\prime}_{ij}(\boldsymbol{k},\tau)+2\mathcal{H}h^{\prime}_{ij}(\boldsymbol{k},\tau)+k^{2}h_{ij}(\boldsymbol{k},\tau)\,=\,\dfrac{2}{M_{P}^{2}}\Pi^{\text{TT}}_{ij}(\boldsymbol{k},\tau)\,,$
(4.2)

where the anisotropic stress projected onto the transverse-traceless component
in Fourier space is [41]

$\Pi^{\text{TT}}_{ij}(\boldsymbol{k},\tau)\,=\left(\,P_{i\ell}(\hat{\boldsymbol{k}})P_{jm}(\hat{\boldsymbol{k}})-\dfrac{1}{2}P_{ij}(\hat{\boldsymbol{k}})P_{\ell
m}(\hat{\boldsymbol{k}})\right)\Pi_{\ell m}(\boldsymbol{k},\tau)\,,$ (4.3)

with $P_{ij}=\delta_{ij}-\hat{k_{i}}\hat{k_{j}}$ and $\hat{k}_{i}=k_{i}/k$ and

$\Pi_{\ell
m}(\boldsymbol{k},\tau)\,=\,\int\dfrac{\mathop{}\\!\mathrm{d}^{3}\boldsymbol{p}}{(2\pi)^{3/2}}\,p_{\ell}\,p_{m}\,\phi(\boldsymbol{p},\tau)\,\phi(\boldsymbol{k}-\boldsymbol{p},\tau)\,.$
(4.4)

###  4.1 Simulated spectrum

We use $\mathcal{C}\text{osmo}\mathcal{L}\text{attice}$ to simulate the
production of gravitational waves for our model. In practice, a TT projector
is implemented in a discretized version on the lattice. The non-uniqueness of
the definition of such a projector can induce small differences in the GW
spectrum. However, it has been argued it only marginally affects the UV part
of the GW spectrum [73, 60]. On the lattice, the GW energy density is computed
as a volume average

$\rho_{\text{GW}}\,=\,\dfrac{M_{P}^{2}}{4a^{4}}\dfrac{1}{V}\int\mathop{}\\!\mathrm{d}^{3}\boldsymbol{k}\,\bar{h}_{ij}^{\prime}(\tau,\boldsymbol{k})\,\bar{h}_{ij}^{\prime*}(\tau,\boldsymbol{k})\,,$
(4.5)

with $\bar{h}_{ij}=ah_{ij}$. This sum is performed by
$\mathcal{C}\text{osmo}\mathcal{L}\text{attice}$, after discretization, over a
finite lattice volume $V$ (details can be found in Ref. [61]). The normalized
GW energy density per logarithmic frequency interval can be expressed as

$\Omega_{\text{GW}}(k)\,=\,\dfrac{1}{\rho_{c}}\dfrac{\mathop{}\\!\mathrm{d}\rho_{\text{GW}}}{\mathop{}\\!\mathrm{d}\log
k}\,,$ (4.6)

where $\rho_{\text{GW}}$ is the total GW energy density and $\rho_{c}$ is the
(time-dependent) critical energy density defined via
$\rho_{c}=3H^{2}M_{P}^{2}$. The GW frequency $f$ at the present epoch can be
related to the comoving Fourier scale $k$ via [41]

$f\,=\,\left(\dfrac{a_{\text{end}}}{a_{0}}\right)\dfrac{k}{2\pi}\,=\,\left(\dfrac{\rho_{\text{rad},0}}{\rho_{\text{end}}}\right)^{1/4}\left(\dfrac{g_{\text{reh}}}{g_{0}}\right)^{1/4}\left(\dfrac{g_{s,0}}{g_{s,\text{reh}}}\right)^{1/3}\dfrac{k}{2\pi}\,\simeq\,1.46\times
10^{8}\,\left(\dfrac{k}{0.7\,m_{\text{end}}}\right)\leavevmode\nobreak\
\text{Hz}\,,$ (4.7)

where we normalized the scale $k$ to the resonant value from the Floquet
analysis $k\simeq 0.7\,m_{\text{end}}$. $\rho_{\text{rad},0}$ is the radiation
energy density at the present time. The GW energy density estimated from our
lattice simulation is shown in Fig. 7, evaluated a different times on the left
panel, as a function of the frequency, and integrated over frequency as a
function of time on the right panel. From this figure one can see that the GW
energy density significantly increases at around $\tau/\tau_{\text{end}}\simeq
60-80$ when scalar inhomogeneities, triggered by the parametric resonance
discussed in Sec. 2.2, start to increase, as seen in Fig. 3. The GW energy
density relative to the critical energy density stabilises at around
$\tau/\tau_{\text{end}}\simeq 175$ when fragmentation is achieved, asymptoting
smoothly towards a constant value
$\rho_{\text{GW}}/\rho_{c}\,\simeq\,7.3\times 10^{-6}$ at larger
$\tau/\tau_{\text{end}}$.

Figure 7: Left: Differential gravitational wave energy density per logarithmic
(present day) frequency interval, normalized to the critical density,
evaluated at selected conformal time. The color code corresponds to the
evaluation time of the right panel. Right: Integrated differential GW energy
density evolution with respect to conformal time, normalized to the critical
energy density.

One of the most striking features of the final GW spectrum, i.e. the red curve
in Fig. 7 on the left panel, is the peak structure. 555Such structure was
already identified for a quartic potential in [33]. In the red curve, the
dominant peak is located at a frequency $f\simeq 1.9\times
10^{8}\leavevmode\nobreak\ \text{Hz}$ slightly larger than the predicted
resonant comoving scale $k/m_{\text{end}}\simeq 0.7$ from the Floquet analysis
from Eq. (4.7). From the left panel of Fig. 7, one can see that the initial GW
spectra dominant peak was located at a slightly lower frequency and
progressively displaced towards higher frequencies, explaining such
difference. Disregarding the smaller peak around $\sim
10^{8}\leavevmode\nobreak\ \text{Hz}$ that seems to appear at the onset of
fragmentation (originated by non-linear effects), we can compare the peak
structure of the spectrum to those of the scalar spectra and predictions from
the Boltzmann approach, cf. Eq. (2.47). A selection of numerical results are
presented in Table 1, and are shown graphically in Fig. 8 for the inflaton
fluctuation PSD and the gravitational wave spectrum. Even if the peak
structure of the scalar spectra are smoothed out by the end of reheating, the
structure of the peaks remains imprinted in the tensor spectra. The Boltzmann
approach allows for a very decent prediction for the location of the peaks in
the scalar spectra. More remarkably, the Boltzmann approach allows to explain
the spacing between the first peaks with a precision at the $<10\%$ level.

Figure 8: Comparison between the peak structure of the scalar phase space distribution at the onset of fragmentation (left), and of the final gravitational wave spectrum (right), with the predictions of the Boltzmann approach (vertical gray lines). The Boltzmann lines are only indicative of the value of the comoving momenta, the ordinate values are matched with the corresponding black curve. | $\hat{q}_{1}$ | $\hat{q}_{1}/\hat{q}_{2}$ | $\hat{q}_{1}/\hat{q}_{3}$ | $\hat{q}_{1}/\hat{q}_{4}$ | $\hat{q}_{2}/\hat{q}_{4}$ | $\hat{q}_{2}/\hat{q}_{5}$
---|---|---|---|---|---|---
Lattice (scalar) | 0.7 | 0.54 | 0.48 | 0.3 | 0.56 | 0.45
Lattice (tensor) | 0.91 | 0.65 | 0.48 | 0.38 | 0.59 | 0.48
Boltzmann | 0.7 | 0.65 | 0.51 | 0.42 | 0.63 | 0.54
Linear/Hartree | 0.7 | - | - | - | - | -
Floquet | 0.7 | - | - | - | - | -
Table 1: Comparison of the peak structure of the scalar and tensor spectra
where $\hat{q}_{i}$ corresponds to the location of the $i^{\text{th}}$ peak in
terms of the rescaled momentum $q\equiv k/m_{\text{end}}$. The Boltzmann
predictions are based on Eq. (2.46) and Eq. (2.47).

###  4.2 Signal and detection prospects

There are essentially two distinct GW signals corresponding to two different
physical processes predicted by this model. First, during inflation quantum
fluctuations of tensor modes are stretched on macroscopic super-horizon
scales. Once such scales re-enter the horizon in a subsequent phase of the
universe, tensor perturbations manifest as a stochastic background of
gravitational waves. Primordial tensor perturbations, customarily parametrized
in terms of the tensor-to-scalar ratio can be expressed for the potential of
Eq. (2.2) as [12]

$r\,\simeq\,\dfrac{12}{N_{*}^{2}}\,\simeq\,3.8\times
10^{-3}\,\left(\dfrac{56}{N_{*}}\right)^{2}\,,$ (4.8)

evaluated at the CMB fiducial scale $k_{*}=0.05\leavevmode\nobreak\
\text{Mpc}^{-1}$. Such value for the tensor-to-scalar ratio appears to be
close to the sensitivity reach for the upcoming Simons Observatory (SO)
$r<6\times 10^{-3}$ [74]. Moreover, future missions such as LiteBIRD and CMB
Stage-4 (CMB-S4) should reach sensitivities as small as $r<2\times 10^{-3}$
[75] and $r<10^{-3}$ [76] respectively, allowing to disprove or confirm this
model.

Figure 9: Simulated amplitude of the present-day GW signal amplitude as a
function of the frequency (solid red line) and sensitivity estimate for
resonant electromagnetic cavities of Ref. [44] (orange dashed line).

The second GW signal is generated during fragmentation from large inflaton
inhomogeneities which we simulated using
$\mathcal{C}\text{osmo}\mathcal{L}\text{attice}$. We can estimate the current
GW energy density from

$\rho_{\text{GW,0}}\,=\,\rho_{\text{GW,reh}}\left(\dfrac{a_{\text{reh}}}{a_{0}}\right)^{4}\,=\,\rho_{\text{GW,reh}}\left(\dfrac{\rho_{\text{rad},0}}{\rho_{c,\text{reh}}}\right)\left(\dfrac{g_{\text{reh}}}{g_{0}}\right)\left(\dfrac{g_{s,0}}{g_{s,\text{reh}}}\right)^{4/3}\,,$
(4.9)

where we took $g_{\text{reh}}=g_{s,\text{reh}}=106.75$, $g_{s,0}=3.909$ and
$g_{0}=3.363$. The present day normalized GW energy density can thus be
expressed as

$\Omega_{\text{GW},0}(f)\,=\,\dfrac{\rho_{\text{rad},0}}{\rho_{c,0}}\left(\dfrac{g_{\text{reh}}}{g_{0}}\right)\left(\dfrac{g_{s,0}}{g_{s,\text{reh}}}\right)^{4/3}\,\Omega_{\text{GW,reh}}(f)\,\simeq\,\,2\times
10^{-5}\,\Omega_{\text{GW,reh}}(f)\,.$ (4.10)

Such gravitational wave signal would induce a deviation to the effective
number of relativistic species $\Delta N_{\text{eff}}\equiv
N_{\text{eff}}-N_{\text{eff}}^{\text{SM}}$ with
$N_{\text{eff}}^{\text{SM}}=3.046$. This contribution can be expressed as

$\Delta
N_{\text{eff}}\,=\,\dfrac{8}{7}\left(\dfrac{11}{4}\right)^{4/3}\left(\dfrac{\rho_{\text{GW,0}}}{\rho_{\text{rad},0}}\right)=\,\dfrac{8}{7}\left(\dfrac{11}{4}\right)^{4/3}\left(\dfrac{g_{\text{reh}}}{g_{0}}\right)\left(\dfrac{g_{s,0}}{g_{s,\text{reh}}}\right)^{4/3}\left(\dfrac{\rho_{\text{GW,reh}}}{\rho_{c,\text{reh}}}\right)\,\simeq\,1.7\left(\dfrac{\rho_{\text{GW,reh}}}{\rho_{c,\text{reh}}}\right)\,.$
(4.11)

We find that for $\tau/\tau_{\text{end}}>300$, as illustrated on the right
panel of Fig. 7, the produced GW energy density asymptotes during reheating to
the value $\rho_{\text{GW}}/\rho_{c}\,\simeq\,7.3\times 10^{-6}$. This allows
to estimate the contribution to the effective number of non-relativistic
degrees of freedom

$\Delta N_{\text{eff}}\,\simeq\,1.2\times 10^{-5}\,,$ (4.12)

which is far below the current sensitivity from Planck
$N_{\text{eff}}=2.98^{+0.39}_{-0.38}$ at $95\%$ confidence level (TT, TE,
EE+lowE+lensing+BAO) [46] as well as future CMB missions CMB-S4 [77] and CMB-
HD [78] with an expected precision for $\Delta N_{\text{eff}}$ of $0.06$ and
$0.027$ respectively at $95\%$ confidence level. In order to facilitate
comparison with sensitivity prospects for future experiments, we can define
the dimensionless characteristic strain as

$h_{c}(f)\,=\,\sqrt{\dfrac{3H_{0}^{2}\Omega_{\text{GW,0}}(f)}{2\pi^{2}}}\dfrac{1}{f}\,.$
(4.13)

In Fig. 9 we compare our estimate of the present day strain as a function of
the frequency, in comparison to the sensitivity estimate for high-frequency
resonant electromagnetic cavities of Ref. [44]. The frequency range of the GW
signal emitted from the fragmentation process falls in the appropriate
sensitivity range of the proposal of Ref. [44]. Remarkably, the reconstitution
of the peak structure of such high-frequency GW signal, in addition to
measurements of the tensor-to-scalar ratio would help break the ambiguities in
determining the realization of inflation.

##  5 Summary and conclusions

In this work we have investigated the post-inflationary dynamics of the
inflaton field and its decay products, under the assumption of a potential
with a quartic minimum. For definiteness, we considered the T-model of
inflation, compatible with current constraints on inflation. After inflation,
the inflaton condensate coherently oscillates about its minimum. Triggered by
parametric resonance effects sourced by the oscillating background,
exponential growth of inflaton inhomogeneities result in the fragmentation of
this condensate. In this work we explored the role and consequences of such
non-linear effects on the post-inflationary history and the successful
reheating of the universe. We summarize in the following the key aspects and
results of our analysis.

Parametric resonances and fragmentation. First, we characterized the
parametric resonant effects triggering fragmentation using a standard Floquet
approach. We found that due to the conformal symmetry of the system, only
Fourier scales located at $k/m_{\text{end}}\simeq 0.7$ would experience a
significant exponential growth. By solving numerically the equation of motion
for the inflaton fluctuations at linear order in perturbation theory, we were
able to observe a growth of modes with $k/m_{\text{end}}\simeq 0.7$, therefore
confirming predictions from the Floquet analysis. To go beyond the linear
analysis, we simulated the evolution of the space-time dependent inflaton
field configuration with the nonperturbative code
$\mathcal{C}\text{osmo}\mathcal{L}\text{attice}$. By computing the occupation
number for the inflaton perturbations, we recovered results from the linear
approach in the first instants after the end of inflation. We find that at
around 5 $e$-folds after the end of inflation, perturbations enter the non-
linear regime, efficiently backreacting on the inflation condensate, leading
to fragmentation.

Before rescatterings fully redistribute the energy density into other modes,
the occupation number of inflaton quanta features several peaks located at
higher momentum than the resonant $k/m_{\text{end}}\simeq 0.7$ Floquet mode.
We estimated analytically the position of these peaks by means of a Boltzmann
approach. Notably, despite the validity of this formalism only in the linear
regime, we find that it not only accurately predicts the dominant Floquet
peak, but it also matches the lattice peak location within a few percent, at
least before they are fully washed-out by rescatterings.

Importantly, we found that the redistribution of the inflaton energy density
does not fully erase the coherent zero mode. Upon the onset of fragmentation,
the energy budget of the universe becomes dominated by a collection of
inflaton quanta redshifting as radiation $\rho_{\delta\phi}\propto a^{-4}$,
but a leftover inflaton condensate persists with $\mkern
1.5mu\overline{\mkern-1.5mu\rho_{\phi}\mkern-1.5mu}\mkern 1.5mu\propto
a^{-5.3}$.

Reheating. The main result of this work consists in the exploration of the
consequences of fragmentation to successfully achieve reheating; that is, the
transition to a universe dominated by radiation in thermal equilibrium. By
using a Boltzmann approach, we computed the contributions to the production
rates of relativistic fermionic states coming from the oscillating inflaton
condensate $R_{\phi}$, and also from the fragmentated population of inflaton
quanta $R_{\delta\phi}$. Prior to backreaction, the rate at which energy is
injected into the thermal bath per unit volume per unit time, scales as
$R_{\phi}\propto a^{-5}$. After fragmentation, we find that the condensate
rate $R_{\phi}\propto a^{-6.6}$, and the inflaton particle rate
$R_{\delta\phi}\propto a^{-5.65}$, both redshift faster, due to the rapid
decrease in the time-dependent induced inflaton effective mass. As a
consequence, reheating is less efficient after fragmentation. The
instantaneous temperature falls as $T\propto a^{-0.91}$ (compared to $T\propto
a^{-3/4}$ without backreaction). For fermionic decays, the reheating
temperature scales as $T_{\rm reh}\propto y^{2}$ for $y\gtrsim 10^{-1}$, and
as $T_{\rm reh}\propto y^{5.68}$ for $y\lesssim 10^{-1}$, where $y$ is the
inflaton-matter Yukawa coupling. Requiring the reheating temperature to be
larger than the BBN temperature $T_{\text{BBN}}\sim\,\text{MeV}$ necessitates
couplings larger than $y>2.8\times 10^{-4}$.

We must emphasize that these results are valid only for a quartic inflaton
minimum. Nevertheless, the lattice+Boltzmann formalism developed in this work
is general, and can be applied to a more diverse collection of inflation
models, with non-quartic potentials. The precise details of the evolution of
energy densities and temperatures will vary on a model-by-model basis, and
will be the subject of future follow-up work.

Tensor perturbations. The large mode-mode couplings for the inflaton
inhomogeneities inducing fragmentation can also source sizable tensor
perturbations. We estimated the production of gravitational waves from a
simulation with $\mathcal{C}\text{osmo}\mathcal{L}\text{attice}$. The
gravitational wave spectrum features several peaks located at high frequencies
$f\sim 10^{8}-10^{9}\leavevmode\nobreak\ \text{Hz}$ whose precise locations
are inherited from the scalar spectrum. The ratios of the location for the
various peaks matches relatively well predictions from the Boltzmann approach.
We find that the GW energy-density contribution to the effective number of
relativistic species $\Delta N_{\text{eff}}\sim 10^{-5}$ is beyond the reach
of any upcoming experiments. However, the frequency range for the GW spectrum
appears in sensitivity-reach estimate for resonant electromagnetic cavities
[44].

## Acknowledgments

We would like to thank Keith Olive, Yann Mambrini, Andreas Ringwald and
Sarunas Verner for helpful discussions. MG is supported by the DGAPA-PAPIIT
grant IA103123 at UNAM, and the CONAHCYT “Ciencia de Frontera” grant
CF-2023-I-17. MP acknowledges support by the Deutsche Forschungsgemeinschaft
(DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC
2121 “Quantum Universe” – 390833306. This work was made possible by the
support of the Institut Pascal at Universite Paris-Saclay during the Paris-
Saclay Astroparticle Symposium 2022, with the support of the P2IO Laboratory
of Excellence (program “Investissements d’avenir” ANR-11-IDEX-0003-01 Paris-
Saclay and ANR-10-LABX-0038), the P2I axis of the Graduate School Physics of
Universite Paris-Saclay, as well as IJCLab, CEA, IPhT, APPEC, the IN2P3 master
projet UCMN and EuCAPT ANR-11-IDEX-0003-01 Paris-Saclay and ANR-10-LABX-0038).
Non-perturbative numerical results in the Hartree approximation were obtained
from a custom Fortran code utilizing the thread-safe arbitrary precision
package MPFUN-For [79].

## References

  * [1] A. H. Guth, “The Inflationary Universe: A Possible Solution to the Horizon and Flatness Problems,” Phys. Rev. D 23 (1981) 347–356.
  * [2] K. Sato, “First Order Phase Transition of a Vacuum and Expansion of the Universe,” Mon. Not. Roy. Astron. Soc. 195 (1981) 467–479.
  * [3] A. H. Guth and E. J. Weinberg, “Could the Universe Have Recovered from a Slow First Order Phase Transition?,” Nucl. Phys. B 212 (1983) 321–364.
  * [4] K. A. Olive, “Inflation,” Phys. Rept. 190 (1990) 307–403.
  * [5] G. N. Felder, L. Kofman, and A. D. Linde, “Instant preheating,” Phys. Rev. D 59 (1999) 123523, arXiv:hep-ph/9812289.
  * [6] G. N. Felder, L. Kofman, and A. D. Linde, “Inflation and preheating in NO models,” Phys. Rev. D 60 (1999) 103505, arXiv:hep-ph/9903350.
  * [7] Y. Shtanov, J. H. Traschen, and R. H. Brandenberger, “Universe reheating after inflation,” Phys. Rev. D 51 (1995) 5438–5455, arXiv:hep-ph/9407247.
  * [8] K. Ichikawa, T. Suyama, T. Takahashi, and M. Yamaguchi, “Primordial Curvature Fluctuation and Its Non-Gaussianity in Models with Modulated Reheating,” Phys. Rev. D 78 (2008) 063545, arXiv:0807.3988 [astro-ph].
  * [9] S. Nurmi, T. Tenkanen, and K. Tuominen, “Inflationary Imprints on Dark Matter,” JCAP 11 (2015) 001, arXiv:1506.04048 [astro-ph.CO].
  * [10] K. Kainulainen, S. Nurmi, T. Tenkanen, K. Tuominen, and V. Vaskonen, “Isocurvature Constraints on Portal Couplings,” JCAP 06 (2016) 022, arXiv:1601.07733 [astro-ph.CO].
  * [11] M. A. G. Garcia, K. Kaneta, Y. Mambrini, and K. A. Olive, “Reheating and Post-inflationary Production of Dark Matter,” Phys. Rev. D 101 no. 12, (2020) 123507, arXiv:2004.08404 [hep-ph].
  * [12] M. A. G. Garcia, K. Kaneta, Y. Mambrini, and K. A. Olive, “Inflaton Oscillations and Post-Inflationary Reheating,” JCAP 04 (2021) 012, arXiv:2012.10756 [hep-ph].
  * [13] A. D. Dolgov and D. P. Kirilova, “On Particle Creation by a Time Dependent Scalar Field,” Sov. J. Nucl. Phys. 51 (1990) 172–177.
  * [14] J. H. Traschen and R. H. Brandenberger, “Particle Production During Out-of-equilibrium Phase Transitions,” Phys. Rev. D 42 (1990) 2491–2504.
  * [15] L. Kofman, A. D. Linde, and A. A. Starobinsky, “Reheating after inflation,” Phys. Rev. Lett. 73 (1994) 3195–3198, arXiv:hep-th/9405187.
  * [16] D. Boyanovsky, M. D’Attanasio, H. J. de Vega, R. Holman, D. S. Lee, and A. Singh, “Reheating the postinflationary universe,” arXiv:hep-ph/9505220.
  * [17] M. Yoshimura, “Catastrophic particle production under periodic perturbation,” Prog. Theor. Phys. 94 (1995) 873–898, arXiv:hep-th/9506176.
  * [18] L. Kofman, A. D. Linde, and A. A. Starobinsky, “Towards the theory of reheating after inflation,” Phys. Rev. D 56 (1997) 3258–3295, arXiv:hep-ph/9704452.
  * [19] J. Garcia-Bellido, M. Garcia Perez, and A. Gonzalez-Arroyo, “Symmetry breaking and false vacuum decay after hybrid inflation,” Phys. Rev. D 67 (2003) 103501, arXiv:hep-ph/0208228.
  * [20] G. N. Felder and L. Kofman, “Nonlinear inflaton fragmentation after preheating,” Phys. Rev. D 75 (2007) 043518, arXiv:hep-ph/0606256.
  * [21] A. V. Frolov, “Non-linear Dynamics and Primordial Curvature Perturbations from Preheating,” Class. Quant. Grav. 27 (2010) 124006, arXiv:1004.3559 [gr-qc].
  * [22] M. A. Amin, M. P. Hertzberg, D. I. Kaiser, and J. Karouby, “Nonperturbative Dynamics Of Reheating After Inflation: A Review,” Int. J. Mod. Phys. D 24 (2014) 1530003, arXiv:1410.3808 [hep-ph].
  * [23] M. A. G. Garcia, K. Kaneta, Y. Mambrini, K. A. Olive, and S. Verner, “Freeze-in from preheating,” JCAP 03 no. 03, (2022) 016, arXiv:2109.13280 [hep-ph].
  * [24] M. A. G. Garcia, M. Pierre, and S. Verner, “Scalar dark matter production from preheating and structure formation constraints,” Phys. Rev. D 107 no. 4, (2023) 043530, arXiv:2206.08940 [hep-ph].
  * [25] P. B. Greene, L. Kofman, A. D. Linde, and A. A. Starobinsky, “Structure of resonance in preheating after inflation,” Phys. Rev. D 56 (1997) 6175–6192, arXiv:hep-ph/9705347.
  * [26] D. I. Kaiser, “Preheating in an expanding universe: Analytic results for the massless case,” Phys. Rev. D 56 (1997) 706–716, arXiv:hep-ph/9702244.
  * [27] J. Garcia-Bellido, D. G. Figueroa, and J. Rubio, “Preheating in the Standard Model with the Higgs-Inflaton coupled to gravity,” Phys. Rev. D 79 (2009) 063531, arXiv:0812.4624 [hep-ph].
  * [28] M. A. Amin, R. Easther, H. Finkel, R. Flauger, and M. P. Hertzberg, “Oscillons After Inflation,” Phys. Rev. Lett. 108 (2012) 241302, arXiv:1106.3335 [astro-ph.CO].
  * [29] M. P. Hertzberg, J. Karouby, W. G. Spitzer, J. C. Becerra, and L. Li, “Theory of self-resonance after inflation. I. Adiabatic and isocurvature Goldstone modes,” Phys. Rev. D 90 (2014) 123528, arXiv:1408.1396 [hep-th].
  * [30] D. G. Figueroa and F. Torrenti, “Parametric resonance in the early Universe—a fitting analysis,” JCAP 02 (2017) 001, arXiv:1609.05197 [astro-ph.CO].
  * [31] K. D. Lozanov and M. A. Amin, “Equation of State and Duration to Radiation Domination after Inflation,” Phys. Rev. Lett. 119 no. 6, (2017) 061301, arXiv:1608.01213 [astro-ph.CO].
  * [32] K. D. Lozanov and M. A. Amin, “Self-resonance after inflation: oscillons, transients and radiation domination,” Phys. Rev. D 97 no. 2, (2018) 023533, arXiv:1710.06851 [astro-ph.CO].
  * [33] C. Fu, P. Wu, and H. Yu, “Production of gravitational waves during preheating with nonminimal coupling,” Phys. Rev. D 97 no. 8, (2018) 081303, arXiv:1711.10888 [gr-qc].
  * [34] S. Antusch, D. G. Figueroa, K. Marschall, and F. Torrenti, “Characterizing the postinflationary reheating history: Single daughter field with quadratic-quadratic interaction,” Phys. Rev. D 105 no. 4, (2022) 043532, arXiv:2112.11280 [astro-ph.CO].
  * [35] O. Lebedev, T. Solomko, and J.-H. Yoon, “Dark matter production via a non-minimal coupling to gravity,” JCAP 02 (2023) 035, arXiv:2211.11773 [hep-ph].
  * [36] C. Cosme, D. G. Figueroa, and N. Loayza, “Gravitational wave production from preheating with trilinear interactions,” Journal of Cosmology and Astroparticle Physics 2023 no. 05, (May, 2023) 023. https://doi.org/10.1088%2F1475-7516%2F2023%2F05%2F023.
  * [37] K. Alam, M. Bastero-Gil, K. Dutta, and H. V. Ragavendra, “Non-thermal moduli production during preheating in $\alpha$-attractor inflation models,” arXiv:2303.17383 [astro-ph.CO].
  * [38] E. J. Copeland, M. Gleiser, and H. R. Muller, “Oscillons: Resonant configurations during bubble collapse,” Phys. Rev. D 52 (1995) 1920–1933, arXiv:hep-ph/9503217.
  * [39] P. Salmi and M. Hindmarsh, “Radiation and Relaxation of Oscillons,” Phys. Rev. D 85 (2012) 085033, arXiv:1201.1934 [hep-th].
  * [40] M. A. Amin, R. Easther, and H. Finkel, “Inflaton Fragmentation and Oscillon Formation in Three Dimensions,” JCAP 12 (2010) 001, arXiv:1009.2505 [astro-ph.CO].
  * [41] J. F. Dufaux, A. Bergman, G. N. Felder, L. Kofman, and J.-P. Uzan, “Theory and Numerics of Gravitational Waves from Preheating after Inflation,” Phys. Rev. D 76 (2007) 123517, arXiv:0707.0875 [astro-ph].
  * [42] J. Garcia-Bellido and D. G. Figueroa, “A stochastic background of gravitational waves from hybrid preheating,” Phys. Rev. Lett. 98 (2007) 061302, arXiv:astro-ph/0701014.
  * [43] A. Ringwald and C. Tamarit, “Revealing the cosmic history with gravitational waves,” Phys. Rev. D 106 no. 6, (2022) 063027, arXiv:2203.00621 [hep-ph].
  * [44] N. Herman, L. Lehoucq, and A. Fúzfa, “Electromagnetic Antennas for the Resonant Detection of the Stochastic Gravitational Wave Background,” arXiv:2203.15668 [gr-qc].
  * [45] D. G. Figueroa and F. Torrenti, “Gravitational wave production from preheating: parameter dependence,” JCAP 10 (2017) 057, arXiv:1707.04533 [astro-ph.CO].
  * [46] Planck Collaboration, N. Aghanim et al., “Planck 2018 results. VI. Cosmological parameters,” Astron. Astrophys. 641 (2020) A6, arXiv:1807.06209 [astro-ph.CO]. [Erratum: Astron.Astrophys. 652, C4 (2021)].
  * [47] Planck Collaboration, Y. Akrami et al., “Planck 2018 results. X. Constraints on inflation,” Astron. Astrophys. 641 (2020) A10, arXiv:1807.06211 [astro-ph.CO].
  * [48] BICEP, Keck Collaboration, P. A. R. Ade et al., “Improved Constraints on Primordial Gravitational Waves using Planck, WMAP, and BICEP/Keck Observations through the 2018 Observing Season,” Phys. Rev. Lett. 127 no. 15, (2021) 151301, arXiv:2110.00483 [astro-ph.CO].
  * [49] R. Kallosh and A. Linde, “Universality Class in Conformal Inflation,” JCAP 07 (2013) 002, arXiv:1306.5220 [hep-th].
  * [50] A. R. Liddle and S. M. Leach, “How long before the end of inflation were observable perturbations produced?,” Phys. Rev. D 68 (2003) 103503, arXiv:astro-ph/0305263.
  * [51] J. Martin and C. Ringeval, “First CMB Constraints on the Inflationary Reheating Temperature,” Phys. Rev. D 82 (2010) 023511, arXiv:1004.5525 [astro-ph.CO].
  * [52] D. J. Fixsen, “The Temperature of the Cosmic Microwave Background,” Astrophys. J. 707 (2009) 916–920, arXiv:0911.1955 [astro-ph.CO].
  * [53] Y. Nambu and A. Taruya, “Evolution of cosmological perturbation in reheating phase of the universe,” Prog. Theor. Phys. 97 (1997) 83–89, arXiv:gr-qc/9609029.
  * [54] B. A. Bassett, D. I. Kaiser, and R. Maartens, “General relativistic preheating after inflation,” Phys. Lett. B 455 (1999) 84–89, arXiv:hep-ph/9808404.
  * [55] B. A. Bassett, F. Tamburini, D. I. Kaiser, and R. Maartens, “Metric preheating and limitations of linearized gravity. 2.,” Nucl. Phys. B 561 (1999) 188–240, arXiv:hep-ph/9901319.
  * [56] K. Jedamzik, M. Lemoine, and J. Martin, “Collapse of Small-Scale Density Perturbations during Preheating in Single Field Inflation,” JCAP 09 (2010) 034, arXiv:1002.3039 [astro-ph.CO].
  * [57] Z. Huang, “The Art of Lattice and Gravity Waves from Preheating,” Phys. Rev. D 83 (2011) 123509, arXiv:1102.0227 [astro-ph.CO].
  * [58] J. T. Giblin and A. J. Tishue, “Preheating in Full General Relativity,” Phys. Rev. D 100 no. 6, (2019) 063543, arXiv:1907.10601 [gr-qc].
  * [59] W. Magnus and S. Winkler, Hill’s Equation.  Dover Books on Mathematics. Dover Publications, Mineola, NY, Feb., 2004\.
  * [60] D. G. Figueroa, A. Florio, F. Torrenti, and W. Valkenburg, “The art of simulating the early Universe – Part I,” JCAP 04 (2021) 035, arXiv:2006.15122 [astro-ph.CO].
  * [61] D. G. Figueroa, A. Florio, F. Torrenti, and W. Valkenburg, “CosmoLattice,” arXiv:2102.01031 [astro-ph.CO].
  * [62] T. Moroi and W. Yin, “Light Dark Matter from Inflaton Decay,” JHEP 03 (2021) 301, arXiv:2011.09475 [hep-ph].
  * [63] T. Moroi and W. Yin, “Particle Production from Oscillating Scalar Field and Consistency of Boltzmann Equation,” JHEP 03 (2021) 296, arXiv:2011.12285 [hep-ph].
  * [64] A. Ghosh and S. Mukhopadhyay, “Momentum distribution of dark matter produced in inflaton decay: Effect of inflaton mediated scatterings,” Phys. Rev. D 106 no. 4, (2022) 043519, arXiv:2205.03440 [hep-ph].
  * [65] M. A. G. Garcia and M. A. Amin, “Prethermalization production of dark matter,” Phys. Rev. D 98 no. 10, (2018) 103504, arXiv:1806.01865 [hep-ph].
  * [66] G. Ballesteros, M. A. G. Garcia, and M. Pierre, “How warm are non-thermal relics? Lyman-$\alpha$ bounds on out-of-equilibrium dark matter,” JCAP 03 (2021) 101, arXiv:2011.13458 [hep-ph].
  * [67] P. B. Greene and L. Kofman, “Preheating of fermions,” Phys. Lett. B 448 (1999) 6–12, arXiv:hep-ph/9807339.
  * [68] G. F. Giudice, M. Peloso, A. Riotto, and I. Tkachev, “Production of massive fermions at preheating and leptogenesis,” JHEP 08 (1999) 014, arXiv:hep-ph/9905242.
  * [69] P. B. Greene and L. Kofman, “On the theory of fermionic preheating,” Phys. Rev. D 62 (2000) 123516, arXiv:hep-ph/0003018.
  * [70] M. Peloso and L. Sorbo, “Preheating of massive fermions after inflation: Analytical results,” JHEP 05 (2000) 016, arXiv:hep-ph/0003045.
  * [71] H. P. Nilles, M. Peloso, and L. Sorbo, “Coupled fields in external background with application to nonthermal production of gravitinos,” JHEP 04 (2001) 004, arXiv:hep-th/0103202.
  * [72] Particle Data Group Collaboration, R. L. Workman et al., “Review of Particle Physics,” PTEP 2022 (2022) 083C01.
  * [73] D. G. Figueroa, J. Garcia-Bellido, and A. Rajantie, “On the Transverse-Traceless Projection in Lattice Simulations of Gravitational Wave Production,” JCAP 11 (2011) 015, arXiv:1110.0337 [astro-ph.CO].
  * [74] Simons Observatory Collaboration, P. Ade et al., “The Simons Observatory: Science goals and forecasts,” JCAP 02 (2019) 056, arXiv:1808.07445 [astro-ph.CO].
  * [75] LiteBIRD Collaboration, E. Allys et al., “Probing Cosmic Inflation with the LiteBIRD Cosmic Microwave Background Polarization Survey,” arXiv:2202.02773 [astro-ph.IM].
  * [76] The Simons Observatory Collaboration, “The simons observatory: Astro2020 decadal project whitepaper,”. https://arxiv.org/abs/1907.08284.
  * [77] K. Abazajian et al., “CMB-S4 Science Case, Reference Design, and Project Plan,” arXiv:1907.04473 [astro-ph.IM].
  * [78] CMB-HD Collaboration, S. Aiola et al., “Snowmass2021 CMB-HD White Paper,” arXiv:2203.05728 [astro-ph.CO].
  * [79] D. H. Bailey, “MPFUN2015: A thread-safe arbitrary precision package,”. http://www.davidhbailey.com/dhbpapers/mpfun2015.pdf.