# Tachyonic Preheating in Plateau Inflation

Eemeli Tomberg  and Hardi Veermäe

###### Abstract

Plateau inflation is an experimentally consistent framework in which the scale
of inflation can be kept relatively low. Close to the edge of the plateau,
scalar perturbations are subject to a strong tachyonic instability. Tachyonic
preheating is realized when, after inflation, the oscillating inflaton
repeatedly re-enters the plateau. We develop the analytic theory of this
process and expand the linear approach by including backreaction between the
coherent background and growing perturbations. For a family of plateau models,
the analytic predictions are confronted with numerical estimates. Our analysis
shows that the inflaton fragments in a fraction of an $e$-fold in all examples
supporting tachyonic preheating, generalizing the results of previous similar
studies. In these scenarios, the scalar-to-tensor ratio is tiny, $r<10^{-7}$.

##  1 Introduction

Cosmic inflation can solve many outstanding issues of the Big Bang model and
successfully predict the spectrum of primordial perturbations [1, 2, 3, 4, 5,
6, 7, 8]. The latest cosmic microwave background (CMB) measurements by the
Planck satellite [8] set tight limits on primordial perturbations, especially
on the spectral index $n_{s}$ and the tensor-to-scalar ratio $r$. The
measurements of $n_{s}$ support plateau inflation, a class of models where the
inflaton field potential has a long, almost flat section at large field values
[9, 8]. Tensor modes are yet to be observed, thus only an upper bound
$r_{0.002}<0.056$ exists for the tensor-to-scalar ratio [8, 10]. Plateau
inflation can accommodate extremely low values of $r$ practically undetectable
even in next-generation CMB experiments, which are projected to have
sensitivities of the order $r\sim 0.001$ [11, 12, 13]. If the upper limit on
$r$ keeps shrinking, plateau models become more and more favoured.

In the hot Big Bang model, the post-inflationary universe must transition into
radiation domination. This transition process is dubbed reheating [14, 15],
and, in typical models, it proceeds through the inflaton field oscillating
around the minimum of its potential and producing relativistic particles
through non-perturbative processes and perturbative decay. The non-
perturbative process, called preheating, often proceeds through parametric
resonance, exciting fields coupled to the inflaton by rapidly and periodically
changing their masses [15]. Another efficient preheating mechanism is
tachyonic preheating, where the squared mass of a perturbed field becomes
negative, and perturbations are produced due to a tachyonic instability [16,
17].

Although tachyonic preheating is often studied in hilltop [16, 17, 18] or
hybrid inflation [19, 20, 21, 22, 23], or in other multi-field setups [24, 25,
26, 27], it is also generic to plateau inflation. In particular, the second
derivative of the inflaton potential, i.e., the effective mass squared of
perturbations, becomes negative at the plateau’s edge, causing fragmentation
of the homogeneous inflaton condensate111A similar effect, dubbed the flapping
resonance, has been studied in the context of axion physics [28, 29].. This
process is most efficient if the inflaton repeatedly returns to the plateau’s
edge during the oscillatory phase instead of being quickly damped by Hubble
friction. As shown in previous studies of Palatini Higgs and $R^{2}$ inflation
[30, 31, 32], this happens when the plateau is extremely flat, the tensor-to-
scalar ratio is small, $r\lesssim 10^{-7}$, and the energy scale of inflation
is relatively low, with a Hubble parameter $H\lesssim 10^{10}$ GeV. In this
regime, tachyonic preheating is, quite generally, the prevalent mechanism for
particle production. Fragmentation is very quick—the field can fragment
entirely in much less than an $e$-fold of expansion. This conclusion is backed
up by lattice computations [33, 34, 35, 36, 37]. A better understanding of the
preheating phase can improve the CMB predictions of these models, predict
gravitational wave (GW) signals, and inspire further model building. For
example, tachyonic preheating has been considered in the context of superheavy
dark matter [31].

In this paper, we study preheating in general models of inflation from an
exponentially flat plateau. Our primary focus is the tachyonic fragmentation
process, for which we outline the necessary conditions for fragmentation to be
efficient. To reach this goal, we apply linear perturbation theory and extend
it by including backreaction between perturbation growth and an adiabatically
evolving background. In this way, we can probe the boundaries of the
semianalytic linearized approach without resorting to full non-linear lattice
simulations. For the inflationary models considered here, we construct simple
(semi)analytic approximations for both the background oscillation dynamics and
the perturbation spectrum. We do not explicitly consider inflaton’s
perturbative decay into Standard Model particles.

We work out preheating in detail for inflaton potentials of the form
$U(\phi)=U_{0}\tanh^{n}\phi/\phi_{0}$ and derive the preheating duration, the
growth rate of the perturbations, and the wavenumber of the leading
perturbation in terms of the model parameters. Consistently with previous
studies [30, 31, 32], we find that preheating is very efficient. Our results
turn out to depend only weakly on the specifics of the model, which hints
towards the existence of universal features underlying preheating in plateau
inflation with exponentially flat potentials.

The paper is organized as follows. In section 2, we outline the general
features of plateau inflation and consider the evolution of the classical
background in the inflationary and the oscillatory epochs. Section 3 is
dedicated to preheating. We give a general theoretical overview of the
tachyonic instability of linear perturbations and numerically and analytically
estimate how the growth of perturbations backreacts on the classical
background. We conclude in section 4. Some technical details are gathered in
the appendices. We use natural units $\hbar=c=M_{\rm P}=1$ and the metric
signature $(-,\\!+,\\!+,\\!+)$.

##  2 Plateau inflation

We consider single-field inflation with a minimally coupled inflaton field
$\phi$ given by the action

$S=\int{\rm
d}^{4}x\sqrt{-g}\left[\frac{1}{2}R-\frac{1}{2}(\partial\phi)^{2}-U(\phi)\right]\,,$
(2.1)

where $g$ is the determinant of the metric $g_{\mu\nu}$, $R$ is the Ricci
scalar, and $U$ is the inflaton potential. In the spatially flat
Friedmann–Robertson–Walker Universe, expansion is controlled by the Friedmann
equation $3H^{2}=\rho$, where $H\equiv\dot{a}/a$ is the Hubble parameter, $a$
is the scale factor, and dot indicates a derivative with respect to time. The
energy density and pressure of the inflaton field are
$\rho=\dot{\phi}^{2}/2+U$ and $P=\dot{\phi}^{2}/2-U$, respectively, and the
field equation reads

$\ddot{\phi}+3H\dot{\phi}+U^{\prime}(\phi)=0\,.$ (2.2)

In the models considered here, the evolution of the background can be split
into two epochs: inflation and a subsequent phase in which the inflaton
oscillates around the minimum of its potential. We will consider these
separately in the following sections.

During the oscillatory phase, perturbations are generated from instabilities.
We are interested in tachyonic preheating, which is triggered by the tachyonic
instability. The existence of the latter sets requirements for the shape of
the potential:

  * •

The potential must have a region in which

$m_{\rm eff}^{2}\equiv U^{\prime\prime}(\phi)<0\,,$ (2.3)

thus making tachyonic instability possible. This condition is automatically
satisfied if the potential has a plateau curving down towards a minimum. In
that case, the tachyonic instability is strongest at the plateau’s edge.

  * •

The inflaton must repeatedly return to this unstable region, e.g. , to the
plateau. This is possible when Hubble friction is not too effective in damping
the amplitude of the post-inflationary oscillations. In particular, inflation
should end before reaching the edge of the plateau. This implies that the
period of oscillations must be shorter than the Hubble time, which can be
recast into the necessary condition

$|U^{\prime\prime}|\gg H^{2}$ (2.4)

somewhere in the oscillatory region. Furthermore, $U^{\prime\prime}$ is
relevant for another timescale – namely, it sets the timescale of the
tachyonic instability. This should also be shorter than the Hubble scale,
giving $-U^{\prime\prime}\gg H^{2}$ again, so that the tachyonic instability
can dominate over Hubble friction resulting in rapid mode growth. In our case,
(2.4) must hold at the edge of the plateau.

Although we aim to address tachyonic preheating in general terms, for the sake
of concreteness, we will work out numerical details for the potentials

$U=U_{0}\tanh^{n}(\phi/\phi_{0})\,,$ (2.5)

with $n$ an even integer. Such Einstein frame potentials arise naturally, e.g.
, in inflationary models in the Palatini formulation [38, 39, 40, 41, 42] and
the string theory motivated $T$-model $\alpha$-attractors [43, 44, 45, 46,
47]. In the context of tachyonic preheating, the cases $n=2,3,4,6$ have been
studied using linear perturbation theory [30, 32, 33, 34] for a wide range of
$\phi_{0}$ values, and using a non-linear lattice setup [33, 34, 36] with
$\phi_{0}\gtrsim 10^{-2}$ (corresponding to the edge of tachyonicity, see
below). In addition, production of GWs during preheating for $n=2$ and
$\phi_{0}\gtrsim 10^{-4}$ has been studied on a lattice [37]. Analytical
conditions for inflaton fragmentation in similar models were previously
studied in [48, 49]. We aim to generalize these studies by providing a
systematic understanding of the linear preheating regime and deriving scaling
relations applicable for a wide range of $\phi_{0}$ and $n$.

The potentials (2.5) have a plateau when $\phi\gtrsim\phi_{0}$ and the
condition (2.3) is satisfied on this plateau. For future convenience, we
define the mass scale

$m_{\rm th}\equiv\frac{\sqrt{U_{0}}}{\phi_{0}}\,,$ (2.6)

characterizing the magnitude of the effective mass of scalar perturbations,
$m^{2}_{\rm eff}=U^{\prime\prime}$. For example, at the edge of the plateau,
the effective squared mass of scalar fluctuations reaches values as low as
$U^{\prime\prime}_{\rm min}\approx-0.65m_{\rm th}^{2}$ almost independently of
$n$. Thus, $m_{\rm th}$ also controls the strength of the tachyonic
instability. The second condition (2.4) then implies that the tachyonic
instability can be active when

$\phi_{0}\ll 1\,,$ (2.7)

where we used that, close to the plateau, inflaton’s energy density will be
potential dominated, i.e., $H^{2}\approx U_{0}/3$, even in the oscillatory
regime.

Interestingly, as we will see below, many of our models’ features, from
inflationary phenomenology to the specifics of tachyonic preheating, can be
deduced from the large $\phi$ behaviour,

$U\sim U_{0}(1-A\exp(-2\phi/\phi_{0}))\,,\qquad\qquad\phi\gtrsim\phi_{0}\,,$
(2.8)

where $A=2n$ for the potentials (2.5). This can be understood intuitively by
noting that the oscillating field does not spend much time around the minimum
of the potential, so both inflationary dynamics and tachyonic preheating are
determined by features of the plateau. We expect the large-$\phi$ behaviour to
be a deciding factor also for other similar classes of potentials. Such models
are fairly common [50, 9] and include, e.g. , the well-known Starobinsky
inflation [1, 2] and the metric Higgs inflation [51, 52, 39], though for them
$\phi_{0}=\sqrt{6}$, which lies outside of the tachyonic parameter region.

###  2.1 Inflationary phase

Plateau inflation can be accurately studied with the slow-roll approximation.
Using the potential slow-roll parameters
$\epsilon_{U}\equiv\frac{1}{2}(U^{\prime}/U)^{2}$ and $\eta_{U}\equiv
U^{\prime\prime}/U$, we can write down the inflationary CMB observables: the
power spectrum of scalar perturbations, the spectral index, and the tensor-to-
scalar ratio,

$A_{s}=\frac{1}{24\pi^{2}}\frac{U}{\epsilon_{U}},\qquad
n_{s}=1-6\epsilon_{U}+2\eta_{U},\qquad r=16\epsilon_{U}\,.$ (2.9)

All quantities are computed when the comoving pivot scale
$k_{*}=0.05\text{Mpc}^{-1}$ exits the horizon.

For the class of potentials (2.5), the number of $e$-folds of inflationary
expansion corresponding to the field value $\phi$ is

$N\approx\int^{\phi}_{0}\frac{{\rm
d}\phi}{\sqrt{\epsilon_{U}}}=\frac{\phi_{0}^{2}}{2n}\sinh^{2}\left(\frac{\phi}{\phi_{0}}\right)\qquad\Rightarrow\qquad\phi\approx\frac{\phi_{0}}{2}\ln\frac{8nN}{\phi_{0}^{2}}\,,$
(2.10)

where the last identity holds when $N\gg 1$.222To estimate the number of
$e$-folds we begin integration from $\phi=0$. This introduces a negligible
error. Note also that plateau inflation can end violating the second slow-roll
condition $|\eta_{U}|\ll 1$ instead of the first one $|\epsilon_{U}|\ll 1$.
Thus, at the leading order in the large $N$ approximation,

$A_{s}\approx\frac{N^{2}m_{\rm th}^{2}}{3\pi^{2}},\qquad n_{s}\approx
1-\frac{2}{N},\qquad r\approx\frac{2\phi_{0}^{2}}{N^{2}}\,.$ (2.11)

Note that the dependence on $n$ has dropped out. For these expressions to
hold, we need $\phi_{0}^{2}\ll 2nN$, which, by Eq. (2.7), is satisfied
whenever we require tachyonic preheating to be active. The inflationary
predictions (2.11) are a generic feature of models in which the potential has
the exponential behaviour (2.8).

The number of $e$-folds $N_{*}$ corresponding to the CMB pivot scale depends
on the details of reheating [53]. Assuming instantaneous reheating, we have
$N_{*}\approx 61+\frac{1}{4}\ln(U_{0})$. In the following we will neglect the
$\ln(U_{0})$ dependence and use $N_{*}\approx 50$ yielding $n_{s}\approx 0.96$
compatible with the Planck measurement [8]. The observed power spectrum
strength $A_{s}=2.1\times 10^{-9}$ fixes the value of our mass parameter,

$m_{\rm th}=5\times 10^{-6}\,.$ (2.12)

The condition for tachyonic preheating (2.7) further implies that $r\ll
8\times 10^{-4}$, $U_{0}\ll 2.5\times 10^{-11}$.

###  2.2 Oscillatory phase

In the oscillatory phase, we closely follow the formalism of Ref. [32] and
work in the adiabatic limit in which the oscillation period is much shorter
than the Hubble scale. In this case, a single oscillation can be studied by
neglecting the expansion. The damping of energy density over multiple
oscillations, on the other hand, can be inferred from the continuity equation

$\dot{\rho}+3H(\rho+\bar{P})=0\,,$ (2.13)

by considering the time-averaged pressure $\bar{P}$ and an effective equation
of state $\bar{P}=\bar{P}(\rho)$ [54]. We remark that, by using (2.13), we
neglect the backreaction from perturbation growth. If the latter is extremely
fast, the coherent field may fragment before completing a single oscillation,
so an oscillatory regime is never realized. We will return to this issue in
section 3.3.

In the adiabatic approximation, we first consider the evolution of the
homogeneous field on a flat background. This can be described as a mechanical
system with the action $\int{\rm d}tP(\phi,\dot{\phi})$ and a conserved energy
density $\rho(\phi,\dot{\phi})=\dot{\phi}\partial P/\partial\dot{\phi}-P$. All
relevant quantities can be computed from the abbreviated action

$W(\rho)\equiv\int^{\phi_{2}}_{\phi_{1}}{\rm d}\phi\left.\frac{\partial
P}{\partial\dot{\phi}}\right|_{\dot{\phi}=\dot{\phi}(\phi,\rho)}=2\int^{\phi_{\text{amp}}}_{0}{\rm
d}\phi\sqrt{2(\rho-U(\phi))}\,,$ (2.14)

between the turning points $\phi_{i}(\rho)$ which are defined as values of the
field at which the velocity vanishes, $\dot{\phi}(\phi_{i},\rho)=0$. The
expression $\dot{\phi}(\phi,\rho)$ is obtained by inverting
$\rho=\dot{\phi}\partial P/\partial\dot{\phi}-P$. While the first expression
in (2.14) is general, the second expression applies in the case of a canonical
kinetic term, as in (2.1), and for symmetric potentials, i.e.
$U(\phi)=U(-\phi)$, which we assume throughout the paper unless specified
otherwise. Then $|\phi_{i}|\equiv\phi_{\text{amp}}$, the oscillation amplitude
of the coherent field component.

The half-period and the time-averaged pressure are [32]

$T=\partial_{\rho}W,\qquad\qquad\bar{P}=W/T-\rho\,.$ (2.15)

By eliminating the time-averaged pressure from the continuity equation (2.13),
we obtain a closed form equation for the energy density,

$\dot{\rho}+3HW/W^{\prime}(\rho)=0\,,$ (2.16)

where the second term describes the average effect of Hubble friction for an
oscillating field in a general potential. Note that the abbreviated action $W$
obeys the considerably simpler equation $\dot{W}+3HW=0$ and is thus conserved
in a comoving volume, that is, $W\propto a^{-3}$. The fraction of energy lost
during a half-oscillation is

$\frac{\Delta\rho}{\rho}=-\frac{W}{H}\,,$ (2.17)

given the scalar dominates the energy density so that $3H^{2}=\rho$.

Figure 1: _Left panel:_ The abbreviated action $W$ (solid) for the potential
(2.5) with $n=2,4,6,8,10$ together with its asymptotics at $\rho\sim 0$ (dot-
dashed) and $\rho\sim U_{0}$ (dashed). _Right panel:_ The dependence of the
half-period on the energy density.

Having these theoretical tools at our disposal, let us consider the potential
(2.5). The abbreviated action corresponding to (2.5) is given by

$W\approx\sqrt{8U_{0}}\phi_{0}f_{n}(\rho/U_{0})\,,\qquad f_{n}(x)\equiv
x^{\frac{1}{n}+\frac{1}{2}}\int^{1}_{0}\frac{\sqrt{1-y^{n}}}{1-x^{2/n}y^{2}}{\rm
d}y\,.$ (2.18)

Close to the plateau ($\rho\approx U_{0}$) and the minimum ($\rho\ll U_{0}$)
we find

$f_{n}(x)\sim\left\\{\begin{array}[]{lc}\frac{\sqrt{\pi}}{2n}x^{\frac{1}{n}+\frac{1}{2}}\left(\frac{\Gamma(1/n)}{\Gamma(3/2+1/n)}+\frac{\Gamma(3/n)/\Gamma(3/2+3/n)}{x^{-\frac{2}{n}}-\frac{\Gamma(3/2+3/n)\Gamma(5/n)}{\Gamma(3/2+5/n)\Gamma(3/n)}}\right)&,\mbox{
when }x\sim 0\\\ f_{n}(1)-\frac{\pi}{2}\sqrt{1-x}&,\mbox{ when }x\sim
1\,.\end{array}\right.$ (2.19)

The abbreviated action together with these asymptotics is shown in Fig. 1 for
selected values of $n$. When $n=2$, the $x\sim 1$ asymptotic matches
$f_{2}(x)$ exactly (see appendix C). For higher $n$, analytic expressions can
be obtained in terms of hypergeometric functions.333For even $n$,
$\displaystyle f_{n}(x)$
$\displaystyle=x^{\frac{1}{n}+\frac{1}{2}}\int^{1}_{0}\sqrt{1-y^{n}}\sum^{n/2-1}_{k=0}\sum^{\infty}_{m=0}(x^{\frac{2}{n}}y^{2})^{\frac{n}{2}m+k}{\rm
d}y$
$\displaystyle=\frac{\sqrt{\pi}}{2n}\sum^{n/2-1}_{k=0}x^{\frac{2k+1}{n}+\frac{1}{2}}\frac{\Gamma\left(\frac{1+2k}{n}\right)}{\Gamma\left(\frac{1+2k}{n}+\frac{3}{2}\right)}{}_{2}F_{1}\left(1,\frac{1+2k}{n},\frac{1+2k}{n}+\frac{3}{2};x\right).$
As a consistency check, we find that the leading term of the $x\sim 0$
asymptotic matches the corresponding result for $\phi^{n}$ potentials [54].

Consider now the first oscillations, as they tend to be the most relevant for
the dynamics of tachyonic preheating. For the tachyonic instability to become
active, the first oscillations must return to the plateau, so they take place
in the regime in which $\rho\approx U_{0}$. From (2.17),

$\frac{\Delta\rho}{\rho}\approx\phi_{0}2\sqrt{6}f_{n}(1)\,,$ (2.20)

so the energy loss per oscillation from Hubble friction is indeed small when
the condition (2.7) holds, i.e., when $\phi_{0}\ll 1$. At the first turning
point, $\rho=U_{0}-\Delta\rho=U(\phi_{\text{amp},i})$, so the initial
oscillation amplitude is

$\phi_{\text{amp},i}\approx-\frac{\phi_{0}}{2}\ln\left(\phi_{0}\frac{\sqrt{6}f_{n}(1)}{n}\right)\,,$
(2.21)

so the initial value of the field lies indeed on the plateau
($\phi_{\text{amp},i}\gtrsim\phi_{0}$) as long as $\phi_{0}$ is small. As a
consistency check, the slow-roll parameters at this point are
$\epsilon_{U}\approx 48f_{n}(1)$ and
$\eta_{U}\approx-8\sqrt{6}f_{n}(1)/\phi_{0}$. Both are large, so oscillations
start after the slow-roll epoch. Specifically, $\eta_{U}$ can be several
orders above one at that point.

By (2.15) and (2.19), the half-period of the first oscillations reads

$T\approx\frac{\pi\phi_{0}}{\sqrt{2(U_{0}-\rho)}}=\frac{\pi m_{\rm
th}^{-1}}{\sqrt{2(1-\rho/U_{0})}}\,.$ (2.22)

Half-periods for the next four even $n$ are shown in Fig. 1. The half-period
is proportional to the inverse mass scale $m_{\rm th}^{-1}$ but diverges for
large oscillation amplitudes when $\rho\to U_{0}$, so the condition (2.4) is
necessary but might not be sufficient to guarantee oscillations faster than
the Hubble rate. However, fixing the initial oscillation amplitude from (2.21)
shows that the initial half-period is

$HT_{i}\approx\frac{0.6\sqrt{\phi_{0}}}{\sqrt{f_{n}(1)}}\,,$ (2.23)

so the adiabaticity condition $T\ll 1/H$ still applies as long as $\phi_{0}\ll
1$, consistent with (2.7). Plugging in $H\sim\sqrt{U_{0}}$, we further see
that $T_{i}\sim(m_{\rm th}\sqrt{\phi_{0}})^{-1}$.

Since $W\propto a^{-3}$, Eq. (2.18) implies that for oscillations reaching the
plateau, $\rho\approx U_{0}$, the background’s energy density transitions from
vacuum energy domination as $\rho/U_{0}\propto 1-(a^{-3}-c_{n})^{2}$, where
$c_{n}$ is a numerical constant depending on $n$. A special case of such
behaviour was observed in Ref. [32]. For oscillations at the bottom of the
potential, $\rho\ll U_{0}$, we recover the well-known dilution $\rho\propto
a^{-6n/(2+n)}$ of $\phi^{n}$ potentials [54]. However, the latter scaling is
not realized in our models due to the fast fragmentation of the homogeneous
field.

Finally, eqs. (2.18) and (2.22) show that the half-period in the $\rho\approx
U_{0}$ region is independent of $n$ in the first approximation, and the
abbreviated action depends on $n$ only through a numerical constant (see Fig.
1 for specific values of $n$). This behaviour is a general consequence of the
exponential large field behaviour (2.8). Indeed, assuming a symmetric
potential with an asymptotic exponential tail of the form (2.8), a
straightforward computation shows the half-period is approximately (2.22)
since most of the time is spent on the plateau. Integrating Eq. (2.15) then
gives the general shape

$W\sim\pi\sqrt{2U_{0}}\phi_{0}(C-\sqrt{1-\rho/U_{0}}),\qquad\mbox{when}\qquad\rho\sim
U_{0}\,,$ (2.24)

where $C$ is a number that depends on the shape of the potential away from the
asymptotic region, e.g. , for the potential considered here with
$n=2,4,6,8,10$ the corresponding values read $C=1,1.22,1.34,1.44,1.51$,
respectively. Thus, the above results for the $U_{0}\approx\rho$ case also
apply for this more general class of potentials to which (2.5) belongs.

##  3 Preheating

While the homogeneous inflaton condensate oscillates, its perturbations start
to grow. This growth proceeds through two channels: parametric resonance [55,
14, 56, 15] and tachyonic instability [16, 17]. Both effects can lead to
preheating, i.e., the complete fragmentation of the coherent background field.
Parametric resonance takes place due to the time dependence of the frequency
of perturbations $\omega_{k}$ and can excite modes in frequency bands
corresponding to integer multipliers of a base frequency. We focus, instead,
on scenarios in which the tachyonic instability, triggered by a negative
squared frequency $\omega_{k}^{2}<0$, dominates preheating.

The scalar perturbations $\delta\phi_{k}$ with wavenumber $k$ follow, to
linear order, the equations of motion

$\delta\ddot{\phi}_{k}+3H\delta\dot{\phi}_{k}+\omega_{k}^{2}\delta\phi_{k}=0\,,\qquad\omega_{k}^{2}\equiv
k^{2}/a^{2}+U^{\prime\prime}(\phi)\,,$ (3.1)

together with the Bunch–Davies initial conditions,
$\delta\phi_{k}=1/(\sqrt{2k}a)$, $\delta\dot{\phi}_{k}=-i(k/a)\delta\phi_{k}$
[57] (see appendix A). The problem then contains two relevant timescales: the
half-period of the background oscillations $T$ sets the timescale for the time
variation of $\omega_{k}$, while $\omega_{k}$ itself sets the relevant
timescale for the growth of the perturbations $\delta\phi_{k}$. As we will
show, these two scales are of similar order.

We will next discuss perturbation growth and the solutions of (3.1) generally.
Focusing on our potentials (2.5), we will then provide simple analytical fits
for the fastest growing modes. Continuing the approach of section 2.2, we will
treat the effect of Hubble friction adiabatically. Assuming the field does not
fragment too rapidly, we will also estimate how perturbation growth backreacts
on the evolution of the background field. In particular, we compute the
duration of the preheating process.

###  3.1 Adiabatic mode growth

Figure 2:  _Left:_ Two oscillations of the background field in the potential
(2.5) with $n=2$ and an amplitude $\phi_{\text{amp}}=3\phi_{0}$ in the $H\to
0$ limit. _Middle:_ Evolution of $\omega_{k}$ for the peak mode with
$k=k_{\text{peak}}$ during the same period. _Right:_ Independent solutions for
the peak mode (red and blue) and the average exponential growth $e^{\mu_{k}t}$
(dashed grey).

In plateau inflation, the characteristic timescales for the perturbations are
much shorter than the Hubble time and the effect of expansion can be studied
perturbatively as we did for the classical background. Thus, at the leading
order, mode growth is computed in the $H\to 0$ limit. In this case, Eq. (3.1)
simplifies to the Hill equation

$\delta\ddot{\phi}_{k}+\omega_{k}^{2}\delta\phi_{k}\approx 0\,,$ (3.2)

where the period of $\omega_{k}^{2}$ matches the background’s half-period
$T$.444This holds for symmetric potentials assumed here. If the potential is
asymmetric, then the period of $\omega_{k}^{2}$ coincides with the full period
of background oscillations instead of the half-period. This is illustrated in
Fig. 2. Floquet’s theorem then invites us to seek solutions that are periodic
up to a growth factor [58], i.e.,

$\delta\phi_{k}(t+T)=e^{\lambda_{k}T}\,\delta\phi_{k}(t)\,,$ (3.3)

where $\lambda_{k}$ is a Floquet exponent. Since (3.2) is a second order
equation, there will be two Floquet exponents. The growth rate of the mode is
given by the largest real part of these exponents,

$\mu_{k}\equiv\max\real\lambda_{k}.$ (3.4)

In this subsection, we will focus on the computation of this growth rate.

Given any two independent solutions $u_{1}(t)$ and $u_{2}(t)$ of the $H\to 0$
mode equation (3.2), their Wronskian matrix $w(t)$ evolves as

$w(t+T)=w(t)G\,,\qquad w(t)\equiv\begin{pmatrix}u_{1}(t)&u_{2}(t)\\\
\dot{u}_{1}(t)&\dot{u}_{2}(t)\\\ \end{pmatrix},$ (3.5)

where $G\equiv w(0)^{-1}w(T)$ is a constant555One can show that ${\rm
d}(w(t)^{-1}w(t+T))/{\rm
d}t=w(t)^{-1}(\dot{w}(t+T)-\dot{w}(t)w(t)^{-1}w(t+T))=0$ by noting that the
mode equation (3.2) implies that the Wronskian matrix satisfies $\dot{w}=Dw$
where $D(t+T)=\pm D(t)$ is a periodic matrix. matrix dubbed the monodromy
matrix. Eigenfunctions of $G$ satsify (3.3) and thus $e^{\lambda_{k}T}$ must
be an eigenvalue of $G$. The most convenient choice of $u_{1,2}$ depends on
the context.666For different choices of the two independent solutions
$u_{1,2}$, the corresponding $G$ matrices are related by a similarity
transformation $G\to LGL^{-1}$, thus the Floquet exponents are independent of
the choice of $u_{1,2}$. For example, when looking for the Floquet exponent
numerically, we solve (3.3) with initial conditions $u_{1}=\dot{u}_{2}=1$,
$u_{2}=\dot{u}_{1}=0$, or equivalently, $w(0)=1$, so that $G=w(T)$.

Figure 3: _Upper panels:_ Growth rate $\mu_{k}$ of modes in the $H\to 0$ limit
for the potential (2.5) with $n=2,4,6,8,10$ assuming a background with
amplitude $\phi_{\rm amp}=2\phi_{0}$ (left panels) and $\phi_{\rm
amp}=4\phi_{0}$ (right panels). _Lower panels:_ ${\rm tr}\,G/2$ corresponding
to the panel above it. The grey band shows the region in which the modes are
stable, i.e., $\mu_{k}=0$.

The equations of motion imply that the Wronskian $\det w(t)$ is constant in
time, so taking the determinant of both sides of (3.5) we find that $\det
G=1$. Therefore, $\lambda_{k,-}=-\lambda_{k,+}$, and, since $G$ is a real
matrix, $\lambda_{k\pm}$ are either real or purely imaginary. Both Floquet
exponents are given by the simple expression

$\lambda_{k,\pm}=\pm T^{-1}{\rm acosh}\left(\frac{1}{2}{\rm tr}\,G\right).$
(3.6)

We drop the $\pm$ subindices and denote by $\lambda_{k}$ the eigenvalue with
positive real part $\mu_{k}=\real\lambda_{k}\geq 0$. The Floquet exponents are
real (imaginary) if $|{\rm tr}\,G|>2$ ($|{\rm tr}\,G|<2$). If $\lambda_{k}$ is
real, one of the solutions is growing and quickly starts to dominate over the
other solution, which is decaying. Such growth marks an instability of the
mode.

An example of mode growth in the $H\to 0$ limit is shown in Fig. 2. Through
$\omega_{k}^{2}$, the growth rate $\mu_{k}$ depends on the shape of the
potential and the oscillation amplitude, in addition to the wavenumber $k$.
Fig. 2 uses the potential (2.5) with $n=2$ where the background field,
depicted in the left panel, oscillates with an amplitude
$\phi_{\text{amp}}=3\phi_{0}$ and a half-period $T=22.4m_{\rm th}^{-1}$. The
right panel of Fig. 2 shows the evolution of two independent solutions to the
mode equation (3.2) for the fastest growing mode $k_{\text{peak}}/a\approx
3.54/T$. Both of these solutions contain the growing mode, so they look
visually identical after a first half-period. In the middle panel, the squared
frequency of the mode $\omega_{k}^{2}$ is seen to take mostly negative values
except when the background field briefly crosses the origin. Thus, most of the
time, the mode follows well the simple exponential curve $e^{\mu_{k}t}$, shown
by a dashed grey line in the left panel. Such exponential growth would be
exact when $\omega_{k}^{2}$ is constant and negative.

Numerical estimates of the growth rate $\mu_{k}$ over a range of $k$-values
are shown in Fig. 3, together with the corresponding trace of the monodromy
matrix (Hill’s discriminant), for potentials (2.5) with $n=2,4,6,8,10$ and a
fixed background amplitude $\phi_{\text{amp}}$.777Floquet charts for
$n=2,3,4,6$ have been computed also in Ref. [33, 34]. As explained in appendix
A, dimensionless quantities such as $\mu_{k}T$ do not depend explicitly on
dimensionful model parameters such as $\phi_{0}$ and $U_{0}$. Note that fixing
$\phi_{\text{amp}}/\phi_{0}$ means that both $\rho$ and $T$ vary with $n$.
With this convention, the first peak only varies slightly for different $n$
when the amplitude is large, $\phi_{\text{amp}}=4\phi_{0}$. The differences
become more visible when $\phi_{\text{amp}}=2\phi_{0}$, although the order of
magnitude remains unchanged. The reduction of differences between different
$n$ with increasing amplitudes is expected as for larger $\phi_{\text{amp}}$
the background dynamics becomes more and more dominated by the $n$-independent
exponential tail. In all, when $\phi_{\text{amp}}\gg\phi_{0}$, the first and
dominant tachyonic peak has a nearly universal shape.

Figure 4: The dependence of maximal value of $\mu_{k}$ (left panel) and the
wavenumber of the corresponding mode $k_{\text{peak}}$ (right panel) on the
oscillation amplitude $\phi_{\text{amp}}$ for the potential (2.5) with
$n=2,4,6,8,10$. The dashed lines correspond to cases where the $\mu_{k}$ does
not peak in the first instability band but in one the secondary bands instead.
The grey dotted lines show the asymptotic behaviour at large amplitudes:
$\mu_{\text{peak}}T-\phi_{\text{amp}}/\phi_{0}\in(1/2,1)$ in the left panel
and $k_{\text{peak}}T=3.54$ in the right panel.

The same cannot be said for higher $k$ peaks—even when the first peaks are
matched. The higher peaks can vary significantly for different $n$ (even for
$\phi_{\text{amp}}=4\phi_{0}$). However, as can be observed from Fig. 3, the
instability bands vary nearly periodically in $k$-space. For
$\phi_{\text{amp}}=4\phi_{0}$, the end of the instability bands is
approximately $k/a\approx j\pi/T$, with $j\geq 2$ an integer. The secondary
instability bands arise predominantly due to parametric resonance rather than
the tachyonic instability as $\omega_{k}^{2}$ is mostly positive in these
bands. They are subleading compared to the dominant tachyonic peak.

The growth rate and wavenumber of the fastest growing mode are shown in Fig.
4. For sufficiently large amplitudes they can be approximated by

$\mu_{\text{peak}}\approx\frac{\phi_{\text{amp}}/\phi_{0}+c}{T}\,,\qquad
k_{\text{peak}}/a\approx\frac{3.54}{T}\qquad\qquad(\phi_{\text{amp}}\gtrsim
3\phi_{0})\,,$ (3.7)

where $c$ is an $\mathcal{O}(1)$ coefficient that lies in the range (0.5,1)
for the cases considered in Fig. 4.888In detail, for $n=2,4,6,8,10$ we find
$c=0.93,0.79,0.68,0.6,0.53$, respectively. Since $\phi_{0}$ is related to the
maximal oscillation amplitude by Eq. (2.21), $\phi_{0}$ must be relatively
small for these approximations to be valid. The condition
$\phi_{\text{amp},i}/\phi_{0}\gtrsim 3$ translates into $\phi_{0}\lesssim
10^{-3}$. For $\phi_{0}\approx 10^{-2}$, $\phi_{\text{amp},i}/\phi_{0}$ goes
from 2 to 2.6 when $n$ varies from 2 to 10, so the above approximation barely
holds even during the initial oscillations, although it still gives the
correct order of magnitude for $\mu_{\text{peak}}$ and $k_{\text{peak}}$. For
even larger $\phi_{0}$ and thus for lower initial amplitudes
$\phi_{\text{amp},i}$, the instability within the first band weakens and the
secondary peaks begin to dominate indicating that tachyonicity becomes less
relevant for preheating.

Relying on (3.7), we can derive simple order-of-magnitude estimates in terms
of the model parameters. $T_{i}\sim m_{\rm th}^{-1}\phi_{0}^{-1/2}$ gives
$k_{\text{peak}}\sim\mu_{\text{peak}}\sim m_{\rm th}\sqrt{\phi_{0}}$. This is
also the order of magnitude of
$\sqrt{-U^{\prime\prime}(\phi_{\text{amp},i})}$, the scale of tachyonicity, so
all relevant time and energy scales are of similar order.

We remark that the half-period $T$ grows rapidly when
$\phi_{\text{amp}}\to\infty$. Thus, by Eq. (3.7), $\mu_{\text{peak}}$ will
generally decrease in this limit. Generally, $\mu_{\text{peak}}$ will peak
when $\phi_{\text{amp}}\approx\mathcal{O}(1)$; the exact value depends on $n$.
When the amplitude decreases below this value as the field relaxes to the
minimum, $\mu_{\text{peak}}$ will approach 0.

According to the numerical estimates in Figs. 3 and 4,
$\mu_{\text{peak}}\gtrsim T^{-1}$ in the tachyonic regime, that is, the
fastest growing mode can increase significantly during a half-oscillation.
Therefore, after these modes begin to contribute noticeably to the total
energy density, we expect that strong feedback effects will lead to a fast
breakdown of adiabatic background evolution. Moreover, since the instability
gets stronger with decreasing $\phi_{0}$, the field may become significantly
fragmented even before the first oscillation is complete, and thus, it is
expected that the adiabatic approach breaks down for very small $\phi_{0}$.
This will be shown in section 3.3.

Finally, consider the asymptotics of $\mu_{k}$. Firstly, the $k=0$ mode
$\delta\phi_{0}$ does not grow, e.g. , in Fig. 3, we see that $\mu_{k}\to 0$
as $k\to 0$. Specifically, the $k\to 0$ asymptotic for symmetric999In case of
asymmetric potentials, one must replace $T$ with the full period, $W$ with the
corresponding integral over the full period and flip the sign of the monodromy
matrix, i.e., ${\rm tr}\,G\to-{\rm tr}\,G$. potentials is (for a derivation,
see appendix B)

$\lambda_{k}=-\frac{i\pi}{T}+(k/a)\frac{\sqrt{W\partial_{\rho}T}}{T}+\mathcal{O}(k/a)^{3}\quad\Leftrightarrow\quad\frac{1}{2}{\rm
tr}\,G=-1-\frac{(k/a)^{2}}{2}W\partial_{\rho}T+\mathcal{O}(k/a)^{4}\,.$ (3.8)

Thus the position of the first stability band is tied to whether the period of
background oscillations decreases with decreasing background energy density or
not: as long as $\partial_{\rho}T>0$, the first instability band begins at
$k=0$ since ${\rm tr}\,G<-2$ for modes surrounding $k=0$. In particular, the
slope of $T$ shown in Fig. (1) implies that the first instability band starts
at $k=0$ for arbitrarily small background energy densities only when $n=2$ –
for larger $n$, however, the infrared modes are stable when the oscillations
are sufficiently damped. In general, for exponentially flat potentials (2.8)
and sufficiently large oscillation amplitudes, $W$ can be approximated by Eq.
(2.24) so that (3.7) gives

$\mu_{k}=(k/a)\sqrt{\frac{\sqrt{2}C}{\pi}m_{\rm
th}T-1}+\mathcal{O}(k/a)^{3}\,.$ (3.9)

In the special case of the potential (2.5) with $n=2$ this expression with
$C=1$ holds for any amplitude (see appendix C).

For large $k$, ${\rm tr}\,G/2\to\cos(kT/a)$ (see appendix B), as can be
observed from the lower panels in Fig. 3. It implies that $\mu_{k}\to 0$ when
$k\to\infty$, as expected. However, the way $\mu_{k}$ approaches 0 depends
non-trivially on the behaviour of background oscillations. General arguments
show that it must approach zero at least as $k^{-1}$ [58]. However, our
numerical estimates for the class of potentials (2.8) show that $|{\rm
tr}\,G|-2$ is suppressed exponentially as $k$ grows. As a result, the
instability bands get exponentially narrower and weaker when $k$ approaches
infinity.

###  3.2 Energy density of perturbations and backreaction

The energy density of perturbations can be estimated as the quantum average of
the corresponding energy density operator (see appendix A)

$\displaystyle\delta\rho(t)=\int\frac{{\rm
d}^{3}k}{(2\pi)^{3}}\delta\rho_{k}(t)\,,\qquad\mbox{where}\qquad\delta\rho_{k}\equiv\frac{1}{2}|\delta\dot{\phi}_{k}|^{2}+\frac{1}{2}\omega_{k}^{2}|\delta\phi_{k}|{}^{2}$
(3.10)

is the energy density of a single mode. Assuming adiabaticity and neglecting
expansion, the time averaged energy density of the perturbations grows as

$\delta\dot{\rho}_{k}(t)=2\mu_{k}(\rho)\delta\rho_{k}(t)$ (3.11)

by Eq. (3.3). As discussed above, mode growth depends on the potential and the
oscillation amplitude $\phi_{\text{amp}}$, or equivalently, the background
energy density $\rho$.101010Note that $\rho$ is reserved for the background
energy density throughout the paper.

We assume initial conditions corresponding to the vacuum state, which we take
to be of the Bunch–Davies form, $\delta\rho_{k}^{\rm vac}=\delta\rho_{k}^{\rm
BD}=ka^{-4}/2$.111111We neglected the contribution from $U^{\prime\prime}$
which induces at most an order one correction. In the numerical approach, we
set $a=1$ at the start of preheating.  As is well known, the energy density is
divergent in the ultraviolet, so the integral (3.10) must be regularized. In
the following, we will use a sharp or an exponential ultraviolet cut-off. In
particular, a sharp cut-off is implicit in our numerical approach below as the
mode equations are solved for a finite set of momenta. Furthermore, we must
renormalize the energy density of perturbations. We do this by demanding that
it vanishes in the initial vacuum state. The physical, renormalized energy
density $\delta\rho^{(R)}$ of a single mode then reads

$\delta\rho^{(R)}_{k}=\delta\rho_{k}-\delta\rho^{\rm vac}_{k}\,,$ (3.12)

where $\delta\rho^{\rm vac}_{k}=\delta\rho_{k}(0)$ is the initial regularized
vacuum energy density for the mode. Essentially, Eq. (3.12) renormalizes the
cosmological constant to its observed value, which is negligible during
preheating.

Above we considered mode growth neglecting the evolution of the background and
the expansion of space. To remedy this, we impose energy conservation for the
combined system of perturbations and background, so that, by Eq. (3.11),
energy is transferred from the background to the perturbations with the rate

$s(t)\equiv\int\frac{{\rm d}^{3}k}{(2\pi)^{3}}2\mu_{k}\delta\rho_{k}(t)\,.$
(3.13)

It matches $\delta\dot{\rho}^{(R)}$ when expansion is negligible.

It is instructive to split the total growth rate into two terms

$s(t)=\bar{\mu}(\rho)\delta\rho^{(R)}(t)+j(\rho)\,,$ (3.14)

where we defined

$\bar{\mu}(\rho)\equiv\int\frac{{\rm
d}^{3}k}{(2\pi)^{3}}2\mu_{k}(\rho)\frac{\delta\rho^{(R)}_{k}(t)}{\delta\rho^{(R)}(t)}\,,\qquad
j(\rho)\equiv\int\frac{{\rm d}^{3}k}{(2\pi)^{3}}2\mu_{k}(\rho)\delta\rho^{\rm
vac}_{k}\,.$ (3.15)

In this formulation, the $j$ term acts as an explicit quantum source for the
perturbations121212The integral defining (3.15) converges due to
$\mu_{k}(\rho)$ being exponentially damped at large $k$., originating from
vacuum fluctuations, and dominates initially when the energy density of
perturbations is negligible. Later, the effective growth rate $\bar{\mu}$
begins to drive the exponential growth of $\delta\rho^{(R)}$. Note that, on
top of $\rho$, $\bar{\mu}(\rho)$ and $j(\rho)$ depend on the spectrum of
perturbations. Therefore, one must still solve the evolution mode by mode.
However, in the idealized limit of a constant background density, one expects
that

$\lim_{t\to\infty}\bar{\mu}=2\mu_{\text{peak}},$ (3.16)

since the $k_{\text{peak}}$-mode grows the fastest and will eventually
dominate the energy density of perturbations. We will return to this limit in
section 3.2.1.

On top of the energy transfer from the background to the fluctuations, both
$\rho$ and $\delta\rho$ are diluted by Hubble friction; for $\rho$, the
dilution follows from (2.16), while $\delta\rho$ dilutes in a $k$-dependent
way. We can then describe the time-averaged evolution of all components of the
scalar’s energy density as

$\displaystyle\dot{\rho}+3HW(\rho)/W^{\prime}(\rho)$ $\displaystyle=-s(t)\,,$
(3.17) $\displaystyle\delta\dot{\rho}_{k}+c_{k}H\delta\rho_{k}$
$\displaystyle=2\mu_{k}(\rho)\delta\rho_{k}\,,$

with expansion determined from the physical energy density by
$3H^{2}=\rho+\delta\rho^{(R)}.$ The factors $c_{k}$ capture the dilution of
the $k$-modes. For massive modes, $c_{k}=3$, while for effectively massless
high-$k$ modes, $c_{k}=4$; we expect $c_{k}$ to vary between these values. To
capture the exact behaviour, $c_{k}$ would need to be solved from the full
mode equation (3.1). In the numerical approach below, we are content with
adequately describing the leading modes, for which the mass term
$U^{\prime\prime}$ is important, and thus we choose $c_{k}=3$ for all
$k$.131313Lattice studies in the parametric resonance regime show that the
fragmented field behaves as radiation when $n>2$ [33, 34] indicating that
$c_{k}=4$ would better capture the dilution after fragmentation. For
consistency, we take the vacuum modes $\delta\rho_{k}^{\rm vac}$ in (3.12) to
dilute the same way, so that they correspond to a non-growing solution with
$\mu_{k}=0$.

As long as the adiabatic approximation is valid, Eqs. (3.17) can be employed
to study more extended time periods, during which Hubble friction can be
relevant. However, as we will see below, in the parameter region allowing for
tachyonic preheating, the feedback from perturbation growth is much stronger
than Hubble friction, i.e., $\bar{\mu}\gg H$, and thus the latter plays an
insignificant role.

All of the above relies on linear perturbation theory and will change when
non-linearities become relevant. Interactions between the highly exited modes
will eventually lead to the thermalization of the fragmented field. As this
process redistributes the energy carried by the modes of the fragmented field,
it does not strongly affect the total energy density in the fragmented field.
However, as it can move modes in or out of instability bands, it could have a
sizable effect on the fragmentation process itself. In particular, even if the
power spectrum develops a dominant peak centred around $k_{\text{peak}}$ as
predicted by (3.16), it will eventually broaden due to non-linear effects.
This process, including the generation of secondary peaks due to rescattering,
has been observed in lattice studies [33, 34].

####  3.2.1 The dominant peak approximation

Time evolution simplifies if mode growth is dominated by a single peak in
$\mu_{k}$, as is the case for strongly tachyonic preheating. Around the peak
mode $k_{\text{peak}}$, we can use the quadratic expansion

$\mu_{k}\approx\mu_{\text{peak}}-\frac{\Delta}{2a^{2}}(k-k_{\text{peak}})^{2}\,$
(3.18)

to describe the first instability band and neglect mode growth in other
instability bands.141414The growth rate is positive as long as
$|k-k_{\text{peak}}|/a<\mu_{\text{peak}}\sqrt{2/\Delta}$. Applying this
approximation away from this range leads to unphysical damping of the modes.
However, as long as $k_{\text{peak}}$ does not evolve too fast, these far-off
modes remain subdominant and can thus be neglected due to the exponential
nature of perturbation growth. In this case, extending the approximation
(3.18) also to the far-off modes does not alter the estimate of the total
energy density of perturbations.

Consider first the simplified case where we neglect both expansion and
backreaction. Then $\mu_{k}$ is time independent and

$\delta\rho_{k}=\delta\rho_{k}^{\rm BD}e^{2\mu_{k}t}\,,$ (3.19)

where $\delta\rho_{k}^{\rm BD}$ is the initial vacuum energy. The quadratic
approximation (3.18) thus produces a Gaussian peak in the perturbation
spectrum with the width $(\Delta t)^{-1/2}$ that decreases in time. Relying on
(3.8), we can use $\mu_{0}=0$ to estimate $\Delta\approx
2\mu_{\text{peak}}a^{2}/k_{\text{peak}}^{2}$. In the limit
$2\mu_{\text{peak}}t\gg 1$, the total energy density of perturbations (3.10)
is then approximately151515Note that in the limit of large perturbations,
$\delta\rho_{k}\gg\delta\rho_{k}^{\mathrm{BD}}$ and thus
$\delta\rho_{k}^{(R)}\approx\delta\rho_{k}$.

$\delta\rho^{(R)}\approx\frac{(k_{\text{peak}}/a)^{4}}{4\pi^{3/2}\sqrt{2\mu_{\text{peak}}t}}e^{2\mu_{\text{peak}}t}.$
(3.20)

In this limit, $\mu_{\text{peak}}$ clearly determines the growth rate as
conjectured in (3.16). The time $\tau$ for the field to half-fragment, i.e.,
to reach $\delta\rho^{(R)}(\tau)\equiv U_{0}/2$ is

$\tau\approx\frac{1}{\mu_{\text{peak}}}\left[\ln\left(\frac{4.0\sqrt{U_{0}}}{(k_{\text{peak}}/a)^{2}}\right)+\frac{1}{4}\ln\ln\left(\frac{4.0\sqrt{U_{0}}}{(k_{\text{peak}}/a)^{2}}\right)\right]\,.$
(3.21)

Fixing $m_{\rm th}$ (2.12), the initial amplitude (2.21) and the corresponding
half-period of oscillations (2.23) and using the large amplitude estimate
(3.7) for $\mu_{\text{peak}}$ and $k_{\text{peak}}$, we obtain that the field
will half-fragment within

$N\approx\tau H\approx\mathcal{O}(1)\sqrt{\phi_{0}}$ (3.22)

$e$-folds. The $\mathcal{O}(1)$ factor depends on $n$ and contains
$\ln\phi_{0}$ corrections, e.g. , for $n=2$ and $\phi_{0}=10^{-4}$ it is 1.1.
These specifics are accounted below in Fig. 6, where this simple approximation
is compared with more accurate estimates for half-fragmentation.

The number of half-oscillations before half-fragmentation is small. Following
the procedure above, it is approximately

$\tau/T\approx\frac{10}{-\log_{10}(\phi_{0})}\,.$ (3.23)

It follows that the field can become significantly fragmented already after
the first oscillation when $\phi_{0}\ll 10^{-5}$. We note, however, that
Floquet theory gives the growth rate after a single half-oscillation and can
thus be applied even if the number of half-oscillations is small (given the
backreaction effects stay irrelevant).

To improve the estimate, let us include the possibility of a time dependent
background. The approximation (3.18) implies that the perturbation spectrum
evolves as a Gaussian peak of the form

$\displaystyle\delta\rho_{k}=\delta\rho_{k}^{\rm BD}\,e^{A-B(k-k_{0})^{2}},$
(3.24)

and the ($H\to 0$) evolution of $\delta\rho$ is given by

$\dot{k}_{0}=-\frac{\Delta}{a^{2}B}(k_{0}-k_{\text{peak}})\,,\quad\dot{B}=\frac{\Delta}{a^{2}}\,,\quad\dot{A}=2\mu_{\text{peak}}-\frac{\Delta}{a^{2}}(k_{0}-k_{\text{peak}})^{2}=2\mu_{k_{0}}\,,$
(3.25)

where $k_{\text{peak}}$, $\mu_{\text{peak}}$ and $\Delta$ depend on the
background energy density $\rho$ and thus on time. The time evolution drives
the centre of the distribution $k_{0}$ towards the fastest-growing mode
$k_{\text{peak}}$. As in the time-independent case, the width $B^{-1/2}$
decreases monotonously—the fastest-growing modes outpace all others and the
distribution gets narrower and narrower as time passes. The parameter $A$
captures the growth of the $k_{0}$ mode. Since the perturbations should start
from the vacuum with $\delta\rho_{k}=\delta\rho_{k}^{\rm BD}$, we choose the
initial conditions $A(0)=0$, $k_{0}(0)=k_{\text{peak},i}$ (the initial
fastest-growing mode), and $B(0)$ sufficiently small (the precise choice for
$B(0)$ is discussed below).

In the limit $B\gg 1/k_{0}^{2}$, we can replace $0\to-\infty$ in the momentum
integral (3.10), so the total perturbation energy density is

$\delta\rho\approx\frac{k_{0}}{4(\pi
B)^{3/2}}\quantity(\frac{3}{2}+Bk_{0}^{2})e^{A}\,.$ (3.26)

Since $B$ increases fast, the induced error is never large in practice. From
Eq. (3.25) we then obtain

$\displaystyle\delta\dot{\rho}$ $\displaystyle\equiv
s(t)\approx\delta\rho\times$ (3.27)
$\displaystyle\Bigg{(}2\mu_{\text{peak}}-\frac{\Delta}{a^{2}}\quantity[\frac{9+2Bk_{0}^{2}}{2B(3+2Bk_{0}^{2})}+\frac{3+6Bk_{0}^{2}}{Bk_{0}(3+2Bk_{0}^{2})}(k_{0}-k_{\text{peak}})+(k_{0}-k_{\text{peak}})^{2}]\Bigg{)}\,.$

The first term in the brackets arises from the fastest-growing mode and the
second term corrects for the finite (and changing) width of the distribution,
and we identified the result with the source term $s(t)$ from (3.13).

With these results, we can numerically integrate the coupled evolution of the
background and perturbations (3.17) by using (3.27) for $s(t)$, replacing the
$k$-equations with (3.25), and using $\delta\dot{\rho}+3H\delta\rho=s(t)$ to
include Hubble friction for the perturbations. We still use
$\delta\rho^{(R)}=\delta\rho-\delta\rho^{\mathrm{vac}}$, where
$\delta\rho^{\mathrm{vac}}=\delta\rho(0)$, to renormalize the perturbation
energy density for the Friedmann equation (note that in this scheme, for a
positive $B$, all the integrated quantities are finite). These equations must
still be solved numerically, but we have reduced the number of equations from
one for each $k$-mode to a total of five, with the cost of some accuracy.

We remark that, when we apply (3.27) at earlier times when the perturbation
spectrum is not narrow, we are effectively using the approximation (3.18) in
the region $|k-k_{\text{peak}}|>2\mu_{\text{peak}}a^{2}/\Delta$ where
$\mu_{k}<0$. If these modes would ever dominate, even $s(t)<0$ can be possible
due to the unphysical damping of vacuum energy at high $k$. Thus, the
quadratic approximation (3.18) can only apply if these modes are never
relevant during the time evolution of the system. In particular, the initial
$B$ should be sufficiently large so the leading modes in the Gaussian peak
(3.24) all reside inside the $\mu_{k}>0$ region; this also guarantees the
validity of (3.26) and (3.27). On the other hand, the initial vacuum state
should have a vanishing $B$. As a compromise, we choose the initial $B$ to be
the minimal value that still produces $s(0)\geq 0$. The initial $\mu_{k}<0$
modes must also stay unimportant at later times; to achieve this, the
background oscillation amplitude, and consequently $k_{\text{peak}}$ and
$\Delta$, should only change a little during the bulk of preheating. This
approximation is then only a slight improvement over the simplified case with
no background evolution considered at the beginning of this section but may be
useful close to the parameter region where the simplified approximation starts
to fail.

###  3.3 Results

We now solve the evolution of the background and the perturbations using the
tools discussed above to compute the duration of preheating and the ensuing
perturbation spectrum. We perform the computation in multiple ways, taking
varying amounts of computational resources. In all cases, we assume adiabatic
background evolution so that the background’s energy density $\rho$ decreases
only a little within a single oscillation. In the order of increasing
simplicity and decreasing accuracy, our approaches are:

  * •

Grid: We numerically integrate (3.17) over a grid of $k$-values. To facilitate
the computation, we first build the function $\mu_{k}(\rho)$ by numerically
solving $\mu_{k}$ as described in section 3.1 for several values of $\rho$ and
interpolating in-between, for all $k$-values on the grid. We use a logarithmic
grid with steps of $\Delta\ln k=0.05$ and covering the range $k=10^{-15}m_{\rm
th}\dots 10^{5}m_{\rm th}$, chosen to capture all important modes that have a
sizeable $\mu_{k}$ and may get excited during tachyonic preheating.161616In
practice, $\mu_{k}=0$ within the accuracy of the computation for all modes
higher than the cut-off value on the grid. All higher $k$ modes stay in their
vacuum state throughout the simulation, and the results are not sensitive to
the cut-off. We tested the convergence of the results by repeating chosen
computations with a doubled grid spacing, which produced only a negligible
difference.

  * •

Dominant peak approximation (DP): We compute the evolution of perturbations
with a time-dependent background as described in section 3.2.1 by solving Eq.
(3.27) with Hubble friction and backreaction included. In this setup,
$\mu_{k}$ is given by the quadratic approximation (3.18). The coefficients
$\mu_{\text{peak}}(\rho)$, $k_{\text{peak}}(\rho)$ and $\Delta(\rho)$ are
obtained by fitting the first tachyonic peak. The contributions of the
secondary peaks are neglected. As above, we build these functions by first
finding $\mu_{k}(\rho)$ them for different values of $\rho$, e.g. , as in Fig.
4.

  * •

Simple approach: We use the simplified version of the dominant peak
approximation which relies on Eqs. (3.20-3.22) and the fits (3.7) to compute
$N$ and the leading $k$. In this approximation, we assume a constant
background energy density, that is, we neglect both backreaction and Hubble
friction. Numerical integration is not needed in this approach.

In the regime of tachyonic preheating, preheating finishes in much less than
one $e$-fold (see Fig. 6). As a consequence, Hubble friction and redshift
effects play a subleading role. Time evolution of $\rho$ is then mainly driven
by the backreaction effects, as anticipated above. We neglect the redshift of
the wavenumber $k$ in $\mu_{k}$ to simplify all the computations.

Figure 5: Evolution of the background energy density $\rho$ and the final
energy spectrum of perturbations for $n=2$, for two example cases
$\phi=10^{-5}$ (left) and $\phi=10^{-3}$ (right). For $\phi_{0}=10^{-3}$, the
dominant peak approximation is also presented. The time evolution is shown in
terms of the number of $e$-folds $N$ and the number of background half-
oscillations $n_{\mathrm{osc}}$. Since the energy density of perturbations
grows exponentially, the backreaction-driven time evolution is slow at first
but accelerates towards the end. For the spectrum, the dotted line correspond
to the Hubble horizon and the dashed lines to the value of $k_{\text{peak}}$
initially and in the end. The corresponding frequencies today are also given.
The spectrum is given in terms of
$\Delta_{k}\equiv\frac{k^{3}}{2\pi^{2}}\delta\rho_{k}$, so that
$\delta\rho=\int_{0}^{\infty}{\rm d}\ln k\,\Delta_{k}$.

Consider the initial and final configurations of the system. After fixing $n$
and setting $m_{\rm th}$ from (2.12), the potential (2.5) contains a single
free parameter $\phi_{0}$. This sets the initial oscillation amplitude
$\phi_{\text{amp},i}$ and the corresponding background energy density $\rho$
through Eq. (2.21). Starting from the vacuum state with $\delta\rho^{(R)}=0$,
we follow the time evolution of the system until half-fragmentation
$\delta\rho^{(R)}=\rho$, that is, until the perturbations have drained half of
the background energy and preheating has essentially completed. Around this
point, our linear analysis breaks down, and subsequent time evolution needs to
be analysed by lattice methods. If Hubble friction is negligible, the final
value of the background oscillation amplitude is
$\phi_{\text{amp}}\approx\phi_{0}$ (with some $n$-dependence), solved from
$U(\phi_{\text{amp}})=U_{0}/2$. At this point, we end the simulation and
record the number of $e$-folds $N$ and the perturbation spectrum. We scan over
a range of $\phi_{0}$ values and repeat the analysis for $n=2$, $4$, and $6$.

Two representative examples with $\phi_{0}=10^{-3}$, $\phi_{0}=10^{-5}$ and
$n=2$ are shown in Fig. 5. Although the duration in $e$-folds grows by two
orders of magnitude when $\phi_{0}$ grows from $=10^{-5}$ to $10^{-3}$, the
background energy density follows a roughly similar exponential decay curve in
both cases. The lower panel of Fig. 5 indicates that the structure of the
instability bands can change considerably during preheating. The position of
the fastest-growing mode $k_{\text{peak}}$ changes by about an order of
magnitude when $\phi_{0}=10^{-3}$ and by almost two orders of magnitude when
$\phi_{0}=10^{-5}$.

In particular, when $\phi_{0}=10^{-3}$, the peak in the perturbation spectrum
from the grid computation matches well the estimate from the dominant peak
method, confirming its usefulness. Initially, $k_{\text{peak}}$ is located at
$k/a=1.1\times 10^{12}$ GeV, but as the background field amplitude decreases,
the peak and growth shift to larger wavenumbers. However, this shift is small:
at the end of the simulation, the maximum of the spectrum lies at
$k/a=3.0\times 10^{12}$ GeV, and the spectrum maintains a sharp, peaked form.
This justifies using the dominant peak approximation.

Figure 6: Duration of preheating in $e$-folds $N$, from the beginning of the
oscillations to the moment at which the energy density of perturbations
overtakes the background, as a function of the model variable $\phi_{0}$, for
different $n$. The corresponding initial oscillation amplitude
$\phi_{\text{amp},i}$ is also shown. _Left:_ Estimate from a numerical grid.
The benchmark points used in Fig. 5 are marked in red. _Right:_ Comparison of
the grid result to the dominant peak and simple approximations in their region
of validity $\phi_{0}\gtrsim 10^{-4}$. In the shaded regions, the results
become unreliable due to breakdown of adiabaticity for $\phi_{0}\lesssim
10^{-6}$ and end of tachyoncity for $\phi_{0}\gtrsim 10^{-2}$. Figure 7:
Different energy scales as functions of $\phi_{0}$. The initial oscillation
amplitude $\phi_{\text{amp},i}$ is also shown. _Left:_ The results using a
numerical grid showing the leading mode $k$ contributing most to the energy
density at the end of preheating, the initial fastest-growing mode
$k_{\text{peak},i}$, and the Hubble scale. The dashed lines indicate the
present frequency of potential GWs generated during preheating, assuming a
quick transition to radiation domination and a standard expansion history
after that. The benchmark points used in Fig. 5 are marked in red. _Right:_
Comparison of the final peak frequency on the grid and in the dominant peak
and simple approximations.

For $\phi_{0}=10^{-5}$, the peak in the power spectrum is much flatter, so the
dominant peak and the grid-based results differ significantly, which is why
Fig. 5 only presents the more accurate grid-based solution. The difference is
due to the larger range of $k$-values amplified during preheating. Initially,
the fastest-growing mode is $k/a=1.1\times 10^{11}$ GeV, while the final peak
in the spectrum sits at $k/a=4.5\times 10^{12}$ GeV, and the spectrum is
flatter than in the previous case. One reason for this is that, for larger
$\phi_{\text{amp},i}/\phi_{0}$, the secondary instability bands have a
relatively larger $\mu_{k}$, as can be seen in Fig. 3, and start to contribute
to the growth of the energy density of perturbations. In addition, the
oscillation amplitude changes more, so the tachyonic resonance band scans over
a wider range of modes. It follows that the modes at the final peak are
initially amplified by the weaker parametric resonance and then by the
stronger tachyonic instability at the end of preheating as the first peak
shifts towards smaller scales. The dominant peak approximation does not
capture the evolution of modes in the secondary instability bands and thus
severely underestimates preheating efficiency.

The simplified approximations fail when $\phi_{0}\lesssim 10^{-4}$. We find an
upper bound $\phi_{0}\lesssim 10^{-2}$ for preheating to be driven by the
tachyonic instability, which is consistent with earlier lattice computations
[33, 34]. When $\phi_{0}\gtrsim 10^{-2}$, the initial amplitude is small, so
the background exits the tachyonic region fast and reheating must complete
through other channels (right shaded region in Figs. 6 and 7). At the other
end of the scale is the lower bound $\phi_{0}\gtrsim 10^{-6}$ due to
adiabaticity – in this case, preheating completes within less than a half-
oscillation 171717We estimate the number of oscillations as
$n_{\mathrm{osc}}=\int{\rm d}t/T(\rho)$ in the presence of an evolving
background density, i.e., when $\dot{\rho}\neq 0$. (left shaded region in
Figs. 6 and 7). In fact, adiabaticity begins to be violated at somewhat larger
values $\phi_{0}\lesssim 10^{-5}$, in the sense that quantities like
$T\dot{\phi}_{\text{amp}}/\phi_{\text{amp}}$ are order one, so the amplitude
changes considerably during a half-oscillation. For
$10^{-5}\lesssim\phi_{0}\lesssim 10^{-2}$, this quantity stays below unity
almost until the end of the simulation, although the total number of half-
oscillations is never high; about $5$ for $\phi_{0}\sim 10^{-2}$. These order-
of-magnitude estimates apply for all our $n$-values.

Figs. 6 and 7 show the duration of preheating and the position of the peak of
the final spectrum for varying $\phi_{0}$, comparing the grid-based estimate,
the dominant peak approximation and the simple estimate. The value of $N$
varies between $10^{-5}$ and $10^{-1}$ in the eligible region and goes down
for smaller $\phi_{0}$. For the leading final $k$, this is also true for
$\phi_{0}\gtrsim 10^{-4}$, but for lower $\phi_{0}$, the peak value starts to
grow again due to the secondary, higher-$k$ peaks becoming more prominent. All
relevant $k$-values are well within the Hubble radius, more so for lower
$\phi_{0}$. These results agree with the earlier findings for $n=4$ [30] and
$n=2$ [32].

Both the dominant peak approximation and its simplified version predict values
of $N$ and the final leading $k$ that are of the correct order of magnitude
for $10^{-4}\lesssim\phi_{0}\lesssim 10^{-2}$, where the initial and final
leading $k$ are close to each other. As expected, the dominant peak
approximation performs better, being within a few percent of the grid result
in most of this region, but even the simple approach produces results within a
factor of two of the grid value for $k$ and $N$. The simple approach slightly
underestimates the leading $k$ and the preheating efficiency. Beyond their
region of validity, the approximations start to deviate from the grid,
severely underestimating $k$ and overestimating $N$.

###  3.4 Discussion

Consider the energy scales for preheating implied by the approximations given
above. Throughout the discussion we will use $m_{\rm th}=1.2\times
10^{13}\,{\rm GeV}$ fixed by the CMB measurement (2.12). Eqs. (3.7) and (2.22)
tell that the position of the dominant peak at the moment of preheating is

$k_{\text{peak}}/a_{\rm preh}\approx\sqrt{\phi_{0}}\,9\times 10^{13}\,{\rm
GeV}\,,$ (3.28)

in the region $10^{-4}\lesssim\phi_{0}\lesssim 10^{-2}$ in which the dominant
peak approximation works well. This expression has a mild dependence on $n$.
The above value was computed for $n=2$. When $10^{-6}\lesssim\phi_{0}\lesssim
10^{-4}$, the perturbation spectrum will peak at higher frequencies as can be
seen from Fig. 7. Assuming instantaneous thermalization, we can approximate
$\rho_{R}\approx U_{0}$, so the temperature after thermalization

$T_{\rm preh}\approx\sqrt{\phi_{0}}\,2\times 10^{15}\,{\rm GeV}$ (3.29)

exceeds the peak momentum by more than an order of magnitude, suggesting that
thermal effects will shift the initial spectrum towards higher modes.

Since tachyonic preheating is such a violent process, one may expect strong
production of secondary GWs. As Figs. 7 and 5 show, a potential GW signal is
created in the 1-10 GHz range when $\phi_{0}=10^{-4}-10^{-2}$. This range is
almost independent of $\phi_{0}$ or $n$. This frequency can be estimated by
redshifting the peak of the spectrum by a factor $a_{\rm now}/a_{\rm
preh}=(g_{*s,\rm preh}/g_{*s,\rm now})^{1/3}T_{\rm preh}/T_{\rm now}$, where
$T_{\rm now}=2.7$ K is the present photon temperature [59]. This gives

$f_{\rm GW}\approx k_{\text{peak}}/a_{\rm now}\approx 5\,{\rm GHz}$ (3.30)

and the dependence on $\phi_{0}$ has dropped out. These numbers are compatible
with the lattice studies [36, 37] performed for similar models. Such
frequencies are not observable with near-future gravitational-wave
interferometers [60, 61, 62], but may be probed through the effects on big
bang nucleosynthesis [63, 64, 65] or 21 cm measurements [65, 66].

Finally, let us take a closer look at the applicability of the linear theory
for small $\phi_{0}$. There are several theoretical arguments that combined
with our numerical results, point towards complications when $\phi_{0}\lesssim
10^{-6}$:

  * •

To estimate whether perturbation theory may be applied, consider the effective
quartic coupling $\lambda\sim U^{(4)}(\phi)\sim U_{0}/\phi_{0}^{4}=m_{\rm
th}^{2}/\phi_{0}^{2}$. Since $m_{\rm th}\approx 5\times 10^{-6}$ is fixed by
(2.12), the quadratic interaction $\lambda$ exceeds unity when
$\phi_{0}\lesssim 10^{-6}$. At this point, perturbativity is violated, and
non-linear contributions should not be neglected even at the beginning of
preheating when the energy density of fluctuations is negligible.

  * •

Due to the shape of the potential (2.5), the model is non-renormalizable even
when gravitational interactions are neglected. However, quantum corrections
can be treated in an effective field theoretic framework, so that all relevant
energy scales must stay below the cut-off scale of the effective field theory,
which for the potential (2.5) is approximately
$\Lambda_{\mathrm{cut}}=\phi_{0}$. From Fig. (7), we see that the peak of the
spectrum does not follow $k_{\text{peak},i}$ when $\phi_{0}\lesssim 10^{-4}$.
In fact, the peak of the perturbation spectrum scales roughly as $k/a\sim
m_{\rm th}$ rather than $k/a\sim m_{\rm th}\sqrt{\phi_{0}}$ as was estimated
in Eq. (3.7). Our numerical estimates show that the position $k$ of the peak
of the spectrum rise above the cut-off when $\phi_{0}\lesssim 10^{-6}$. The
same unitarity bound, $U_{0}\gtrsim\Lambda_{\mathrm{cut}}^{4}$, was derived in
Ref. [67], although not in the context of tachyonic preheating.

  * •

When the (non-renormalized) vacuum energy density of modes that can contribute
significantly to perturbations is larger than the background energy density,
then only a mild amplification of the vacuum fluctuations can result in the
total fragmentation of the field. An analogous argument was used in Ref. [36],
to derive the upper bound $4\times 10^{10}$ Hz for the frequency of GWs from
preheating.

Technically this means that $jT\gtrsim\rho$, so that the quantum source term
$j$ defined in Eq. (3.15) dominates. This is possible due to the contribution
of higher instability bands.181818Note that, due to Eq. (3.7), the (non-
renormalized) vacuum energy density around $k_{\text{peak},i}$ is always much
smaller than the initial background energy density, $\delta\rho^{\rm
vac}_{k_{\text{peak}}}\sim k_{\text{peak}}^{4}\sim U_{0}m_{\rm th}^{2}\ll
U_{0}$. As the condition implies significant energy transfer during a half-
oscillation, it is essentially equivalent to the breakdown of adiabaticity.
Indeed, neglecting backreaction and computing the evolution of linear
perturbations on the grid from the mode equation (3.1), $\delta\rho^{(R)}$
exceeds $\rho$ by the first zero-crossing of $\phi$ when $\phi_{0}\lesssim
10^{-6}$, implying, at the very least, an exit from the oscillatory regime. We
remark that the dominant peak estimate (3.23) for the number of oscillations
would severely underestimate the speed of fragmentation in this case.

Due to these complications, the parameter region $\phi_{0}\lesssim 10^{-6}$
seems to behave in a fundamentally non-linear way and thus, a more involved
analysis is indeed needed. The fact that only a mild amplification, or
equivalently, the production of a few quanta leads to fragmentation can also
be interpreted as the onset of a dominantly quantum regime. The implication of
the latter is that the $\phi_{0}\lesssim 10^{-6}$ region cannot be probed by
classical lattice studies.

##  4 Conclusions

We have studied preheating in inflationary scenarios in which the inflaton
possesses an exponentially flat plateau. We focused on the parameter region in
which the fragmentation of the inflaton is mainly due to a strong tachyonic
instability. Although perturbation growth was studied at the linear level, the
usual linear approach was extended by including backreaction effects between
the coherent oscillating background and the growing perturbations.

Tachyonic preheating in plateau inflation turns out to be a rapid process,
lasting only a fraction of an $e$-fold. This permitted us to treat cosmic
expansion adiabatically. Thus, at the leading order, the quantitative
characteristics of preheating can be determined in a Minkowski background. We
analytically derived characteristics of inflation, the oscillating background
and linear mode growth, which are universal to exponentially flat potentials.
In particular, the growth of individual modes lends itself to Floquet
theoretic treatment, and using this method, we numerically computed the
wavenumber and the growth rate of the fastest growing mode. In the tachyonic
regime, these quantities depend only weakly on the exact shape of the
exponentially flat potential.

During tachyonic preheating, the energy density of the oscillating background
gets damped mainly because its energy is transferred to the fragmented
component of the inflaton. The damping due to Hubble friction is always
subdominant. As the background’s oscillations proceed with ever decreasing
amplitude, the structure of the instability bands, which determines the growth
of each individual mode, changes. To account for the interplay between the
evolving background and mode growth, we constructed a set of coupled
continuity equations that appropriately account for the energy transfer
between these components in the adiabatic limit.

Mode growth in tachyonic preheating is characterised by a dominant instability
band at the lowest wavenumbers. We showed that further simplification is
possible when considering mode growth only from modes surrounding the fastest
growing mode. Within this approximation, it is possible to construct analytic
order of magnitude estimates that relate the global characteristics of
preheating, e.g. , its duration in $e$-folds, to the parameters of the model.
Using more accurate numerical estimates, we confirmed that this approximation
works well as long as adiabaticity is not violated.

In particular, we considered potentials of the form

$U(\phi)=U_{0}\tanh^{n}(\phi/\phi_{0})\,$

with even $n$, as representative cases of the general class of exponentially
flat potentials. For such potentials, CMB measurements fix $m_{\rm
th}\equiv\sqrt{U_{0}}/\phi_{0}=5\times 10^{-6}$, the half-period of background
oscillations is $T\approx m_{\rm th}^{-1}\phi_{0}^{-1/2}\ll H^{-1}$, and the
fastest growing mode as well as its growth rate are
$k_{\text{peak}}\approx\mu_{\text{peak}}\approx\pi/T$. Since
$\mu_{\text{peak}}\sim T^{-1}$, preheating is expected to complete in a few
oscillations. All of these expressions are nearly independently of $n$. The
qualitative features are completely controlled by the parameter $\phi_{0}$. In
particular, we find that

  * •

When $\phi_{0}\gtrsim 10^{-2}$, then mode growth is driven by parametric
resonance and the tachyonic instability is subdominant.

  * •

When $10^{-6}\lesssim\phi_{0}\lesssim 10^{-2}$, preheating is dominated by the
tachyonic instability, and the background evolves adiabatically even if
backreaction effects are taken into account. In this regime, preheating
completes in approximately $N\approx\sqrt{\phi_{0}}\lesssim 0.1$ $e$-folds. We
find that the analytic estimates work best in the parameter region
$10^{-4}\lesssim\phi_{0}\lesssim 10^{-2}$, while, when
$10^{-6}\lesssim\phi_{0}\lesssim 10^{-4}$, the adiabaticity begins to fail.
The latter leads to deviations between the simplified estimates considering
only the dominant peak in the mode growth spectrum and the more accurate
numerical estimates accounting for the exact growth spectrum.

  * •

When $\phi_{0}\lesssim 10^{-6}$, the background does not evolve adiabatically,
and the field can fragment significantly before completing a single half-
oscillation. Although the tachyonic instability would dominate preheating
after inflation, the damping of the background would rapidly activate the
parametric regime. In all, since the linear approach becomes meaningless,
$\phi_{0}\lesssim 10^{-6}$ is the region of non-linear preheating in which
quantum effects become relevant. As a related point, we estimate that
inflaton’s self-interactions are strong enough to violate unitarity in this
regime.

Our description inevitably fails near the end of preheating, when non-linear
effects become important, and for $\phi_{0}\lesssim 10^{-6}$ due to the loss
of adiabaticity. Rapid preheating is still expected, but the details have to
be resolved by a more comprehensive analysis. For $\phi_{0}\gtrsim 10^{-6}$,
this can be done with a more resource-intensive classical lattice computation,
like in [33, 34, 35, 36, 37]. In the future, our results can inform such
studies regarding the relevant time and energy scales. For $\phi_{0}\lesssim
10^{-6}$, classical lattice simulations are not applicable, since non-linear
effects come into play already in the quantum regime, after the production of
only a few quanta. New techniques are needed to resolve the dynamics in this
case.

Finally, preheating’s violent fragmentation process is expected to produce
high-frequency GWs in the 1-10 GHz range. This frequency range is roughly
determined by the fastest growing mode $k_{\text{peak}}\approx\pi/T$. Although
we do not address the GW spectrum and energy density in this paper, looking
for such signals can provide valuable information about the very first moments
in the early universe.

## Acknowledgments

This work was supported by the Estonian Research Council grants PRG803,
PRG1055, MOBTP135, MOBJD381 and MOBTT5 and by the EU through the European
Regional Development Fund CoE program TK133 “The Dark Side of the Universe.”

##  Appendix A Treatment of linear perturbations

We treat the linear inflaton perturbations quantum mechanically, given by the
field operator

$\delta\hat{\phi}(x)=\int\frac{{\rm
d}k^{3}}{(2\pi)^{3/2}}\quantity(\delta\phi_{k}\hat{a}_{\vec{k}}e^{i\vec{k}\cdot\vec{x}}+\delta\phi^{*}_{k}\hat{a}^{\dagger}_{\vec{k}}e^{-i\vec{k}\cdot\vec{x}})\,,\quad\commutator{\hat{a}_{\vec{k}}}{\hat{a}^{\dagger}_{\vec{p}}}=\delta^{3}(\vec{k}-\vec{p})\,,\quad\hat{a}_{\vec{k}}\ket{0}=0\,.$
(A.1)

The ladder operators $\hat{a}_{\vec{k}}$ define the vacuum sate $\ket{0}$
related to the mode functions $\delta\phi_{k}$. Time evolution is contained in
the mode functions $\delta\phi_{k}$, which follow the equations of motion
(3.1),

$\delta\ddot{\phi}_{k}+3H\delta\dot{\phi}_{k}+\omega_{k}^{2}\delta\phi_{k}=0\,,\qquad\omega_{k}^{2}\equiv\frac{k^{2}}{a^{2}}+U^{\prime\prime}(\phi)\,.$
(A.2)

and we use the Bunch–Davies initial conditions

$\delta\phi_{k}=\frac{1}{\sqrt{2k}a}\,,\qquad\delta\dot{\phi}_{k}=-i\frac{k}{a}\delta\phi_{k}\,,$
(A.3)

corresponding to the Minkowski-like adiabatic vacuum state. The expectation
value of the energy density of perturbations in the vacuum state is now

$\displaystyle\delta\rho$
$\displaystyle=\frac{1}{2}\expectationvalue{\delta\dot{\hat{\phi}}{}^{\dagger}\delta\dot{\hat{\phi}}+\frac{1}{a^{2}}\nabla_{i}\delta\hat{\phi}{}^{\dagger}\nabla_{i}\delta\hat{\phi}+U^{\prime\prime}(\phi)\delta\hat{\phi}{}^{\dagger}\delta\hat{\phi}}$
(A.4) $\displaystyle=\int\frac{{\rm
d}^{3}k}{(2\pi)^{3}}\quantity(\frac{1}{2}|\delta\dot{\phi}_{k}|^{2}+\frac{1}{2}\omega_{k}^{2}|\delta\phi_{k}|{}^{2})\,,$

which we use in (3.10).

Note that we have only included inflaton perturbations and neglected metric
scalar perturbations, even though these two are coupled at linear order. This
approximately corresponds to working in the spatially flat gauge, where the
effect of metric perturbations is minimized. In this gauge, the remaining
metric fluctuations play a subleading role for the leading sub-Hubble modes in
the tachyonic regime [30, 32], and we neglect them for simplicity.

In (A.2), the function $\omega_{k}^{2}$ depends on the background, which
(disregarding backreaction) follows the equations

$\ddot{\phi}+3H\dot{\phi}+U_{,\phi}=0\,,\qquad
3H^{2}=\frac{1}{2}\dot{\phi}^{2}+U\,,\qquad
U=U_{0}\tanh^{n}\frac{\phi}{\phi_{0}}\,$ (A.5)

introduced in section 2, where we specialized to our model with the potential
(2.5). To better understand the dependence of the dynamics on the input
parameters $\phi_{0}$ and $U_{0}$, it is useful to introduce the rescaled
variables

$\tilde{\phi}\equiv\frac{\phi}{\phi_{0}}\,,\quad\tilde{t}\equiv t\,m_{\rm
th}\,,\quad
E\equiv\frac{H}{\sqrt{U_{0}}}\,,\quad\tilde{U}\equiv\frac{U}{U_{0}}\,,\quad{\tilde{k}}\equiv\frac{k}{m_{\rm
th}}\,,$ (A.6)

where $m_{\rm th}\equiv\sqrt{U_{0}}/\phi_{0}$ from (2.6). In terms of these,
the equations read

$\displaystyle\tilde{\phi}^{\prime\prime}+3E\tilde{\phi}^{\prime}\phi_{0}+\tilde{U}_{,\tilde{\phi}}=0\,,\qquad
3E^{2}=\frac{1}{2}\tilde{\phi}^{\prime}{}^{2}+\tilde{U}\,,$ (A.7)
$\displaystyle\delta\phi^{\prime\prime}_{\tilde{k}}+3E\delta\phi^{\prime}_{\tilde{k}}\phi_{0}+\quantity(\frac{{\tilde{k}}^{2}}{a^{2}}+\tilde{U}_{,\tilde{\phi}\tilde{\phi}})\delta\phi_{\tilde{k}}=0\,,$
(A.8)

where a prime denotes a derivative w.r.t. the rescaled time $\tilde{t}$. The
potential height $U_{0}$ cancels out; it can be used to correctly fix the
amplitude $A_{s}$ of CMB, but does not affect the dynamics. The parameter
$\phi_{0}$ also only features in the friction terms, making them small, and
does not enter in the $H\to 0$ limit, e.g. , when computing the growth rate
$\mu_{k}$ (3.4) in section 3.1. Therefore, as $\phi$ can be measured in units
of $\phi_{0}$ and energy and time scales are measured in units of $m_{\rm th}$
(or a derived quantity like the oscillation time $T$), the explicit dependence
on $\phi_{0}$, $U_{0}$, and $m_{\rm th}$ can be removed by an appropriate
rescaling of dimensionful quantities. In particular, combinations like
$\mu_{k}T$ can be computed only by fixing a single background parameter, e.g.
, $\phi_{\text{amp}}/\phi_{0}$.

##  Appendix B Asymptotics for $H\to 0$ mode growth

In the following we will show that the growth of small $k$ modes for symmetric
potentials is determined by191919As expansion is irrelevant for the results of
this section, we set $a=1$.

$\frac{1}{2}\Tr
G=-1-\frac{k^{2}}{2}W\partial^{2}_{\rho}W+\mathcal{O}(k^{4})\quad\Rightarrow\quad\mu_{k}=k\real\frac{\sqrt{W\partial^{2}_{\rho}W}}{\partial_{\rho}W}+\mathcal{O}(k^{2})$
(B.1)

where $W$ is given by (2.14). For asymmetric potentials one must replace $\Tr
G\to-\Tr G$ and $W$ by the corresponding integral over the full period.

The large $k$ behaviour follows from treating the background adiabatically,
which gives

$\frac{1}{2}\Tr
G=\cos\left(\int^{T}_{0}\sqrt{k^{2}+U^{\prime\prime}}+\mathcal{O}(k^{-2})\right)\quad\Rightarrow\quad\lim_{k\to\infty}\mu_{k}=0$
(B.2)

with the damping of $\mu_{k}$ depending non-trivially on the shape of the
potential.

### The $k=0$ mode

The equation of motion (3.2) for the $k=0$ mode is

$\delta\ddot{\phi}_{0}-\frac{\dddot{\phi}}{\dot{\phi}}\delta\phi_{0}=0\,,$
(B.3)

where we used that $U^{\prime\prime}=-\dddot{\phi}/\dot{\phi}$ in a Minkowski
background. It follows that the system has two independent solutions

$u_{1}=\dot{\phi},\qquad u_{2}=\dot{\phi}\int\frac{{\rm
d}t}{\dot{\phi}^{2}}\,,$ (B.4)

where the second solution follows from the fact that the time derivative of
the Wronskian $\det w=u_{1}\dot{u}_{2}-u_{2}\dot{u}_{1}$ vanishes. This set of
solution satisfies $\det w=1$. Hill’s discriminant, i.e., the trace of the
monodromy matrix (3.5) is then

$\displaystyle\frac{1}{2}\Tr G^{(0)}$
$\displaystyle=\frac{1}{2}\left(u_{1}(0)\dot{u}_{2}(T)-\dot{u}_{1}(0)u_{2}(T)+u_{1}(T)\dot{u}_{2}(0)-\dot{u}_{1}(T)u_{2}(0)\right)$
(B.5)
$\displaystyle=\pm\frac{1}{2}\left(u_{1}(T)\dot{u}_{2}(T)-\dot{u}_{1}(T)u_{2}(T)+u_{1}(0)\dot{u}_{2}(0)-\dot{u}_{1}(0)u_{2}(0)\right)=\pm
1\,,$

where on the second line we used $\det w=1$ and $u_{1}(t+T)=\pm u_{1}(t+T)$,
depending whether $T$ is the period or half-period of $\phi$ (depending
whether the potential $U$ is asymmetric or symmetric, respectively). In both
cases, this implies that the ground state mode does not grow,

$\mu_{k=0}=0\,.$ (B.6)

We remark that (B.5) holds whenever the mode equation has a periodic solution
– the periodicity of the second solution $u_{2}$ is not necessary. This is
also clear from $\det G=1$, which implies that if one solution is oscillating,
then the other solution must have a vanishing growth rate. Despite that,
$u_{2}$ is not bounded and grows linearly. To see this we need to take a
closer look at the singularity in its defining integral (B.4) at the turning
point, i.e., when $\dot{\phi}=0$.

In the following, we will focus on symmetric potentials with $T$ the half-
period, so that $\dot{\phi}(t+T)=-\dot{\phi}(t)$. Let us pick $t=0$ at the
bottom of the potential, so that
$\ddot{\phi}(0)=-U^{\prime}(\phi_{\text{amp}})=0$. The turning point
$\phi_{\text{amp}}$ is reached at $t=T/2$, so $\dot{\phi}(T/2)=0$, and after
that, the background field rolls down the potential with
$\dot{\phi}(t)=-\dot{\phi}(T-t)$. In order to avoid integrating over the
singularity, we define

$u_{2}(t)=\left\\{\begin{array}[]{lc}\tilde{u}_{2}(t)&,\quad 0\leq t\leq
T/2\\\ \tilde{u}_{2}(T-t)+\beta\dot{\phi}(t)&,\quad T/2\leq t\leq
T\end{array}\right.$ (B.7)

with $\beta$ a constant and

$\tilde{u}_{2}(t)\equiv\dot{\phi}(t)\int^{t}_{0}\frac{{\rm
d}t^{\prime}}{\dot{\phi}(t^{\prime})^{2}}$ (B.8)

defined in the range $0\leq t\leq T/2$. With this choice of $u_{2}$, we have
that

$w(0)=\begin{pmatrix}\dot{\phi}(0)&0\\\
0&\dot{\phi}(0)^{-1}\end{pmatrix}\,,\qquad
w(T)=\begin{pmatrix}\dot{\phi}(T)&\dot{\phi}(T)\beta\\\
0&\dot{\phi}(T)^{-1}\end{pmatrix}$ (B.9)

and therefore, with $\dot{\phi}(T)=-\dot{\phi}(0)$, the monodromy matrix is
given by

$G_{k=0}\equiv w(0)^{-1}w(T)=\begin{pmatrix}-1&-\beta\\\ 0&-1\end{pmatrix}$
(B.10)

indicating a linear growth of a $k=0$ mode when $\beta\neq 0$, i.e., after $n$
half-oscillations $(G^{n})_{12}=(-1)^{n}n\beta$.

Our remaining task is to find the constant $\beta$. As long as $u_{2}$ is
continuous and differentiable, (B.7) will give a solution of (B.3) satisfying
$u_{1}\dot{u}_{2}-u_{2}\dot{u}_{1}=1$. The only instance that needs to be
checked separately is at $t=T/2$. First, at $t=T/2$,
$\tilde{u}_{2}(T/2)=-1/\ddot{\phi}(t)=1/U^{\prime}(\phi_{\text{amp}})$ and
$\dot{\phi}=0$ and thus $u_{2}$ is continuous as long as
$U^{\prime}(\phi_{\text{amp}})\neq 0$. The derivative $\dot{u}_{2}$ at $t=T/2$
is $\dot{\tilde{u}}_{2}(T/2)$ when $t\to T/2^{-}$ and
$-\dot{\tilde{u}}_{2}(T/2)+\beta\ddot{\phi}(T/2)$ when $t\to T/2^{+}$. Thus,
differentiability demands that

$\displaystyle\beta$
$\displaystyle=\frac{2\dot{\tilde{u}}_{2}(T/2)}{\ddot{\phi}(T/2)}=\lim_{t\to
T/2}\frac{2}{\ddot{\phi}(T/2)}\left[\ddot{\phi}(t)\int^{t}_{0}\frac{{\rm
d}t^{\prime}}{\dot{\phi}(t^{\prime})^{2}}+\frac{1}{\dot{\phi}(t)}\right]=-\partial_{\rho}T=-\partial_{\rho}^{2}W\,.$
(B.11)

The last two identities follow from (2.15) and

$\displaystyle\partial_{\rho}T$
$\displaystyle=2\partial_{\rho}\int^{\phi_{\text{amp}}}_{0}\frac{{\rm
d}\phi^{\prime}}{(2(\rho-U(\phi^{\prime}))^{1/2}}$ (B.12)
$\displaystyle=\lim_{\phi\to\phi_{\text{amp}}}2\left[-\int^{\phi}_{0}\frac{{\rm
d}\phi^{\prime}}{(2(\rho-U(\phi))^{3/2}}+\frac{\partial_{\rho}\phi_{\text{amp}}}{(2(\rho-U(\phi))^{1/2}}\right]$
$\displaystyle=\lim_{t\to T/2}2\left[-\int^{t}_{0}\frac{{\rm
d}t^{\prime}}{\dot{\phi}(t^{\prime})^{2}}-\frac{1}{\dot{\phi}(t)\ddot{\phi}(t)}\right]\,.$

### Small $k$ modes

To find the growth rates for small $k$, we look for perturbative solutions to
the mode equation (3.2) and expand

$\delta\phi_{k}=\sum_{i}(-k^{2})^{n}u^{(n)}$ (B.13)

with each subsequent order determined iteratively by

$\delta u^{(n)}-\frac{\dddot{\phi}}{\dot{\phi}}u^{(n)}=u^{(n-1)}\,.$ (B.14)

The leading order solution $u^{(0)}$ is given by any linear combination of
$u_{1,2}$. With boundary conditions $u^{(n)}(0)=\dot{u}^{(n)}(0)=0$ when
$n\geq 1$, the higher order solutions are found iteratively from

$u^{(n)}(t)=-u_{1}(t)\int^{t}_{0}{\rm
d}t^{\prime}\,u_{2}(t^{\prime})u^{(n-1)}(t^{\prime})+u_{2}(t)\int^{t}_{0}{\rm
d}t^{\prime}\,u_{1}(t^{\prime})u^{(n-1)}(t^{\prime}).$ (B.15)

Denoting by $u^{(n)}_{i}(t)$ the correction generated to $u^{(0)}=u_{i}$ and
using the fact that the correction as well as its derivative vanishes at
$t=0$, we can expand $\Tr G=\sum_{n}(-k^{2})\Tr G^{(n)}$.

The first correction is then

$\displaystyle\Tr G^{(1)}$
$\displaystyle=u_{1}(0)\dot{u}^{(1)}_{2}(T)-\dot{u}_{1}(0)u^{(1)}_{2}(T)+u^{(1)}_{1}(T)\dot{u}_{2}(0)-\dot{u}^{(1)}_{1}(T)u_{2}(0)$
(B.16)
$\displaystyle=(u_{2}(T)\dot{u}_{2}(0)-u_{2}(0)\dot{u}_{2}(T))\int^{T}_{0}{\rm
d}t^{\prime}\,u_{1}(t^{\prime})^{2}$

We used the periodicity of $\phi$ and the conservation of the Wronskian at the
$k^{2}$-order, $\det(w^{(0)}+w^{(1)})=1+\mathcal{O}(k^{4})$ which implies
$u_{1}\dot{u}^{(1)}_{2}-\dot{u}_{1}u^{(1)}_{2}=u_{2}\dot{u}^{(1)}_{1}-\dot{u}_{2}u^{(1)}_{1}$.
Plugging in $u_{1}=\dot{\phi}$, $u_{2}$ from Eq. (B.7) and using the
definition (2.14), we obtain

$\Tr G^{(1)}=W\partial_{\rho}T\,.$ (B.17)

#### Mode growth in the $k\to\infty$ limit

For large $k$, the changes in the potential can be treated adiabatically. In
particular, we require that $\dot{\omega}_{k}\ll\omega_{k}^{2}$, where
$\omega_{k}=\sqrt{k^{2}+U^{\prime\prime}}$. At leading order in $k$, this is
equivalent to $|U^{\prime\prime\prime}\dot{\phi}|\ll k^{3}$. In this case, the
two independent solutions of the mode equation (3.2) are approximately

$u_{1,2}=\omega_{k}^{-1/2}\exp(\pm i\int{\rm d}t\,\omega_{k})\,,$ (B.18)

so that

$\frac{1}{2}\Tr G=\cos\left(\int^{T}_{0}{\rm d}t\,\omega_{k}\right)\,,$ (B.19)

since $\omega_{k}$ is periodic in $T$ when the potential is symmetric (or $T$
when it is asymmetric). As $|\Tr G|\leq 2$, we find that $\mu_{k}=0$
confirming the well-known result that modes with an adiabatic evolution will
not grow.

##  Appendix C Some exact results for the $n=2$ potential

In $H\to 0$ limit, the field equation (2.2) admits an exact solution for the
potential (2.5) with $n=2$, i.e., when $U=U_{0}\tanh^{2}(\phi/\phi_{0})$. This
solution reads

$\phi(t)=\phi_{0}\,{\rm asinh}\left[\frac{\cos(\pi
t/T)}{\sqrt{U_{0}/\rho-1}}\right]\,,$ (C.1)

and

$W=\pi\sqrt{2U_{0}}\phi_{0}\left(1-\sqrt{1-\rho/U_{0}}\right),\qquad
T=\frac{\pi\phi_{0}}{\sqrt{2(U_{0}-\rho)}}$ (C.2)

and thus satisfies (2.24) with $C=1$.

The two independent solutions to the $k=0$ mode equation read

$u_{1}=-\frac{\sin\left(\frac{\pi
t}{T}\right)}{\sqrt{\frac{U_{0}}{\rho}-\sin^{2}\left(\frac{\pi
t}{T}\right)}},\qquad u_{2}=-\frac{\frac{U_{0}}{\rho}\cos\left(\frac{\pi
t}{T}\right)+\frac{\pi t}{T}\sin\left(\frac{\pi
t}{T}\right)}{\sqrt{\frac{U_{0}}{\rho}-\sin^{2}\left(\frac{\pi
t}{T}\right)}}\,.$ (C.3)

In this exact solution, the linear growth of $u_{2}$ is explicit. From (B.16)
it directly follows, that

$\frac{1}{2}\Tr G=-1-\frac{k^{2}T^{2}}{2}\left((1-\rho/U_{0})^{-1/2}-1\right)$
(C.4)

consistent with (3.9) when $C=1$.

## References

  * [1] A. A. Starobinsky, A New Type of Isotropic Cosmological Models Without Singularity, Phys. Lett. B 91 (1980) 99–102.
  * [2] A. A. Starobinsky, The Perturbation Spectrum Evolving from a Nonsingular Initially De-Sitter Cosmology and the Microwave Background Anisotropy, Sov. Astron. Lett. 9 (1983) 302.
  * [3] A. H. Guth, The Inflationary Universe: A Possible Solution to the Horizon and Flatness Problems, Phys. Rev. D 23 (1981) 347–356.
  * [4] A. D. Linde, A New Inflationary Universe Scenario: A Possible Solution of the Horizon, Flatness, Homogeneity, Isotropy and Primordial Monopole Problems, Phys. Lett. B 108 (1982) 389–393.
  * [5] A. Albrecht and P. J. Steinhardt, Cosmology for Grand Unified Theories with Radiatively Induced Symmetry Breaking, Phys. Rev. Lett. 48 (1982) 1220–1223.
  * [6] A. D. Linde, Chaotic Inflation, Phys. Lett. B 129 (1983) 177–181.
  * [7] D. H. Lyth and A. Riotto, Particle physics models of inflation and the cosmological density perturbation, Phys. Rept. 314 (1999) 1–146, [hep-ph/9807278].
  * [8] Planck Collaboration, Y. Akrami et al., Planck 2018 results. X. Constraints on inflation, Astron. Astrophys. 641 (2020) A10, [arXiv:1807.06211].
  * [9] J. Martin, C. Ringeval, R. Trotta, and V. Vennin, The Best Inflationary Models After Planck, JCAP 03 (2014) 039, [arXiv:1312.3529].
  * [10] BICEP2, Keck Array Collaboration, P. A. R. Ade et al., BICEP2 / Keck Array x: Constraints on Primordial Gravitational Waves using Planck, WMAP, and New BICEP2/Keck Observations through the 2015 Season, Phys. Rev. Lett. 121 (2018) 221301, [arXiv:1810.05216].
  * [11] T. Matsumura et al., Mission design of LiteBIRD, J. Low Temp. Phys. 176 (2014) 733, [arXiv:1311.2847].
  * [12] M. Hazumi et al., LiteBIRD: A Satellite for the Studies of B-Mode Polarization and Inflation from Cosmic Background Radiation Detection, J. Low Temp. Phys. 194 (2019), no. 5-6 443–452.
  * [13] CORE Collaboration, J. Delabrouille et al., Exploring cosmic origins with CORE: Survey requirements and mission design, JCAP 04 (2018) 014, [arXiv:1706.04516].
  * [14] L. Kofman, A. D. Linde, and A. A. Starobinsky, Reheating after inflation, Phys. Rev. Lett. 73 (1994) 3195–3198, [hep-th/9405187].
  * [15] L. Kofman, A. D. Linde, and A. A. Starobinsky, Towards the theory of reheating after inflation, Phys. Rev. D 56 (1997) 3258–3295, [hep-ph/9704452].
  * [16] G. N. Felder, J. Garcia-Bellido, P. B. Greene, L. Kofman, A. D. Linde, and I. Tkachev, Dynamics of symmetry breaking and tachyonic preheating, Phys. Rev. Lett. 87 (2001) 011601, [hep-ph/0012142].
  * [17] G. N. Felder, L. Kofman, and A. D. Linde, Tachyonic instability and dynamics of spontaneous symmetry breaking, Phys. Rev. D 64 (2001) 123517, [hep-th/0106179].
  * [18] J. Garcia-Bellido and E. Ruiz Morales, Particle production from symmetry breaking after inflation, Phys. Lett. B 536 (2002) 193–202, [hep-ph/0109230].
  * [19] E. J. Copeland, S. Pascoli, and A. Rajantie, Dynamics of tachyonic preheating after hybrid inflation, Phys. Rev. D 65 (2002) 103517, [hep-ph/0202031].
  * [20] N. Barnaby and J. M. Cline, Nongaussian and nonscale-invariant perturbations from tachyonic preheating in hybrid inflation, Phys. Rev. D 73 (2006) 106012, [astro-ph/0601481].
  * [21] N. Barnaby and J. M. Cline, Nongaussianity from Tachyonic Preheating in Hybrid Inflation, Phys. Rev. D 75 (2007) 086004, [astro-ph/0611750].
  * [22] J. Garcia-Bellido, D. G. Figueroa, and A. Sastre, A Gravitational Wave Background from Reheating after Hybrid Inflation, Phys. Rev. D 77 (2008) 043517, [arXiv:0707.0839].
  * [23] J.-F. Dufaux, G. Felder, L. Kofman, and O. Navros, Gravity Waves from Tachyonic Preheating after Hybrid Inflation, JCAP 03 (2009) 001, [arXiv:0812.2917].
  * [24] B. R. Greene, T. Prokopec, and T. G. Roos, Inflaton decay and heavy particle production with negative coupling, Phys. Rev. D 56 (1997) 6484–6507, [hep-ph/9705357].
  * [25] J. F. Dufaux, G. N. Felder, L. Kofman, M. Peloso, and D. Podolsky, Preheating with trilinear interactions: Tachyonic resonance, JCAP 07 (2006) 006, [hep-ph/0602144].
  * [26] A. A. Abolhasani, H. Firouzjahi, and M. M. Sheikh-Jabbari, Tachyonic Resonance Preheating in Expanding Universe, Phys. Rev. D 81 (2010) 043524, [arXiv:0912.1021].
  * [27] D. Battefeld, T. Battefeld, and J. T. Giblin, On the Suppression of Parametric Resonance and the Viability of Tachyonic Preheating after Multi-Field Inflation, Phys. Rev. D 79 (2009) 123510, [arXiv:0904.2778].
  * [28] N. Kitajima, J. Soda, and Y. Urakawa, Gravitational wave forest from string axiverse, JCAP 10 (2018) 008, [arXiv:1807.07037].
  * [29] H. Fukunaga, N. Kitajima, and Y. Urakawa, Efficient self-resonance instability from axions, JCAP 06 (2019) 055, [arXiv:1903.02119].
  * [30] J. Rubio and E. S. Tomberg, Preheating in Palatini Higgs inflation, JCAP 04 (2019) 021, [arXiv:1902.10148].
  * [31] A. Karam, M. Raidal, and E. Tomberg, Gravitational dark matter production in Palatini preheating, JCAP 03 (2021) 064, [arXiv:2007.03484].
  * [32] A. Karam, E. Tomberg, and H. Veermäe, Tachyonic preheating in Palatini R 2 inflation, JCAP 06 (2021) 023, [arXiv:2102.02712].
  * [33] K. D. Lozanov and M. A. Amin, Equation of State and Duration to Radiation Domination after Inflation, Phys. Rev. Lett. 119 (2017), no. 6 061301, [arXiv:1608.01213].
  * [34] K. D. Lozanov and M. A. Amin, Self-resonance after inflation: oscillons, transients and radiation domination, Phys. Rev. D 97 (2018), no. 2 023533, [arXiv:1710.06851].
  * [35] T. Krajewski, K. Turzyński, and M. Wieczorek, On preheating in $\alpha$-attractor models of inflation, Eur. Phys. J. C 79 (2019), no. 8 654, [arXiv:1801.01786].
  * [36] K. D. Lozanov and M. A. Amin, Gravitational perturbations from oscillons and transients after inflation, Phys. Rev. D 99 (2019), no. 12 123504, [arXiv:1902.06736].
  * [37] A. Bhoonah, J. Bramante, S. Nerval, and N. Song, Gravitational Waves From Dark Sectors, Oscillating Inflatons, and Mass Boosted Dark Matter, JCAP 04 (2021) 043, [arXiv:2008.12306].
  * [38] F. Bauer and D. A. Demir, Inflation with Non-Minimal Coupling: Metric versus Palatini Formulations, Phys. Lett. B 665 (2008) 222–226, [arXiv:0803.2664].
  * [39] J. Rubio, Higgs inflation, Front. Astron. Space Sci. 5 (2019) 50, [arXiv:1807.02376].
  * [40] V.-M. Enckell, K. Enqvist, S. Rasanen, and L.-P. Wahlman, Inflation with $R^{2}$ term in the Palatini formalism, JCAP 02 (2019) 022, [arXiv:1810.05536].
  * [41] I. Antoniadis, A. Karam, A. Lykkas, and K. Tamvakis, Palatini inflation in models with an $R^{2}$ term, JCAP 11 (2018) 028, [arXiv:1810.10418].
  * [42] T. Tenkanen, Tracing the high energy theory of gravity: an introduction to Palatini inflation, Gen. Rel. Grav. 52 (2020), no. 4 33, [arXiv:2001.10135].
  * [43] S. Ferrara, R. Kallosh, A. Linde, and M. Porrati, Minimal Supergravity Models of Inflation, Phys. Rev. D 88 (2013), no. 8 085038, [arXiv:1307.7696].
  * [44] R. Kallosh, A. Linde, and D. Roest, Superconformal Inflationary $\alpha$-Attractors, JHEP 11 (2013) 198, [arXiv:1311.0472].
  * [45] J. J. M. Carrasco, R. Kallosh, and A. Linde, Cosmological Attractors and Initial Conditions for Inflation, Phys. Rev. D 92 (2015), no. 6 063519, [arXiv:1506.00936].
  * [46] M. Galante, R. Kallosh, A. Linde, and D. Roest, Unity of Cosmological Inflation Attractors, Phys. Rev. Lett. 114 (2015), no. 14 141302, [arXiv:1412.3797].
  * [47] R. Kallosh and A. Linde, Universality Class in Conformal Inflation, JCAP 07 (2013) 002, [arXiv:1306.5220].
  * [48] J. Kim and J. McDonald, Inflaton Condensate Fragmentation: Analytical Conditions and Application to $\alpha$-Attractor Models, Phys. Rev. D 95 (2017), no. 12 123537, [arXiv:1702.08777].
  * [49] A. Lloyd-Stubbs and J. McDonald, Sub-Planckian $\phi^{2}$ inflation in the Palatini formulation of gravity with an $R^{2}$ term, Phys. Rev. D 101 (2020), no. 12 123515, [arXiv:2002.08324].
  * [50] J. Martin, C. Ringeval, and V. Vennin, Encyclopædia Inflationaris, Phys. Dark Univ. 5-6 (2014) 75–235, [arXiv:1303.3787].
  * [51] F. L. Bezrukov and M. Shaposhnikov, The Standard Model Higgs boson as the inflaton, Phys. Lett. B 659 (2008) 703–706, [arXiv:0710.3755].
  * [52] F. Bezrukov, The Higgs field as an inflaton, Class. Quant. Grav. 30 (2013) 214001, [arXiv:1307.0708].
  * [53] A. R. Liddle and S. M. Leach, How long before the end of inflation were observable perturbations produced?, Phys. Rev. D 68 (2003) 103503, [astro-ph/0305263].
  * [54] M. S. Turner, Coherent Scalar Field Oscillations in an Expanding Universe, Phys. Rev. D 28 (1983) 1243.
  * [55] J. H. Traschen and R. H. Brandenberger, Particle Production During Out-of-equilibrium Phase Transitions, Phys. Rev. D 42 (1990) 2491–2504.
  * [56] Y. Shtanov, J. H. Traschen, and R. H. Brandenberger, Universe reheating after inflation, Phys. Rev. D 51 (1995) 5438–5455, [hep-ph/9407247].
  * [57] N. D. Birrell and P. C. W. Davies, Quantum Fields in Curved Space.  Cambridge Monographs on Mathematical Physics. Cambridge Univ. Press, Cambridge, UK, 2, 1984.
  * [58] R. Reissig, Magnus, W. und S. Winkler, Hill’s Equation (Interscience Tracts in Pure and Applied Mathematics, Number 20). VIII + 127 S. New York/London/Sydney 1966. Interscience Publishers. Preis geb. 68 s.net, Zeitschrift Angewandte Mathematik und Mechanik 48 (Jan., 1968) 138–139.
  * [59] D. J. Fixsen, The Temperature of the Cosmic Microwave Background, Astrophys. J. 707 (2009) 916–920, [arXiv:0911.1955].
  * [60] KAGRA Collaboration, K. Somiya, Detector configuration of KAGRA: The Japanese cryogenic gravitational-wave detector, Class. Quant. Grav. 29 (2012) 124007, [arXiv:1111.7185].
  * [61] VIRGO Collaboration, F. Acernese et al., Advanced Virgo: a second-generation interferometric gravitational wave detector, Class. Quant. Grav. 32 (2015), no. 2 024001, [arXiv:1408.3978].
  * [62] B. P. Abbott et al., Sensitivity of the Advanced LIGO detectors at the beginning of gravitational wave astronomy, Phys. Rev. D 93 (2016), no. 11 112004, [arXiv:1604.00439]. [Addendum: Phys.Rev.D 97, 059901 (2018)].
  * [63] L. Pagano, L. Salvati, and A. Melchiorri, New constraints on primordial gravitational waves from Planck 2015, Phys. Lett. B 760 (2016) 823–825, [arXiv:1508.02393].
  * [64] A. Ejlli, D. Ejlli, A. M. Cruise, G. Pisano, and H. Grote, Upper limits on the amplitude of ultra-high-frequency gravitational waves from graviton to photon conversion, Eur. Phys. J. C 79 (2019), no. 12 1032, [arXiv:1908.00232].
  * [65] V. Domcke and C. Garcia-Cely, Potential of radio telescopes as high-frequency gravitational wave detectors, Phys. Rev. Lett. 126 (2021), no. 2 021104, [arXiv:2006.01161].
  * [66] A. Ringwald, J. Schütte-Engel, and C. Tamarit, Gravitational Waves as a Big Bang Thermometer, JCAP 03 (2021) 054, [arXiv:2011.04731].
  * [67] Y. Ema, R. Jinno, K. Nakayama, and J. van de Vis, Preheating from target space curvature and unitarity violation: Analysis in field space, Phys. Rev. D 103 (2021), no. 10 103536, [arXiv:2102.12501].