# (P)reheating and gravitational waves in $\alpha$-attractor models

Tomasz Krajewski  Krzysztof Turzyński

###### Abstract

We study post-inflationary evolution in $\alpha$-attractor T-models of
inflation. We consider the dynamics of both scalar fields present in these
models: the inflaton and the spectator, as a negative field-space curvature
may lead to geometrical destabilization of the spectator. We perform state-of-
the-art lattice simulations with a dedicated numerical code optimized for
those models. We corroborate earlier findings that the perturbations of the
spectator field are much more unstable than the perturbations of the inflaton
field, so the dynamics of early stages of preheating is dominated by the
evolution of spectator perturbations. We also calculate the spectrum of
gravitational waves originating from scalar fluctuations and find that the
$\alpha$-attractor T-models can be constrained or even ruled out by present
cosmological observations, but not by direct searches of gravitational waves.

##  1 Introduction

By now, cosmological inflation has become a firmly rooted part of the standard
cosmological model (see e.g. [1] for a pedagogical introduction). While the
increasing accuracy of the determination of cosmological parameters has
significantly narrowed down the range of empirically allowed inflationary
models, there are still a number of theoretically motivated possibilities. In
addition, there are also many ideas about the _reheating_ , i.e. the mechanism
of populating the Universe with degrees of freedom collectively satisfying the
equation of state characteristic of radiation. Although there exist possible
scenarios for reheating, well embedded in the general framework of quantum
field theory [2, 3], there are typically no observables that could at present
distinguish between various proposed scenarios.

One of the actively studied mechanisms of reheating is the so-called
_preheating_. It relies on the fact that the fluctuations of the inflaton and
other scalar fields – in the background of the homogeneous inflaton
oscillating around the minimum of its potential – can be unstable, which
potentially leads to mode amplification and particle production via a
parametric resonance [3, 4, 5].

A notable example of a class of models that are naturally equipped with the
preheating mechanism is the T-model subset of $\alpha$-attractor models of
inflation [6], which continue to attract a lot of theoretical interest. It is
known that in a single-field effective theory of the inflaton field stemming
from these models, the inflaton experiences self-resonance [7] at the end of
inflation and its perturbations may become highly unstable. Once they dominate
the Universe, the effective equation of state quickly approach that of
radiation [8, 9]. As a consequence, the radiation-dominated era effectively
begins very soon after inflation, which greatly reduces the theoretical
uncertainty stemming from interpreting cosmological observations in the
context of our ignorance about the reheating era [10, 11, 12, 13, 14, 15, 16,
17].

However, even in the minimal supergravity construction, one should in
principle consider both real degrees of freedom present in the scalar part of
the chiral multiplet. It turns out that there are interesting regions of the
parameter space of $\alpha$-attractor models which cannot be described by a
single-field effective theory. The second scalar field, residing in the chiral
multiplet of the supergravity construction typicaly plays the role of a
passive _spectator_. It is heavy when the modes that we now observe through
temperature anisotropies of the cosmic microwave background (CMB) radiation
leave the Hubble radius; its presence can, therefore, be safely neglected for
the calculation of the power spectrum of the curvature perturbations. This
assumption breaks down at the end of inflation, when the spectator field
becomes transiently tachyonic and unstable, so its perturbations may
eventually dominate the Universe. This is possible because the non-canonical
form of the kinetic part of the Lagrangian gives rise to geometrical
destabilization [18].

It is therefore interesting to investigate numerically the impact of the
spectator field on the dynamics of reheating, using the full two-field
approach. In the past, we made the first step in this direction [19], firmly
demonstrating that the spectator inhomogeneities quickly dominate over those
of the inflaton. As a consequence, we showed that preheating lasts even
shorter than previously thought. Moreover, due to nonlinear interactions, the
dramatically fast increase of the amplitude of spectator fluctuations is
translated to inflaton perturbations and the entire system is destabilized.
However, we have not been able to resolve the instability in sufficient detail
to track the effective equation of state of the Universe, because of the
limited lattice size.

In the present paper, we attempt to address this shortcoming of our previous
analysis by increasing the lattice size. We also study the spectrum of
gravitational waves that are produced during preheating out of scalar
fluctuations. We find that the instability is so strong that increasing the
linear lattice size by a factor of 4 (and thus increasing the lattice volume
by a factor of 64) compared to our previous study does not allow for resolving
all aspects of the dynamics of the unstable modes. We thereby confirm that the
instability of the spectator field is a crucial element of the evolution of
the Universe and that the preheating occurs almost instantaneously. We also
analyze the spectrum of the gravitational waves produced out of the scalar
perturbations and find that the amplitude of the gravitational waves grows
large enough that this class of models can be ruled out by the present
observations.

Our paper is organized as follows. In the section 2, we recapitulate the main
ingredients of the $\alpha$-attractor T-models that we consider. In the
section 3, we describe our methods of lattice simulations for scalar degrees
of freedom and present the results of those simulations. In the section 4, we
present our results for the spectrum of gravitational waves. Technical remarks
and supplementary results which help understand our methods and results, but
which are not necessary for inferring the main message of our paper are
relegated to appendices.

##  2 $\alpha$-attractor T-models of inflation

###  2.1 Presentation of the model

In this paper, we will consider the simplest and most thoroughly studied class
of $\alpha$-attractor models of inflation, the so-called T-models. They can be
naturally implemented in the supergravity with hyperbolic geometry as follows
[20, 21, 22, 23, 24, 6, 25, 26]. In so-called half-plane variables, $T$ and
$\bar{T}$, the Kähler potential takes the form

$K_{H}=-\frac{3\alpha}{2}\log{\left(\frac{(T-\bar{T})^{2}}{4T\bar{T}}\right)}+S\bar{S}$
(2.1)

with parameter $\alpha>0$ and the superpotential is

$W_{H}=\sqrt{\alpha}\mu S{\left(\frac{T-1}{T+1}\right)^{n}},$ (2.2)

where $n>0$ and $\mu$ is a constant parameter.

As shown in [25], the superfield $S$ can be stabilized during and after
inflation. Therefore, we will neglect its contribution into evolution of the
Universe and set $S=0$. Then, the scalar sector of the model can be expressed
in terms of two real scalar fields $\phi$ and $\chi$, which are related to the
scalar component of the superfield $T$ by

$\displaystyle\left|\frac{T-1}{T+1}\right|^{2}$
$\displaystyle=\left(\frac{\cosh(\beta\phi)\cosh(\beta\chi)-1}{\cosh(\beta\phi)\cosh(\beta\chi)+1}\right),$
where $\displaystyle\beta$ $\displaystyle:=\sqrt{\frac{2}{3\alpha}}.$ (2.3)

This choice is beneficial from the point of view of lattice simulations since
it leads to a particularly simple form of the field-space metric in the
kinetic term of the scalar Lagrangian:

$\mathcal{L}=-\frac{1}{2}\Big{(}\partial_{\mu}\chi\partial^{\mu}\chi+\textrm{e}^{2b(\chi)}\partial_{\mu}\phi\partial^{\mu}\phi\Big{)}-V(\phi,\chi).$
(2.4)

where $b(\chi)\equiv\log(\cosh(\beta\chi))$. The potential of the model in
this variables reads

$V(\phi,\chi)=M^{4}\left(\frac{\cosh(\beta\phi)\cosh(\beta\chi)-1}{\cosh(\beta\phi)\cosh(\beta\chi)+1}\right)^{n}{\left(\cosh(\beta\chi)\right)^{2/\beta^{2}}},$
(2.5)

with $M^{4}=\alpha\mu^{2}$.

With such a field-space metric, the reparametrization invariant field-space
curvature can be easily calculated

$\mathbb{R}=-2(b^{\prime 2}+b^{\prime\prime})=-2\beta^{2}=-\frac{4}{3\alpha},$
(2.6)

where ′ denotes the derivative of the function with respect to its argument.
We note that $\mathbb{R}$ is nonzero, thus one cannot canonically normalize
both degrees of freedom at once. Moreover, $\mathbb{R}$ is negative, which is
crucial for the possibility of geometrical destabilization. Finally, the value
of $\mathbb{R}$ is determined by $\alpha$, and the inverse proportionality
means that the size of $\mathbb{R}$ increases with decreasing value of
$\alpha$ parameter, thus one expect the stronger destabilization for smaller
values of $\alpha$.

###  2.2 Inflationary trajectory

It was shown in [6] that T-models presented in Section 2.1 admit an inflating
solution with the inflationary trajectory proceeding along $\chi=0$. Moreover,
in [27] the case $\chi\neq 0$ was investigated and it was shown that there are
two stages of inflation; the first one proceeds nearly along $\phi=const$ and
the second one with $\chi$ being exponentially close to $0$. However, for
reasonable values of the ratio $\frac{\chi_{0}}{\phi_{0}}$ of initial values
of the field strengths, the number of e-folds of the second stage is much
larger than that of the first stage. Therefore, the generic observational
prediction of these models for the CMB anistotropies is that of the usual
single-field inflation. Remarkably, the model is consistent with Planck data
for a wide range of its parameters (see e.g. [23]).

The inflationary trajectory can be effectively described in terms of the
canonically normalized inflaton field $\phi$ with the potential

$V(\phi,0)=M^{4}\tanh^{2n}\left(\frac{\beta|\phi|}{2}\right).$ (2.7)

Setting $\chi=0$ during inflation, the trajectory is the solution of following
equations of motion for the space homogeneous quantities $H:=\dot{a}/a$ and
$\phi$:

$\displaystyle H^{2}$
$\displaystyle=\frac{1}{3}\left(\frac{1}{2}\dot{\phi}^{2}+V(\phi,0)\right),$
$\displaystyle\ddot{\phi}+3H\dot{\phi}+V_{\phi}(\phi,0)$ $\displaystyle=0,$
(2.8)

where $V_{\phi}$ stands for a derivative in the $\phi$ direction
($V_{\phi}\equiv\frac{\partial V}{\partial\phi}$) and $\dot{}$ denotes the
derivative with respect to cosmic time $t$. We introduced $a(t)$ for the scale
factor of Friedmann-Robertson-Walker metric.

As usual in single-field models, the inflaton accelerates during inflation and
eventually leaves the slow-roll regime defined as
$\epsilon_{H}\equiv-\dot{H}/H^{2}\ll 1$ and
$\eta_{H}\equiv-\ddot{H}/(2\dot{H}H)\ll 1$. In our considerations, we assume
that the end of the inflation is indicated by the condition $\epsilon_{H}=1$.

At same time when the velocity of the inflaton grows large enough to go beyond
the slow-roll regime, a negative value of field space curvature can cause so-
called geometrical destabilization [18] of perturbations of the field $\chi$.
Therefore, in order to track the dynamics of perturbations accurately, both
fields should be taken into account around the end of the inflation. In the
section 2.3 we will briefly present linear analysis of perturbations of
inflaton $\phi$ and spectator field $\chi$. Our results of nonlinear numerical
lattice simulations will be described in the section 3.3.

###  2.3 First-order perturbations

Equations of motion for perturbations in two-field models described by (2.4)
can be found e.g. in [28], without slow-roll approximation or any additional
assumptions. The perturbed Friedmann-Robertson-Walker metric (in longitudinal
gauge, with only scalar degrees of freedom included and constraints taken into
account) reads

$\mathrm{d}s^{2}=-(1+2\Psi)\,\mathrm{d}t^{2}+a^{2}(1-2\Psi)\,{\mathrm{d}\mathbf{x}^{2}}.$
(2.9)

Linear perturbations are described in terms of gauge-invariant Mukhanov-Sasaki
variables:

$\displaystyle Q_{\phi}$ $\displaystyle:=\delta\phi+\frac{\dot{\phi}}{H}\Psi$
and $\displaystyle Q_{\chi}$
$\displaystyle:=\delta\chi+\frac{\dot{\chi}}{H}\Psi,$ (2.10)

which, for $\chi=0$, obey very similar equations of motion

$\ddot{Q}_{\phi/\chi}+3H\dot{Q}_{\phi/\chi}+\left(\frac{k^{2}}{a^{2}}+m^{2}_{\phi/\chi}\right)Q_{\phi/\chi}=0,$
(2.11)

with the only difference in effective masses

$\displaystyle m^{2}_{\phi}$ $\displaystyle=V_{\phi\phi},$ $\displaystyle
m_{\chi}^{2}$
$\displaystyle=V_{\chi\chi}+\frac{1}{2}\dot{\phi}^{2}\mathbb{R}.$ (2.12)

Writing eqs. (2.11) and (2.12), we neglected all contributions suppressed by
the Planck scale, as the energy scale of inflation is much smaller. We also
used the assumption $\chi=0$ which immediately implies that
$V_{\chi}(\phi,0)=0$ and $V_{\phi\chi}(\phi,0)=0$ for $V$ given by (2.5).

Taking into account the fact that the field space curvature is negative,
according to eq. (2.6), one deduces from eq. (2.12) that for large values of
$\beta$, i.e. for small values of $\alpha$, perturbations $Q_{\chi}$ are
expected to exhibit an _intermittent tachyonic instability_ as $|\dot{\phi}|$
increases towards the end of inflation. The linear analysis of evolution of
perturbations of both the inflaton $Q_{\phi}$ and the spectator $Q_{\chi}$
during preheating, performed on the basis of the Floquet theorem in [19] and
extended in [27], confirmed this prediction. Moreover, the so-called Floquet
exponents, which measure speed of the exponential growth of perturbations
during instability in oscillating background of the inflaton $\phi$, turned
out to be much larger for the spectator than for the inflaton, revealing that
the described _intermittent tachyonic instability_ is an important factor for
mode amplification and particle production, thereby competing with the
parametric resonance of the inflaton perturbations, which was found to be
effective for reheating even for single-field model [8, 9].

In [27], the authors described an intriguing simple scaling of the instability
of linear spectator perturbations $Q_{\chi}$ (especially the Floquet
exponents) with changes of $\alpha$ in regime of high curvature of the field
space. Their numerical calculation of the Floquet exponents as functions of
$\phi_{0}/\sqrt{\alpha}$ and $k/\mu=k\sqrt{\alpha}M_{Pl}/M^{2}$ give very
similar results for $\alpha=10^{-3},10^{-4}$; those results are also quite
similar for those for $\alpha=10^{-2}$.

Analytic arguments explaining the observed universality can be briefly
summarized, as follows. The potential (2.5) is a fixed function of
$\delta=\phi/\sqrt{\alpha}$ for $\chi=0$. The calculation in slow-roll
approximation shows that the fraction
$\delta_{\text{end}}=\phi_{\text{end}}/\sqrt{\alpha}$ of the value of the
inflaton field at the end of inflation $\phi_{\text{end}}$ (and the beginning
of oscillations) has only logarithmic dependence on $\alpha$ (and on $n$). In
fact, this dependence leads to residual one of other quantities. Furthermore,
it turns out that the value of parameter $\mu=M^{2}/\sqrt{\alpha}$ is
determined independently of $\alpha$ from normalization of scalar power
spectrum. Hence, the dynamics of background oscillations of $\delta$ is
independent of $\alpha$ up to dependence of Hubble parameter $H$ which can be
neglected in the first approximation since frequency of oscillations $\omega$
for small $\alpha$ is much higher then Hubble time. As a result, the
oscillations of $\delta$ have the amplitude which is logarithmically dependent
on $\alpha$ and they are weakly damped (proportionally to $H$). This leads to
conclusion that the effective mass $m_{\chi}^{2}$ of the spectator
perturbations is practically independent of $\alpha$ when expressed in terms
of $\delta$ and $\mu$, since the dependence on $\alpha$ cancels out between
$\mathbb{R}$ and $\dot{\phi}^{2}$ (similarly, dependence from derivatives with
respect to $\chi$ cancels the one from $M^{4}$ in $V_{\chi\chi}$ term).
Neglecting the time dependence of the scale factor $a$ one obtains the Floquet
exponents expressed as functions of $\delta_{0}$ and $\mu$, without an
explicit dependence on $\alpha$.

The dependence of the Floquet exponents on potential parameter $n$ is also
discussed in [27], but it is not as straightforward as the dependence on
$\alpha$. The similarities are visible only for highly curved field spaces
with $\alpha\lesssim 10^{-4}$ when the reheating is essentially instantaneous
(and with $\phi/\phi_{\text{end}}$ used instead of $\delta$).

The authors of [27] also estimated the duration of reheating, calculating the
time required to transfer the entire inflaton energy density into the energy
of spectator perturbations, extrapolating linear analysis outside its validity
range. Therefore, that estimation has obvious limitations. The perturbative
expansion will break down as soon as the produced fluctuations will backreact
on the background evolution. Moreover, the modes with large wavenumbers will
be produced in re-scatterings due to nonlinearity of the potential (2.5)
leading to broadening of the energy spectrum of fluctuations. These effects
can be studied to some extent in numerical lattice simulations, as we discuss
in the section 3.3.

##  3 Numerical lattice simulations

###  3.1 Description of methods

We developed own code directly designed for numerical studies of the evolution
of the scalar perturbations in $\alpha$-attractor models on the lattice. Our
main rationale was that the existing codes written to simulate preheating
after inflation on the lattice (see e.g [29] for a review) are not suitable
for non-canonical kinetic terms. The only exception is GABE, which is,
however, non-symplecic and thus poorly reproduce energy-conservation,

The theoretical background of our approach is described in detail in the
appendix A of [19]. We used second order central finite difference scheme to
discretize in space the continuum theory given by the Lagrangian density
(2.4). Introducing canonical momenta $\pi_{\phi}$, $\pi_{\chi}$, $p_{a}$
conjugate to $\phi$, $\chi$ and $a$, respectively, we formed a discrete
Hamiltonian describing the dynamics of finite number of degrees of freedom in
the discretized theory. It turns out that such a discrete Hamiltonian
(analogously as its continuum counterpart) is separable and can be divided
into four terms in such a way that none of the parts depends on both the
canonical variable and its conjugate momentum. Hamilton’s equation for each
term separately can be integrated exactly. Using collocation method for these
four exact flows we formed symplectic integrator of a second-order precision
in time.

In order to control the error made during numerical computation and guarantee
the reliability of results we used an error indicator based on a priori error
estimation to detect the breakdown of used discretization due to production of
fluctuations with wavelengths shorted then lattice spacing $h$ (i.e. physical
distant between neighbouring lattice points). Moreover, we performed 5
independent simulations for each set of parameters which we compared and
verified their consistency. This attempt guarantee full reliability of the
presented results in contrast to consistency verification alone which cannot
detect some systematic error.111In addition we unit tested our numerical code,
thus the existence of unknown bugs in it is unlikely.

The main development in our code with respect to [19] is the usage of GPU
acceleration to speed up our computations, allowing us to increase the size of
lattice from $128^{3}$ used previously to $512^{3}$. Taking advantage of the
obtained decrease in computation time, we were able to study deeper the
parametric dependence of geometrical destabilisation. We simulated the process
of the preheating for 12 different parameters sets, with $\alpha=10^{-3}$,
$10^{-3.5}$ and $10^{-4}$, and $n=1,1.5,2$ and $3$. For each set of
$(\alpha,n),$ we performed simulations with three different lattice spacings,
which we parametrize by the maximal wavevector $k_{\text{max}}=h/\sqrt{3}\pi$
that fits on a cubic lattice of spacing $h$. Comparing results of simulations
with values of $k_{\text{max}}$ that differ by the factor of $2$ we were able
to study the effects of spatial discretization of equations of motion and the
role of IR and UV cutoffs that are inevitably associated with lattice
simulations. For each parameter set, we chose a value of $k_{\text{max}}$ that
ensured a good trade-off between granularity in the Floquet instability
regions and factoring in higher frequency modes. Then, simulations with lower
cutoff confirm that the production of fluctuations due to first instability
band is properly reproduced and ones with higher cutoff that modes with
wavelengths lower than UV cutoff do not have significant influence on the
evolution during simulated period of time.

The evolution of fields $\phi$ and $\chi$ obtained from our lattice
simulations was also used to determine the spectrum of gravitational waves
emitted during (p)reheating. The details of this calculation are described in
the section 4.

###  3.2 Initial conditions for perturbations

At this point we take an advantage of the eq. (2.11) and describe initial
conditions for our numerical simulations.

One usually chooses the Bunch-Davies initial conditions for perturbations
which are appropriate for quantum fields in time-dependent, de Sitter
background. This procedure is well-known for fields with trivial field-space
metric and was generalized for non-trivial cases, e.g. for two-field models
with kinetic terms as in the eq. (2.4), the prescription for initial
conditions of the perturbations can be found e.g. in [28]. It simplifies when
$\chi=0$ and when one can neglect the contributions suppressed by the Planck
scale.

Defining $u_{\phi}:=aQ_{\phi}$ and $u_{\chi}:=aQ_{\chi}$, and using conformal
time $\tau$ (related to the cosmic time $t$ as
$\mathrm{d}\tau=a^{-1}\mathrm{d}t$), eqs. (2.11) can be written as equations
of motion of two independent harmonic oscillators with time dependent masses
(for details, see e.g. [30]). Therefore, a usual single-field quantization
procedure (see, e.g. [31]) can be performed. In the so-called adiabatic
approximation, it provides the following initial conditions for perturbations
at $\tau=\tau_{0}$:

$\displaystyle u_{I}(k,\tau_{0})$
$\displaystyle=\frac{1}{\sqrt{2\omega_{I,k}}}\textrm{e}^{-i\omega_{I,k}\tau_{0}},$
$\displaystyle u^{\prime}_{I}(k,\tau_{0})=$
$\displaystyle-i\sqrt{\frac{\omega_{I,k}}{2}}\textrm{e}^{-i\omega_{I,k}\tau_{0}},$
(3.1)

where $I\in\\{\phi,\chi\\}$, the prime ′ denotes a derivative with respect to
the conformal time and the frequencies $\omega_{I,k}$ are defines as:

$\omega_{I,k}^{2}=k^{2}+a^{2}\left(m^{2}_{I}+a^{-3}\partial_{\tau}^{2}a\right)\approx
k^{2}+a^{2}\left(m^{2}_{I}+2H^{2}\right),$ (3.2)

where in the last approximate equality, the terms proportional to slow-roll
parameters were neglected.

We use expressions (3.1) to set Gaussian initial conditions for perturbations
in our lattice simulations in the similar way as it is done for LATTICEEASY
code [32] (however later we transform velocities into conjugate momenta).

The adiabatic approximation can be used when the following condition is
satisfied:

$\frac{\omega_{I}^{-1}}{2}\frac{\mathrm{d}^{2}\omega_{I}^{-1}}{\mathrm{d}\tau^{2}}-\frac{1}{4}\left(\frac{\mathrm{d}\omega_{I}^{-1}}{\mathrm{d}\tau}\right)^{2}\ll
1\,.$ (3.3)

The initial time $\tau_{0}$ for the lattice simulations must be chosen so that
the condition (3.3) holds true.

###  3.3 Results of simulations

We start with the description of results of simulations with $n=1$ and
$n=1.5$, and for $\alpha$ ranging from $10^{-4}$ to $10^{-3}$, since for those
sets of parameters, lattice simulations were previously performed by us and
discussed in [19].

From the point of view of comparing CMB measurements with predictions of
inflationary models, the barotropic parameter $w=p/\rho$ is a very useful
quantity, as it is directly connected to the way the Universe expands. In our
model, the barotorpic parameter can be written in terms of fields as:

$w=\frac{\frac{1}{2}\left(\textrm{e}^{2b(\chi)}\dot{\phi}^{2}+\dot{\chi}^{2}\right)-\frac{1}{6a^{2}}\left(\textrm{e}^{2b(\chi)}(\nabla\phi)^{2}+(\nabla\chi)^{2}\right)-V(\phi,\chi)}{\frac{1}{2}\left(\textrm{e}^{2b(\chi)}\dot{\phi}^{2}+\dot{\chi}^{2}\right)+\frac{1}{2a^{2}}\left(\textrm{e}^{2b(\chi)}(\nabla\phi)^{2}+(\nabla\chi)^{2}\right)+V(\phi,\chi)}\,,$
(3.4)

where spatial averaging is assumed for both the numerator and the denominator.

The evolution of $w$ is presented in the figure 1 for $n=1$ and $n=1.5$. It
can be seen that the value of barotropic parameter oscillates, which begins
within the end of the inflation, corresponding to $N=0$, due to the periodic
transfer of the energy between the kinetic and the potential energy of the
inflaton. The initial amplitude of the oscillations is from $-1$ to $1$.
However, due to instability the amplitude decreases, because the energy of the
homogeneous background is transferred into fluctuations with gradient
component of the energy. For all investigated sets of parameters, this
decrease is visible and, with the exception of $n=1$ and $\alpha=10^{-3}$, it
starts rapidly after a few oscillations of the background inflaton field, in
just a fraction of the e-fold after the end of the inflation. As we argued in
[19], and confirm in the present paper, rapid evolution of the barotropic
parameter is caused by efficient production of fluctuations of fields due to
geometrical destabilization.

Figure 1: Evolution of the barotropic parameter $w$ (blue solid line) as a
function of the number of e-folds $N$ from the end of inflation ($\epsilon=1$)
for $n=1$ (left column) and $n=1.5$ (right column) and $\alpha=10^{-3}$ (top),
$\alpha=10^{-3.5}$ (middle), $\alpha=10^{-4}$ (bottom). Dashed orange vertical
line correspond to $w=1/3$.

Comparing the figure 1 with the figure 6 in [19], it can be seen that the
improved spatial resolution allows for longer and more accurate simulations,
but our lattice size is not yet sufficient to proclaim the fate of the
barotropic parameter after inflation, as we are unable to determine the
asymptotic value of $w$.

More insight about this matter can be gained from studying the evolution of
various components of the energy density, which we present in the figure 2.
The initial evolution of energy density is the same as in single-field models.
From the beginning of the simulation, the potential energy $E_{pot}$ decreases
slowly, translating into the kinetic energy $E_{kin,\phi}$ of the inflaton.
From the end of inflation, the potential energy density starts to oscillate
and the energy is exchanged between the potential and the kinetic components.
However, around that time, geometrical destabilization becomes important and
both the kinetic $E_{kin,\chi}$ and gradient $E_{grad,\chi}$ energy densities
of the spectator start increasing. Such a simultaneous growth at similar rate
indicates that not only the homogeneous mode, but mainly the spatial
fluctuations of the spectator are produced. Finally, after a fraction of an
e-fold after the end of inflation, the gradient energy density $E_{grad,\phi}$
of the inflaton starts to grow rapidly. This is a symptom of the production of
fluctuations of the inflaton caused by the nonlinear interactions with the
induced modes of the spectator field, the effect that cannot be deduced from
linear Floquet analysis. During further evolution, the energy density stored
in fluctuations of both fields increases. With the exception of $n=1$ with
$\alpha=10^{-3}$, we can follow the evolution of the fields until the kinetic
and gradient energy density of the spectator is within one order of magnitude
from the average energy density of the homogeneous inflaton field.

Unfortunately, the investigated instability is so strong that the used
discretization of equations of motion on the constant lattice breaks down due
to production of short wavelength modes (shorter than the lattice spacing $h$)
and we have to stop our simulations around the time when the energy density of
fluctuations is dominant component since the further calculations would give
unreliable results. Yet our calculations clearly demonstrate that very shortly
after the simulated phase of the evolution of the Universe the single field
description of the considered model breaks down due to the dynamics of the
spectator.

Figure 2: Evolution of energy components as a function of the number of
e-folds $N$ from the end of inflation for $n=1$ (left column) and $n=1.5$
(right column) and $\alpha=10^{-3}$ (top), $\alpha=10^{-3.5}$ (middle),
$\alpha=10^{-4}$ (bottom).

The evolution of fluctuations of fields can be also studied with Fourier
transform of fields strength. In the figure 3, we present the evolution of
amplitudes of the Fourier transform of the spectator field $\chi$ for selected
wave vectors $k$. Simultaneously with the end of inflation the amplitude of
the lowest $k$ mode available on the lattice (deep red) starts to grow by
orders of magnitude. On the other hand, high $k$ modes stay nearly unaffected
for some time until they rapidly increase close to the end of simulation.
Modes for some intermediate values of $k$ develop in between.

Figure 3: Evolution of amplitudes of selected modes of fluctuations
$\delta\chi$ of the spectator field as a function of the number of e-folds $N$
from the end of inflation ($\epsilon=1$) for $n=1$ (left column) and $n=1.5$
(right column) and $\alpha=10^{-3}$ (top), $\alpha=10^{-3.5}$ (middle),
$\alpha=10^{-4}$ (bottom).

Some further characteristics of the instability can be recognized from spectra
of fluctuations at various fixed times. We present those spectra in the figure
4. We can see that at the beginning of the destabilization only modes with
wave vectors $k$ smaller than the scale of the negative mass squared of the
spectator are excited. Such pattern is predicted by the linear Floquet
analysis, because these modes are excited by the tachyonic instability caused
by the negative contribution to the effective mass coming from negative
curvature of the field space. Moreover, a ‘ringing’ pattern resulting from
rescattering of amplified modes is clearly visible. Near the end of our
simulations, these low-$k$ modes produced due to tachyonic instability
backreact and cause the growth of modes with larger wave vector $k$ including
finally all modes available on the lattice. In contrast to the initial stage
of the destabilization, which can be modelled via linear approximation,
amplification of high $k$ modes is a purely nonlinear effect. However, this
effect plays a very important role in reheating, since it leads to the
fragmentation of fields and, consequently, to the transition to the radiation
domination era. It is worth stressing that the maxima of the spectra of
spectator fluctuations for all simulations are well separated from the $k=0$,
i.e. the value at the maximum is always orders of magnitude higher than the
value for the lowest $k$ available on the lattice. This proves that the
dynamics of low-$k$ modes is reliably reproduced in all our simulations during
the considered period of time, thus the IR cutoff is optimally chosen.

Figure 4: Evolution of the spectrum of fluctuations $\delta\chi$ of the
spectator field as a function of the number of e-folds $N$ from the end of
inflation for $n=1$ (left column) and $n=1.5$ (right column) and
$\alpha=10^{-3}$ (top), $\alpha=10^{-3.5}$ (middle), $\alpha=10^{-4}$
(bottom). Lines of different colors correspond to different time points
described by the number of e-folds after the end of inflation, as indicated on
the bar next to each panel.

For completeness, we also show analogous plots for fluctuations $\delta\phi$
of the inflaton in the figure 5. Again, we can see for most of the simulated
period of time, only low-$k$ modes of the inflaton fluctuations are excited.
This observation is consistent with the Floquet analysis, since the inflaton
is affected only by the parametric resonance [8, 9, 19], which is much weaker
than the tachyonic instability of the spectator. However, towards the end of
the simulation higher wave vector modes of $\delta\phi$ follow those of
$\delta\chi$, since the two fields $\phi$ and $\chi$ are tightly coupled
through the non-canonical kinetic term. Although the amplitudes of modes of
the fluctuations $\delta\phi$ of the inflaton are smaller, the contributions
to the energy density coming from them dominate over those of $\delta\chi$ due
to the factor of $\cosh^{2}(\beta\chi)$ multiplying the kinetic term of
$\phi$.

Figure 5: Evolution of the spectrum of fluctuations $\delta\phi$ of the
spectator field as a function of the number of e-folds $N$ from the end of
inflation for $n=1$ (left column) and $n=1.5$ (right column) and
$\alpha=10^{-3}$ (top), $\alpha=10^{-3.5}$ (middle), $\alpha=10^{-4}$
(bottom). Lines of different colors correspond to different time points
described by the number of e-folds after the end of inflation, as indicated on
the bar next to each panel.

So far, we have discussed the dynamics of geometrical destabilization at the
end of inflation, taking as examples models with $n=1$ and $n=1.5$. Now we
shall also address models with $n=2$ and $n=3$, in part to study the parameter
dependence of that phenomenon. In the figure 6, we show the evolution of the
barotropic parameter for models with $n=2$ and $n=3$. Comparing these plots
with those in the figure 1, we see that decreasing the value of $\alpha$
parameter decreases the number of e-folds needed for a prospective transition
to radiation-domination era. With the notable exception of the case $n=1$,
$\alpha=10^{-3}$, this number of e-folds is smaller than $1$ for $\alpha\leq
10^{-3}$; thus we are corroborating our claim from [19] that for small values
of $\alpha\leq 10^{-3}$, reheating in two-field T-models is practically
instantaneous. The duration of reheating decreases with increasing value of
$n$. In all cases, a shorter transition time is associated with a decreasing
number of oscillations of background inflaton field needed to transfer energy
to fluctuations.222The number of oscillations of the inflaton can be directly
counted from zeros of potential energy density in figures 2 and 10. However,
each oscillation before the transition is associated with the oscillation of
barotropic parameter $w$. These observations are in agreement with predictions
of approximate analysis of ref. [27]. This apparent concordance hints at a
certain subtlety: in our simulations we can see the onset of nolinear
backreaction, which, in particular, leads to excitations of the inflaton modes
through rescattering and is vital for the evolution of the barotropic
parameter. Our numerical simulations indicate, therefore, that the asymptotic
value $w=1/3$, characteristic of radiation, can be attained right after our
simulated period, which could not be reliably foreseen from approximate,
linear analyses.

Figure 6: Evolution of the barotropic parameter $w$ (blue solid line) as a
function of the number of e-folds $N$ from the end of inflation ($\epsilon=1$)
for $n=2$ (left column) and $n=3$ (right column) and $\alpha=10^{-3}$ (top),
$\alpha=10^{-3.5}$ (middle), $\alpha=10^{-4}$ (bottom). Dashed orange vertical
line correspond to $w=1/3$.

Another consequence of geometrical destabilization in models under
consideration, which can be compared with analysis presented in [27], are the
spectra of excited fluctuations of the spectator field $\chi$. Obviously, the
linear approximation used in [27] cannot capture the non-linear effect of
rescattering of produced modes, which leads to broadening of the fluctuations
spectra and production of inflaton fluctuations from those of the spectator.
Hence, such a comparison should be made at the beginning of the geometrical
destabilization. We can see that the initial production of modes is compatible
with the results of the linear theory, and that the spectator spectrum
displays a similar structure as presented in the figure 14 of [27]. Both
spectra are peaked for $k/\mu<1$ as predicted by the linear theory, and
geometrical destabilization initially affects only the modes with
$k/\mu\lesssim 1$ for all considered models.

Figure 7: Evolution of the spectrum of fluctuations $\delta\chi$ of the
spectator field as a function of the number of e-folds $N$ from the end of
inflation for $n=2$ (left column) and $n=3$ (right column) and
$\alpha=10^{-3}$ (top), $\alpha=10^{-3.5}$ (middle), $\alpha=10^{-4}$
(bottom).

##  4 Spectrum of gravitational waves

###  4.1 Numerical calculation of the spectrum

The inhomogeneities of scalar fields produced during preheating can act as a
source of gravitational waves (GWs). A numerical algorithm often used for
computation of spectrum of GWs in numerical lattice simulations was proposed
in [33]. According to [33], the spectrum can be estimated by solving the
following ordinary differential equation in the Fourier space:

${\partial_{\tau}}^{2}\widehat{\bar{h}}_{ij}+\left(|k|^{2}-\frac{{\partial_{\tau}}^{2}a}{a}\right)\widehat{\bar{h}}_{ij}=\frac{2a}{{M_{Pl}}^{2}}\widehat{T^{TT}}_{ij},$
(4.1)

where $\widehat{\bar{h}}_{ij}$ is a Fourier transform of

$\bar{h}_{ij}=ah_{ij},$ (4.2)

which is a transverse ($\partial_{i}h_{ij}=0$) and traceless ($h_{ii}=0$) part
of the fluctuation around Friedman-Robertson-Walker metric:

$\mathrm{d}s^{2}=\mathrm{d}t^{2}-a^{2}(t)\left(\delta_{ij}+h_{ij}\right)\mathrm{d}x^{i}\mathrm{d}x^{j}=a^{2}(\tau)\left(\mathrm{d}\tau^{2}-(\delta_{ij}+h_{ij})\mathrm{d}x^{i}\mathrm{d}x^{j}\right)\,.$
(4.3)

The source term on the right hand side of eq. (4.1) is the Fourier transform
of a transverse-traceless part $T^{TT}$ of the stress tensor $T$. It can be
expressed in a compact form as:

$\widehat{T^{TT}}_{ij}(k)=\mathcal{O}_{ijlm}(k)\widehat{T}_{ij}(k)\,,$ (4.4)

using the projection operators

$\mathcal{O}_{ijlm}(k)=P_{il}(k)P_{jm}(k)-\frac{1}{2}P_{ij}(k)P_{lm}(k)\,,$
(4.5a) where $P_{ij}(k)=\delta_{ij}-\frac{k_{i}k_{j}}{|k|^{2}}\,.$ (4.5b)

In models described by the Lagrangian (2.4), such as T-models of
$\alpha$-attractors, the lowest order contribution to $T^{TT}$ in expansion in
powers of $h$ is given by:

$T_{ij}=e^{2b(\chi)}\partial_{i}\phi\partial_{j}\phi+\partial_{i}\chi\partial_{j}\chi+\operatorname{\mathcal{O}}\left(h\right)+\text{homogeneous
terms}$ (4.6)

According to [33], the energy density of GWs per mode is given by:

$\varrho_{GW}(\tau_{f},k)=\frac{{M_{Pl}}^{2}}{4\pi{a(\tau_{f})}^{4}V}\sum_{i,j}\left(\left|\left(\partial_{\tau}\widehat{\bar{h}}_{ij}\right)(\tau_{f},k)\right|^{2}+|k|^{2}\left|\widehat{\bar{h}}_{ij}(\tau_{f},k)\right|^{2}\right),$
(4.7)

where $\tau_{f}$ is a conformal time at which the source disappears.

In this paper, we are mostly interested in the spectrum of the energy density
of GWs per unit logarithmic frequency interval:

$\frac{d\rho_{GW}}{d\log|k|}(\tau,k)=|k|^{3}\int_{S^{2}}d\Omega\varrho_{gw}(\tau,k),$
(4.8)

where $\int_{S^{2}}d\Omega$ denotes the integration over the direction of the
wave vector $k$.

In [33], the time evolution of the scale factor $a$ was approximated as
proportional to the conformal time, as during radiation domination. In such a
case, the eq. (4.1) can be solved using the retarded Green’s function:

$\widehat{\bar{h}}_{ij}(\tau,k)=\frac{2}{{M_{Pl}}^{2}}\int_{\tau_{i}}^{\tau}d\tau^{\prime}\frac{\sin\left(|k|\left(\tau-\tau^{\prime}\right)\right)}{|k|}a(\tau^{\prime})\widehat{T^{TT}}_{ij}(\tau^{\prime},k),$
(4.9)

where $\tau_{i}$ is the value of conformal time before which the source
$\widehat{T^{TT}}_{ij}$ appeared. However, in our numerical simulations, the
scale factor $a$ is a dynamical variable, whose evolution deviates from simple
linear dependence $a\propto\tau$. In order to take into account this effect,
we developed an algorithm for direct integration of the eq. (4.1). We present
it in the appendix A.

The spectrum of GWs calculated from a lattice simulation needs to be properly
redshifted in order to obtain the spectrum at the present time for comparison
with sensitivities of GWs detectors and exclusion limits. According to [33], a
monochromatic GW with comoving wave-vector $k$ will correspond to wave with
the present-day frequency

$f={a_{f}}^{-1}\left(\frac{k}{\rho_{f}^{1/4}}\right)\left(\frac{a_{f}}{a_{*}}\right)^{1-\frac{3}{4}(1+w)}4\times
10^{10}\,\textrm{Hz},$ (4.10)

where $a_{f}$, $\rho_{f}$ are values of, respectively, the scale factor and
the critical energy density at the time $\tau_{f}$ of end of the simulation
and $a_{*}$ is the value of scale factor at time $\tau_{*}$ when thermal
equilibrium is established. This expression was derived under the assumption
that the Universe evolved from $\tau_{f}$ to $\tau_{*}$ with a mean effective
barotropic parameter equal to $w$. As we showed in the section 3.3, we can
expect that the reheating in most $\alpha$-attractor T-models with considered
parameters is instantaneous and that the effective barotropic parameter
reaches value $\frac{1}{3}$ within one e-fold from the end of the inflation.
Thus, we shall presumably make a negligible error, ignoring the factor
$\frac{a_{f}}{a_{*}}$ in (4.10). The energy density spectrum today is given by

$\Omega_{GW}h^{2}(k)=9.3\times
10^{-6}\rho_{f}^{-4}\frac{d\rho_{GW}}{d\log|k|}(\tau_{f},k).$ (4.11)

Frequencies (4.10) of GWs produced during preheating are generically predicted
to be very large in comparison with sensitivity ranges of ground-based
interferometric detectors (LIGO/Virgo [34, 35, 36, 37] or ET [38, 39]) for GWs
generated around inflation. Although certain ways of detecting GWs at very
high frequencies were proposed [40], their predicted reach in terms of
abundance is still above the current lower bound on the amplitude of GWs
coming from estimation of number of relativistic degrees of freedom
$N_{\text{eff}}$ during Big Bang nucleosynthesis (BBN) based on CMB
observations [41, 42]. Current constraint from the Planck 2018 results,
$N_{\text{eff}}=3.27\pm 0.15$ [43], requires that $\Omega_{GW}h^{2}\leq
1.85\times 10^{-6}$ and we use this bound in the plots of spectra of GWs in
the remainder of our analysis.

###  4.2 Results from numerical calculations

We have computed the spectra of gravitational waves emitted during preheating
in lattice simulations using a numerical algorithm described in the section 4.
We directly integrated in time equation (4.1), using the integrator based on
Magnus expansion presented in the appendix A. Our results for
$\alpha=10^{-3},10^{-3.5},10^{-4}$ are presented in the figure 8 for $n=1,1.5$
and in the figure 9 for $n=2,3$. For more informative character, we plotted
the spectra in various moments of time, i.e. we show how the spectrum of GWs
would look like today if its production out of scalar fluctuations stopped at
given time. Similarly as for spectra of fluctuations of fields $\phi$ and
$\chi$, the color of the line encodes the time. The red band visible at the
top of some figures correspond to BBN exclusion limit discussed in the section
4.1.

Figure 8: Evolution of the spectrum of gravitational waves as a function of
the number of e-folds $N$ from the end of inflation for $n=1$ (left column)
and $n=1.5$ (right column) and $\alpha=10^{-3}$ (top), $\alpha=10^{-3.5}$
(middle), $\alpha=10^{-4}$ (bottom). Lines of different colors correspond to
different time points described by the number of e-folds after the end of
inflation, as indicated on the bar next to each panel. The red band visible at
the top of some figures correspond to BBN exclusion limit discussed in the
section 4.

Figure 9: Evolution of the spectrum of gravitational waves as a function of
the number of e-folds $N$ from the end of inflation for $n=2$ (left column)
and $n=3$ (right column) and $\alpha=10^{-3}$ (top), $\alpha=10^{-3.5}$
(middle), $\alpha=10^{-4}$ (bottom). Lines of different colors correspond to
different time points described by the number of e-folds after the end of
inflation, as indicated on the bar next to each panel. The red band visible at
the top of some figures correspond to BBN exclusion limit discussed in the
section 4.

As we showed in the section 3.3, the process of reheating in two-fields
$\alpha$-attractor T-models has a certain pattern dictated by the geometrical
destabilization of the spectator field. At the end of inflation when
$\epsilon\sim 1$ the negative curvature of the field space produce the
negative contribution to the effective mass of the spectator, driving this
mass tachyonic and leading to the production of spacial fluctuations of this
field. During reheating the energy density of the homogeneous inflaton
background is transferred into excitations of the spectator. Finally, when the
significant fraction of the energy is siphoned off by the spectator, the
nonlinear coupling of fields leads to production of inflaton fluctuations and
rescattering of spectator ones into higher momenta modes. This pattern has its
implications for productions of GWs from fluctuations of fields.

As is clearly visible from plots in figures 8 and 9, the spectra computed for
times before the end of inflation (negative $N$) saturate at extremely low
level of $10^{-24}$—$10^{-20}$ which probably correspond to the numerical
precision with which we reproduce the amplitudes of Bunch–Davis vacuum
fluctuations.

After the end of inflation ($N=0$) the spectra start to grow and develop well
established maxima which are located at the wave vectors $k$ of the order of
two times the location of maxima of spectra of fluctuations of the spectator
field, since the two scalar modes are needed for production of GW. The
amplitudes of the spectra at this stage are roughly of the order of
$10^{-10}$–$10^{-9}$, typical for the resonance driven preheating in inflation
at these scales.

For times around the end of the simulation, high-frequency GWs are massively
produced exceeding the BBN bound for nearly all simulations. This rapid growth
of high energy GWs is associated with production of high-$k$ modes of the
inflaton field due to nonlinear interaction of spectator and inflaton
fluctuations, discussed in the section 3.3. However, there are reasons that
our results for the GWs should be treated as approximate. First, the
underlying instability of the scalar fluctuations cannot be fully resolved in
fixed-lattice simulations, as was previously pointed out in [44]. The
amplitude of the spectrum can be overestimated, since the modes with
wavenumbers larger than the UV cutoff cannot be simulated on the lattice; the
energy that should be deposited in those modes through rescattering stays with
lower-$k$ modes. This effect may affect the computed spectra, boosting their
values, especially at high frequency domain which is responsible for the
violation of the BBN. We looked at this possibility, comparing spectra
computed in simulations with different cutoffs, and found no qualitative
differences. Therefore, we believe that we can treat our results as
trustworthy, especially that we assumed a conservative condition for stopping
our simulations. On one hand, such a conservative condition is meant for
ensuring that our results are reliable during the entire simulated time
period. On the other hand, this leads to a possibility of underestimating the
GW spectra, since the production of GWs will proceed longer than the time span
of our simulation. With all these caveats in mind, we conclude that our
calculations of the spectra GWs point to the possible exclusion of two-field
$\alpha$-attractor T-models due to overproduction of GWs during reheating.

##  5 Conclusions and outlook

In this paper, we extended our earlier investigation of (p)reheating in two-
field $\alpha$-attractor models, reported in [19] in two directions. We
extended our lattice simulations to a larger $512^{3}$ lattice, refining
certain technical aspects to make sure that our results are reliable. We
corroborated the fact that fluctuations of the spectator field are highly
unstable after inflation and quickly grow, with their energy density catching
up with that of the homogeneous inflaton. These fluctuations enter nonlinear
regime and become a driving force for the growth of inflaton fluctuations. We
also find that the growth of the fluctuations is so large and encompasses such
a wide range of modes that it becomes a formidable computational challenge to
resolve this nonlinear dynamics of the lattice.

We also calculated the power spectra of gravitational waves produced out of
the fluctuations of the scalar fields. We found that large amplitudes of the
scalar fields translate into large amplitudes of gravitational waves at
frequencies above the range that is probed or will soon be with ‘direct
detection’ experiments. The energy density of the gravitational waves turned
out to be so large that it is incompatible with constraints on the number of
the extra relativistic degrees of freedom, inferred from the CMB and BBN data.
This would mean that two-field $\alpha$-attractor T-models are ruled out due
to overproduction of gravitational waves during reheating.

The research reported in this paper can be extended in several directions. In
spite of a significant increase of the lattice size and the use of a powerful
computational facility, we have been unable to fully resolve the instability.
This suggests that it is necessary to go beyond the fixed lattice calculations
and to perform multi-scale analysis of the instability. Such a new numerical
tool would allow giving more accurate predictions about the dynamics of the
fluctuations of both the inflaton and the spectator, and about the effective
equation of state of the Universe during post-inflationary geometrical
destabilization. It would also enable a more refined scan over the parameter
space, greatly reducing theoretical uncertainty of the inflationary
predictions for the $\alpha$-attractor T-models. Last but not least, a
possibility of producing primordial black holes out of exponentially growing
curvature perturbations could also be addressed quantitatively. Albeit
interesting, all these goals require a novel numerical approach and
construction of beyond state-of-the-art numerical tools. We shall therefore
consider them a challenge to be addressed in future research.

## Acknowledgments

T.K. is supported by grant 2019/32/C/ST2/00248 from the National Science
Centre (Poland). K.T. is partially supported by grant 2018/30/Q/ST9/00795 from
the National Science Centre (NCN). This research was supported in part by
PL–Grid Infrastructure.

##  Appendix A Numerical scheme for calculation of spectra of gravitational
waves in non-trivial background metric evolution

###  A.1 Construction of the integrator

The method presented here is implementation of the second order scheme based
on Magnus expansion presented in [45].

Eq. (4.1) can be rewritten as a first order non-homogeneous ordinary
differential equation:

$\begin{split}{\partial_{\tau}}Y:=&\,\,\partial_{\tau}\begin{bmatrix}\widehat{\bar{h}}_{ij}\\\
\partial_{\eta}\widehat{\bar{h}}_{ij}\end{bmatrix}=\begin{bmatrix}0&1\\\
-|k|^{2}+\frac{{\partial_{\tau}}^{2}a}{a}&0\end{bmatrix}\begin{bmatrix}\widehat{\bar{h}}_{ij}\\\
\partial_{\eta}\widehat{\bar{h}}_{ij}\end{bmatrix}+\begin{bmatrix}0\\\
\frac{2a}{{M_{Pl}}^{2}}\widehat{T^{TT}}_{ij}\end{bmatrix}\\\
=&\,\,MY+F.\end{split}$ (A.1)

According to [45], a formal solution of (A.1) is of the form:

$Y(\tau+\delta\tau)=\Phi(\tau+\delta\tau,\tau)Y(\tau)+\int_{\tau}^{\tau+\delta\tau}\Phi(\tau+\delta\tau,\xi)F(\xi)d\xi\,,$
(A.2)

where $\Phi$ is the fundamental solution matrix of the associated homogeneous
equation:

$\partial_{\tau}\Phi(\tau,\tau_{0})=M(\tau)\Phi(\tau,\tau_{0}).$ (A.3)

Using Magnus approximation in the lowest order, the fundamental solution of
eq. (A.3) can be represented as:

$\Phi(\tau+\delta\tau,\eta)=\exp{\left(\int_{\tau}^{\tau+\delta\tau}M(\xi)d\xi+\operatorname{\mathcal{O}}\left(\delta\tau^{4}\right)\right)}.$
(A.4)

One can use the trapezoid rule to approximate integrals in eqs. (A.4) and
(A.2) in order to obtain an explicit, time-reversible, second-order
discretization scheme:

$Y_{n+1}=\exp{\left(\frac{\delta\tau}{2}(M_{n+1}+M_{n})\right)}\left(Y_{n}+\frac{\delta\tau}{2}F_{n}\right)+\frac{\delta\tau}{2}F_{n+1}$
(A.5)

given in [45]. Applying (A.5) to the equation of motion of metric
perturbations (A.1) gives the following set of equations:

$\begin{split}\widehat{\bar{h}}_{ij,n+1}=&\cos{(\omega\delta\tau)}\widehat{\bar{h}}_{ij,n}+\frac{i}{\omega}\sin{(\omega\delta\tau)}\left(\partial_{\tau}\widehat{\bar{h}}_{ij,n}+\delta\tau\frac{a_{n}}{{M_{Pl}}^{2}}\widehat{T^{TT}}_{ij,n}\right)\\\
\partial_{\tau}\widehat{\bar{h}}_{ij,n+1}=&\cos{(\omega\delta\tau)}\left(\partial_{\tau}\widehat{\bar{h}}_{ij,n}+\delta\tau\frac{a_{n}}{{M_{Pl}}^{2}}\widehat{T^{TT}}_{ij,n}\right)+i\omega\sin{(\omega\delta\tau)}\widehat{\bar{h}}_{ij,n}\\\
&+\delta\tau\frac{a_{n+1}}{{M_{Pl}}^{2}}\widehat{T^{TT}}_{ij,n+1}\,,\end{split}$
(A.6)

where the frequency $\omega$ is defined as:

$\omega=\sqrt{|k|^{2}-\frac{{\partial_{\tau}}^{2}a_{n+1}}{2a_{n+1}}-\frac{{\partial_{\tau}}^{2}a_{n}}{2a_{n}}}\,,$
(A.7)

and $n$ numerates time steps.

###  A.2 Error estimation for the integrator

In order to use the theory of error estimation for Magnus expansion integrator
developed in [46] one need to rewrite eq. (A.1) as a homogeneous equation as
was described in [45]:

${\partial_{\tau}}Z:=\partial_{\tau}\begin{bmatrix}Y\\\
1\end{bmatrix}=\begin{bmatrix}M&F\\\ 0&0\end{bmatrix}\begin{bmatrix}Y\\\
1\end{bmatrix}=AZ.$ (A.8)

The errors in the integration with discretization scheme given by the eq.
(A.6) come from two sources: from the truncation of the Magnus expansion and
from the trapezoid rule.

We use the error introduced by the first source to adaptatively choose
integration time step, since the estimation of the second requires computation
of time derivatives of the stress tensor. The next term in the Magnus
expansion which we neglected constructing second order scheme is of the form:

$\Omega^{[4]}-\Omega^{[2]}=\left[A^{(1)},A^{(0)}\right],$ (A.9)

where $[\cdot,\cdot]$ is a commutator of matrices and $A^{(i)}$ is defined as:

$A^{(i)}:=\frac{1}{\delta\tau^{i}}\int_{\tau}^{\tau+\delta\tau}\left(\xi-\tau-\frac{\delta\tau}{2}\right)^{i}A(\xi)d\xi.$
(A.10)

In our case, the commutator gives:

$\Omega^{[4]}-\Omega^{[2]}=\begin{bmatrix}-h\left(\frac{{\partial_{\tau}}^{2}a}{a}\right)^{(1)}&0&-h\left(\frac{2a}{{M_{Pl}}^{2}}\widehat{T^{TT}}_{ij}\right)^{(1)}\\\
0&h\left(\frac{{\partial_{\tau}}^{2}a}{a}\right)^{(1)}&0\\\
0&0&0\end{bmatrix}.$ (A.11)

We can use the trapezoidal rule to obtain the lowest order in $\delta\tau$
approximation to the integrals in the error estimation:

$\displaystyle\left(\frac{{\partial_{\tau}}^{2}a}{a}\right)^{(1)}$
$\displaystyle=\frac{\delta\tau}{4}\left(\frac{{\partial_{\tau}}^{2}a_{n+1}}{a^{n+1}}-\frac{{\partial_{\tau}}^{2}a_{n}}{a_{n}}\right)+\operatorname{\mathcal{O}}\left(\delta\tau^{2}\right),$
(A.12)
$\displaystyle\left(\frac{2a}{{M_{Pl}}^{2}}\widehat{T^{TT}}_{ij}\right)^{(1)}$
$\displaystyle=\frac{\delta\tau}{2{M_{Pl}}^{2}}\left(\widehat{T^{TT}}_{ij,n+1}-\widehat{T^{TT}}_{ij,n}\right)+\operatorname{\mathcal{O}}\left(\delta\tau^{2}\right)$
(A.13)

Finally, we need to translate the error estimation for the exponent $\Omega$
into expressions for errors of the components of $Y$ vector. Let us denote
$Z(\tau+\delta\tau)=\exp{(\Omega(\tau+\delta\tau))}Z_{0}$ as an exact solution
to eq. (A.1) and $\widetilde{Z}_{n}=\exp{(\widetilde{\Omega}_{n})}Z_{0}$ as
our approximation calculated using discretization (A.6). Expressing initial
$Z_{0}$ by $\widetilde{Z}_{n}$ with the inverse of the evolution operator, one
gets:

$Z(\tau+\delta\tau)-\widetilde{Z}_{n+1}=\left(\exp{(\Omega(\tau+\delta\tau))}\exp{(-\widetilde{\Omega}_{n})}-1\right)\widetilde{Z}_{n+1}.$
(A.14)

An expansion based on Baker-Campbell-Hausdorff formula leads to:

$Z(\tau+\delta\tau)-\widetilde{Z}_{n+1}=\left(\exp{\left(\Omega^{[4]}_{n}-\Omega^{[2]}_{n}\right)}-1\right)\widetilde{Z}_{n+1}+\operatorname{\mathcal{O}}\left(\delta\tau^{4}\right).$
(A.15)

Due to small dimension of $Z$ the exponent
$\exp{\left(\Omega^{[4]}_{n}-^{[2]}_{n}\right)}$ can be calculated exactly
giving explicit expressions for error estimates, as follows:

$\displaystyle\delta\widehat{\bar{h}}_{ij,n+1}=$
$\displaystyle\left(\exp{\left(-\frac{\delta\tau^{2}}{4}\left(\frac{{\partial_{\tau}}^{2}a_{n+1}}{a_{n+1}}-\frac{{\partial_{\tau}}^{2}a_{n}}{a_{n}}\right)\right)}-1\right)\widehat{\bar{h}}_{ij_{,}n+1}$
(A.16)
$\displaystyle+\frac{\frac{2}{{M_{Pl}}^{2}}\left(\widehat{T^{TT}}_{ij,n+1}-\widehat{T^{TT}}_{ij,n}\right)}{\left(\frac{{\partial_{\tau}}^{2}a_{n+1}}{a_{n+1}}-\frac{{\partial_{\tau}}^{2}a_{n}}{a_{n}}\right)}$
(A.17)
$\displaystyle\left(\exp{\left(-\frac{\delta\tau^{2}}{4}\left(\frac{{\partial_{\tau}}^{2}a_{n+1}}{a_{n+1}}-\frac{{\partial_{\tau}}^{2}a_{n}}{a_{n}}\right)\right)}-1\right)$
(A.18) $\displaystyle\delta\partial_{\tau}\widehat{\bar{h}}_{ij,n+1}=$
$\displaystyle\left(\exp{\left(\frac{\delta\tau^{2}}{4}\left(\frac{{\partial_{\tau}}^{2}a_{n+1}}{a_{n+1}}-\frac{{\partial_{\tau}}^{2}a_{n}}{a_{n}}\right)\right)}-1\right)\partial_{\tau}\widehat{\bar{h}}_{ij,n+1}$
(A.19)

##  Appendix B Supplementary results

In this section, we provide further results from simulations with $n=2$ and
$n=3$ which support our hypothesis, but are qualitatively analogous to cases
$n=1$ and $n=1.5$ discussed in the section 3.3. We believe that they can be
useful for other researchers working in the field.

Figure 10: Evolution of energy components as a function of the number of
e-folds $N$ from the end of inflation for $n=2$ (left column) and $n=3$ (right
column) and $\alpha=10^{-3}$ (top), $\alpha=10^{-3.5}$ (middle),
$\alpha=10^{-4}$ (bottom).

Figure 11: Evolution of the spectrum of fluctuations $\delta\phi$ of the
spectator field as a function of the number of e-folds $N$ from the end of
inflation for $n=2$ (left column) and $n=3$ (right column) and
$\alpha=10^{-3}$ (top), $\alpha=10^{-3.5}$ (middle), $\alpha=10^{-4}$
(bottom).

## References

  * [1] V. Mukhanov, _Physical Foundations of Cosmology_ , Cambridge University Press (2005), 10.1017/CBO9780511790553.
  * [2] J.H. Traschen and R.H. Brandenberger, _Particle Production During Out-of-equilibrium Phase Transitions_ , _Phys. Rev. D_ 42 (1990) 2491.
  * [3] L. Kofman, A.D. Linde and A.A. Starobinsky, _Reheating after inflation_ , _Phys. Rev. Lett._ 73 (1994) 3195 [hep-th/9405187].
  * [4] Y. Shtanov, J.H. Traschen and R.H. Brandenberger, _Universe reheating after inflation_ , _Phys. Rev. D_ 51 (1995) 5438 [hep-ph/9407247].
  * [5] L. Kofman, A.D. Linde and A.A. Starobinsky, _Towards the theory of reheating after inflation_ , _Phys. Rev. D_ 56 (1997) 3258 [hep-ph/9704452].
  * [6] J.J.M. Carrasco, R. Kallosh and A. Linde, _Cosmological Attractors and Initial Conditions for Inflation_ , _Phys. Rev. D_ 92 (2015) 063519 [1506.00936].
  * [7] M.A. Amin, R. Easther, H. Finkel, R. Flauger and M.P. Hertzberg, _Oscillons After Inflation_ , _Phys. Rev. Lett._ 108 (2012) 241302 [1106.3335].
  * [8] K.D. Lozanov and M.A. Amin, _Equation of State and Duration to Radiation Domination after Inflation_ , _Phys. Rev. Lett._ 119 (2017) 061301 [1608.01213].
  * [9] K.D. Lozanov and M.A. Amin, _Self-resonance after inflation: oscillons, transients and radiation domination_ , _Phys. Rev. D_ 97 (2018) 023533 [1710.06851].
  * [10] 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].
  * [11] P. Adshead, R. Easther, J. Pritchard and A. Loeb, _Inflation and the Scale Dependent Spectral Index: Prospects and Strategies_ , _JCAP_ 02 (2011) 021 [1007.3748].
  * [12] P. Creminelli, D. López Nacir, M. Simonović, G. Trevisan and M. Zaldarriaga, _$\phi^{2}$ or Not $\phi^{2}$: Testing the Simplest Inflationary Potential_, _Phys. Rev. Lett._ 112 (2014) 241303 [1404.1065].
  * [13] L. Dai, M. Kamionkowski and J. Wang, _Reheating constraints to inflationary models_ , _Phys. Rev. Lett._ 113 (2014) 041302 [1404.6704].
  * [14] J. Martin, C. Ringeval and V. Vennin, _Observing Inflationary Reheating_ , _Phys. Rev. Lett._ 114 (2015) 081303 [1410.7958].
  * [15] J.L. Cook, E. Dimastrogiovanni, D.A. Easson and L.M. Krauss, _Reheating predictions in single field inflation_ , _JCAP_ 04 (2015) 047 [1502.04673].
  * [16] Y. Ueno and K. Yamamoto, _Constraints on $\alpha$-attractor inflation and reheating_, _Phys. Rev. D_ 93 (2016) 083524 [1602.07427].
  * [17] M. Eshaghi, M. Zarei, N. Riazi and A. Kiasatpour, _CMB and reheating constraints to $\alpha$-attractor inflationary models_, _Phys. Rev. D_ 93 (2016) 123517 [1602.07914].
  * [18] S. Renaux-Petel and K. Turzyński, _Geometrical Destabilization of Inflation_ , _Phys. Rev. Lett._ 117 (2016) 141301 [1510.01281].
  * [19] T. Krajewski, K. Turzyński and M. Wieczorek, _On preheating in $\alpha$-attractor models of inflation_, _Eur. Phys. J. C_ 79 (2019) 654 [1801.01786].
  * [20] S. Ferrara, R. Kallosh, A. Linde and M. Porrati, _Minimal Supergravity Models of Inflation_ , _Phys. Rev. D_ 88 (2013) 085038 [1307.7696].
  * [21] R. Kallosh and A. Linde, _Non-minimal Inflationary Attractors_ , _JCAP_ 10 (2013) 033 [1307.7938].
  * [22] R. Kallosh and A. Linde, _Multi-field Conformal Cosmological Attractors_ , _JCAP_ 12 (2013) 006 [1309.2015].
  * [23] R. Kallosh, A. Linde and D. Roest, _Superconformal Inflationary $\alpha$-Attractors_, _JHEP_ 11 (2013) 198 [1311.0472].
  * [24] R. Kallosh and A. Linde, _Planck, LHC, and $\alpha$-attractors_, _Phys. Rev. D_ 91 (2015) 083528 [1502.07733].
  * [25] J.J.M. Carrasco, R. Kallosh, A. Linde and D. Roest, _Hyperbolic geometry of cosmological attractors_ , _Phys. Rev. D_ 92 (2015) 041301 [1504.05557].
  * [26] A. Achúcarro, R. Kallosh, A. Linde, D.-G. Wang and Y. Welling, _Universality of multi-field $\alpha$-attractors_, _JCAP_ 04 (2018) 028 [1711.09478].
  * [27] O. Iarygina, E.I. Sfakianakis, D.-G. Wang and A. Achucarro, _Universality and scaling in multi-field $\alpha$-attractor preheating_, _JCAP_ 06 (2019) 027 [1810.02804].
  * [28] Z. Lalak, D. Langlois, S. Pokorski and K. Turzynski, _Curvature and isocurvature perturbations in two-field inflation_ , _JCAP_ 07 (2007) 014 [0704.0212].
  * [29] M.A. Amin, M.P. Hertzberg, D.I. Kaiser and J. Karouby, _Nonperturbative Dynamics Of Reheating After Inflation: A Review_ , _Int. J. Mod. Phys. D_ 24 (2014) 1530003 [1410.3808].
  * [30] D. Langlois and S. Renaux-Petel, _Perturbations in generalized multi-field inflation_ , _JCAP_ 04 (2008) 017 [0801.1085].
  * [31] V. Mukhanov and S. Winitzki, _Introduction to Quantum Effects in Gravity_ , Cambridge University Press (2007), 10.1017/CBO9780511809149.
  * [32] G.N. Felder and I. Tkachev, _LATTICEEASY: A Program for lattice simulations of scalar fields in an expanding universe_ , _Comput. Phys. Commun._ 178 (2008) 929 [hep-ph/0011159].
  * [33] J.F. Dufaux, A. Bergman, G.N. Felder, L. Kofman and J.-P. Uzan, _Theory and Numerics of Gravitational Waves from Preheating after Inflation_ , _Phys. Rev. D_ 76 (2007) 123517 [0707.0875].
  * [34] E. Thrane and J.D. Romano, _Sensitivity curves for searches for gravitational-wave backgrounds_ , _Phys. Rev._ D88 (2013) 124032 [1310.5300].
  * [35] LIGO Scientific collaboration, _Advanced LIGO_ , _Class. Quant. Grav._ 32 (2015) 074001 [1411.4547].
  * [36] LIGO Scientific, Virgo collaboration, _GW150914: Implications for the stochastic gravitational wave background from binary black holes_ , _Phys. Rev. Lett._ 116 (2016) 131102 [1602.03847].
  * [37] LIGO Scientific, Virgo collaboration, _Search for the isotropic stochastic background using data from Advanced LIGO’s second observing run_ , _Phys. Rev. D_ 100 (2019) 061101 [1903.02886].
  * [38] M. Punturo et al., _The Einstein Telescope: A third-generation gravitational wave observatory_ , _Class. Quant. Grav._ 27 (2010) 194002.
  * [39] S. Hild et al., _Sensitivity Studies for Third-Generation Gravitational Wave Observatories_ , _Class. Quant. Grav._ 28 (2011) 094013 [1012.0908].
  * [40] N. Aggarwal et al., _Challenges and Opportunities of Gravitational Wave Searches at MHz to GHz Frequencies_ , 2011.12414.
  * [41] T.L. Smith, E. Pierpaoli and M. Kamionkowski, _A new cosmic microwave background constraint to primordial gravitational waves_ , _Phys. Rev. Lett._ 97 (2006) 021301 [astro-ph/0603144].
  * [42] S. Henrot-Versille et al., _Improved constraint on the primordial gravitational-wave density using recent cosmological data and its impact on cosmic string models_ , _Class. Quant. Grav._ 32 (2015) 045003 [1408.5299].
  * [43] Planck collaboration, _Planck 2018 results. VI. Cosmological parameters_ , _Astron. Astrophys._ 641 (2020) A6 [1807.06209].
  * [44] J. Li, H. Yu and P. Wu, _Production of gravitational waves during preheating in $\alpha$-attractor inflation_, _Phys. Rev. D_ 102 (2020) 083522.
  * [45] S. Blanes and E. Ponsoda, _Time-averaging and exponential integrators for non-homogeneous linear ivps and bvps_ , _Applied Numerical Mathematics_ 62 (2012) 875.
  * [46] A. Iserles, A. Marthinsen and S. Norsett, _On the implementation of the method of magnus series for linear differential equations_ , _BIT_ 39 (1999) 281.