# Anatomy of geometrical destabilization of inflation

Tomasz Krajewski  Krzysztof Turzyński

###### Abstract

We study geometrical destabilization of inflation with the aim of determining
the fate of excited unstable modes. We use numerical lattice simulations to
track the dynamics of both the inflaton and spectator field. We find that
geometrical destabilization is a short-lived phenomenon and that a negative
feedback loop prevents field fluctuations from growing indefinitely. As a
result, fields undergoing geometrical destabilization are merely shifted to a
new classical configuration corresponding to a uniform value of the spectator
field within a Hubble patch.

##  1 Introduction

Cosmological inflation [1, 2, 3, 4, 5, 6] has by now become a natural
ingredient of the standard cosmological model (see e.g. [7] for a pedagogical
introduction). It describes a phase of the evolution of the Universe in which
its expansion accelerates and vacuum quantum fluctuations of the gravitational
and matter fields are amplified to cosmological perturbations [8, 9, 10, 11,
12, 13], which later seed the cosmic microwave background (CMB) anisotropies
and the large-scale structure of our Universe.

As inflation is rather a broad theoretical framework than a well-defined
theory, there are many open questions about specific points of history of the
inflationary Universe. One of them is related to reheating, i.e. the passage
from the inflationary era to radiation-dominated era, because there is no
definitive prescription for coupling the inflationary sector to the Standard
Model of particle physics or its extension thereof.

Nonetheless, already the pioneering works of Starobinsky describe the model of
inflation together with reheating by gravitational particle creation in the
regime of the weak narrow parametric resonance [14, 15] (see, e.g., [16] for a
review). At present, the observations can be effectively described with a
minimal setup: a single scalar inflaton field with canonical kinetic term,
minimally coupled to gravity, and evolving in a sufficiently flat potential
[17, 18, 19, 20, 21]. Inflation is then typically terminated when the inflaton
leaves the slow-roll regime. The subsequent advent of the radiation domination
era is not related do any distinct observational signatures, which has led
many authors to include the ignorance about that era into theoretical
uncertainties in the predictions of inflationary models.

However, many inflationary models come naturally equipped with the possibility
that an effective classical force associated with the inflaton can lead to
non-adiabatic excitations of scalar field fluctuations through parametric
resonance, thereby modifying the evolution of the Universe. This view was
first advocated (for narrow parametric resonance) in [22], but the effect was
later shown negligible [23], whereas the amplification of those fluctuations
by a broad parametric resonance [23, 24, 25] is still a viable candidate for
reheating.

This phenomenon of preheating is typically considered after the slow-roll
phase (see, e.g., [26, 27] for a review), but it can also affect the dynamics
of the Universe during inflation in the multi-field setup [28], including a
premature end of inflation. There are also other mechanisms which can end
inflation prematurely or change its course. In hybrid inflation [29, 30, 31,
32], there is a ‘spectator’ field with an inflaton-dependent mass. This field
develops a tachyonic instability which either quickly terminates inflation or
the field takes over the role of the original inflaton field [33, 34, 35].

In the context of multi-field inflationary models with non-canonical kinetic
terms, a possibility of geometrical destabilization has been proposed [36].
With a negative field-space curvature, an inflationary trajectory can be
destabilized by a ‘geometrical’ force dominating over forces originating from
the potential. It has been hinted that geometrical destabilization may either
end inflation prematurely [37] or trigger a new phase of inflation [39, 38],
especially that geometrical destabilization is a self-constraining phenomenon
[40].

The goal of our analysis is to provide a definitive answer to this dilemma by
means of state-of-the-art numerical lattice simulations of the dynamics of the
scalar fields.

Our work is organized as follows. In Section 2, we briefly describe the basics
of geometrical destabilization. In Section 3 we introduce specific,
representative models of geometrical destabilization and present results of
lattice simulations of scalar field dynamics in these models. The discussion
of our numerical results is presented in Section 4. Technical remarks about
our numerical lattice simulations are relegated to Appendix A. Throughout the
paper we adopt natural units with Planck mass $M_{P}=1$, unless indicated
otherwise.

##  2 Rudiments of geometrical destabilization

The idea of geometrical destabilization is based on the possibility that the
action for scalar fields minimally coupled to gravity contains kinetic terms
that can be written as:

${\cal L}_{\rm
kin}=-\frac{1}{2}G_{IJ}(\phi^{K})g^{\mu\nu}\partial_{\mu}\phi^{I}\partial_{\nu}\phi^{J}\,,$
(2.1)

where the manifold described by the field-space metric $G_{IJ}$ has nonzero
curvature.

The dynamics of models of inflation described by the Lagrangian (2.1) have
been extensively studied in the past two decades (see, e.g., [41, 42, 44]). On
a spatially flat Friedmann-Lemaître-Robertson-Walker universe, with metric

$\mathrm{d}s^{2}=-\mathrm{d}t^{2}+a^{2}(t)\mathrm{d}{\vec{x}}^{2}\,,$ (2.2)

where $t$ is cosmic time and $a(t)$ denotes the scale factor, and with
homogeneous scalar fields $\phi^{I}$, the equations of motion take the form:

$\displaystyle 3H^{2}M_{Pl}^{2}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\dot{\sigma}^{2}+V\,,$ (2.3)
$\displaystyle\dot{H}M_{P}^{2}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\dot{\sigma}^{2}\,,$ (2.4) $\displaystyle{\cal
D}_{t}\dot{\phi}^{I}+3H\dot{\phi}^{I}+G^{IJ}V_{,J}$ $\displaystyle=$
$\displaystyle 0\,.$ (2.5)

In these expressions, dots denote derivatives with respect to $t$,
$H\equiv\dot{a}/a$ is the Hubble parameter,
$\tfrac{1}{2}\dot{\sigma}^{2}\equiv\tfrac{1}{2}G_{IJ}\dot{\phi}^{I}\dot{\phi}^{J}$
is the kinetic energy of the fields, and, hereafter, ${\cal
D}_{t}A^{I}\equiv\dot{A^{I}}+\Gamma^{I}_{JK}\dot{\phi}^{J}A^{K}$ for a field-
space vector $A^{I}$ (field-space indices are lowered and raised with the
field-space metric and its inverse respectively).

The behavior of linear fluctuations about such a background is described by
the second-order action

$\displaystyle
S_{(2)}=\int\mathrm{d}t\,\mathrm{d}^{3}x\,a^{3}\left(G_{IJ}\mathcal{D}_{t}Q^{I}\mathcal{D}_{t}Q^{J}-\frac{1}{a^{2}}G_{IJ}\partial_{i}Q^{I}\partial^{i}Q^{J}-M_{IJ}Q^{I}Q^{J}\right),$
(2.6)

where $Q^{I}$’s are the Mukhanov-Sasaki variables and $M_{IJ}$ is a mass
(squared) matrix. The equations of motion following from (2.6) read:

$\displaystyle{\cal D}_{t}{\cal D}_{t}Q^{I}+3H{\cal
D}_{t}Q^{I}+\frac{k^{2}}{a^{2}}Q^{I}+M^{I}_{\,J}Q^{J}=0\,$ (2.7)

with

$\displaystyle M^{I}_{\,J}$ $\displaystyle=$ $\displaystyle
V^{I}_{;J}-{\mathcal{R}^{I}_{\,KLJ}\dot{\phi}^{K}\dot{\phi}^{L}}-\frac{1}{a^{3}M_{Pl}^{2}}\mathcal{D}_{t}\left(\frac{a^{3}}{H}\dot{\phi}^{I}\dot{\phi}_{J}\right).$
(2.8)

In a curved field space, i.e. for non-vanishing $\mathcal{R}^{I}_{\,KLJ}$, the
second term in (2.8) can be negative. As a result, at least some of the field
fluctuactions can become tachyonic, making the inflationary trajectory
unstable. More specifically, we can rewrite equations of motion (2.7) in the
adiabatic-entropic basis $(e_{\sigma}^{I},e_{s}^{I})$ [43, 44], where
$e_{\sigma}^{I}\equiv\dot{\phi}^{I}/{\dot{\sigma}}$ is the unit vector tangent
to the background trajectory in field space, and $e_{s}^{I}$ is orthonormal to
$e_{\sigma}^{I}$. The adiabatic perturbation $Q_{\sigma}\equiv e_{\sigma
I}Q^{I}$ is proportional to the comoving curvature perturbation ${\cal
R}=\frac{H}{{\dot{\sigma}}}Q_{\sigma}$, while the entropic fluctuation
$Q_{s}$, perpendicular to the background trajectory, exhibits genuine multi-
field effects.

The equation of motion for superhorizon modes of the entropic fluctuation
simplifies as well to

$\ddot{Q}_{s}+3H\dot{Q}_{s}+m^{2}_{s{\rm(eff)}}Q_{s}=0\,,$ (2.9)

where we denote by $m^{2}_{s{\rm(eff)}}$ the effective entropic mass on
superhorizon scales

$\frac{m^{2}_{s{\rm(eff)}}}{H^{2}}\equiv\frac{V_{;ss}}{H^{2}}+3\eta_{\perp}^{2}+\epsilon_{1}\,\mathbb{R}M_{Pl}^{2}\,.$
(2.10)

It contains three contributions: the first one is the usual Hessian of the
potential, the second reflects bending of the inflationary trajectory and the
third, proportional to the field-space Ricci scalar curvature $\mathbb{R}$,
encodes field-space geometrical effects [45, 46].

From (2.9) and (2.10), the mechanism of geometrical destabilization of
inflation is readily identified: it corresponds to situations in which the
geometrical contribution is negative and dominates the sum of the two other
contributions, so that the entropic fluctuation is tachyonic, and the
underlying background trajectory is unstable. As $\epsilon_{1}$ is a positive
quantity, the geometrical destabilization in two-field models can only arise
in setups with a scalar curvature that is negative, which is related to the
fact that this makes neighbouring geodesics diverge from one another. When the
field-space curvature is positive, it renders the entropic fluctuations even
more massive, and does not modify the standard picture. Hence, from now in, we
shall consider only negatively curved field spaces.

##  3 Results

###  3.1 A minimal realization of geometrical destabilization

For our numerical simulations, we use a model of slow-roll inflation driven by
a scalar field $\phi$ with canonical kinetic term and potential $V(\phi)$,
with the Lagrangian ${\cal L}_{\phi}=-\frac{1}{2}(\partial\phi)^{2}-V(\phi)$,
where

$V(\phi)=\Lambda^{4}\left(1-e^{-\sqrt{\frac{2}{3}}\frac{\phi}{M_{Pl}}}\right)^{2}$
(3.1)

is the Starobinsky potential normalized to reproduce the central value of the
observed normalization of the curvature perturbations. In addition, we
consider an extra scalar field $\chi$, stabilized at the bottom of its
potential by a large mass, larger than the Hubble scale. This is described by
the simple Lagrangian ${\cal
L}_{\chi}=-\frac{1}{2}(\partial\chi)^{2}-\frac{1}{2}m^{2}\chi^{2}$, where $m$
stands for the heavy mass, i.e. $m^{2}\gg H^{2}$.

The dimension six operator describing interactions with the two sector
outlined above is ${\cal L}_{{\rm
int}}\propto-(\partial\phi)^{2}\chi^{2}/M^{2}$, where $M$ is a scale of new
physics that lies well above the Hubble scale, $M\gg H$. Such an operator
respects the (approximate) shift-symmetry of the inflaton and is therefore
expected from an effective field theory point of view. Our total Lagrangian
thus reads

${\cal
L}=-\frac{1}{2}(\partial\phi)^{2}\left(1+2\frac{\chi^{2}}{M^{2}}\right)-V(\phi)-\frac{1}{2}(\partial\chi)^{2}-\frac{1}{2}m^{2}\chi^{2}\,.$
(3.2)

The dimension six operator generates a curved field space with metric
$(1+2\chi^{2}/M^{2})({\rm d}\phi)^{2}+({\rm d}\chi)^{2}$, whose Ricci scalar
is negative and reads $\mathbb{R}=-4/M^{2}\times(1+2\chi^{2}/M^{2})^{-2}$.
Along the inflationary valley $\chi=0$, the entropic fluctuation $Q_{s}$,
which then simply coincides with the fluctuation of $\chi$, thus acquires the
effective mass (2.10), i.e.

$m^{2}_{s{\rm(eff)}}=m^{2}-4\epsilon_{1}H^{2}\left(M_{Pl}/M\right)^{2}\,,$
(3.3)

as we have here $V_{;ss}=m^{2},\left(\mathbb{R}\right)|_{\chi=0}=-4/M^{2}$,
and the inflationary trajectory along $\chi=0$ is a field-space geodesic, so
that $\eta_{\perp}=0$.

$\frac{m_{h}^{2}}{H_{\rm c}^{2}}$ | $\frac{M^{2}}{M_{Pl}^{2}}$ | $\frac{k_{\rm cutoff}}{M_{Pl}}$ | $\frac{m_{h}^{2}}{H_{\rm c}^{2}}$ | $\frac{M^{2}}{M_{Pl}^{2}}$ | $\frac{k_{\rm cutoff}}{M_{Pl}}$
---|---|---|---|---|---
Model 1A | Model 1B
$10^{2}$ | $10^{-3}$ | $320$ | $10^{2}$ | $10^{-4}$ | $320$
Model 2A | Model 2B
$10^{3}$ | $10^{-4}$ | $80$ | $10^{3}$ | $10^{-5}$ | $80$
Model 3A | Model 3B
$10^{4}$ | $10^{-5}$ | $40$ | $10^{4}$ | $10^{-6}$ | $40$
$\epsilon_{1,{\rm c}}=2.5\cdot 10^{-2}$ | $\epsilon_{1,{\rm c}}=2.5\cdot 10^{-3}$
$\Delta N_{\rm c}=4.0$ | $\Delta N_{\rm c}=15.4$
Table 1: Six parameter sets used in the lattice simulations. Models denoted by
different numbers correspond to different mass parameters of $\chi$. Models
denoted by different letters correspond to different values of
$\epsilon_{1,{\rm c}}$, i.e., different onsets of geometrical destabilization,
and to different values of $\Delta N_{\rm c}$, i.e., different numbers of
e-folds between the onset of geometrical destabilization and the end of
single-field inflation.

As $\epsilon_{1}H^{2}$ grows during inflation, at the critical point such that

$\epsilon_{1,{\rm c}}=\frac{1}{4}\left(\frac{m}{H_{\rm
c}}\right)^{2}\left(\frac{M}{M_{Pl}}\right)^{2}\,.$ (3.4)

the effective mass (3.3) becomes negative, which triggers geometrical
destabilization of $\chi$. The model has thus two tunable dimensionless
parameters, $m^{2}/H_{\rm c}^{2}$ and $M^{2}/M_{Pl}^{2}$. We shall consider
six sets of parameters, listed in Table 1. For further convenience, the layout
of Table 1 directly corresponds to the layout of subsequent Figures.

###  3.2 Presentation of the lattice simulations results

Figure 1: Evolution of various components of the energy density for six models
listed in Table 1. There are six components, corresponding to the kinetic,
gradient and potential energy of each of fields $\phi$ and $\chi$.

Perhaps the most instructive way of presenting our results consists in showing
the time evolution of various components of the total energy density, as
illustrated in the Figure 1. We show kinetic and gradient energy density of
each of fields $\phi$ and $\chi$, as well as the total potential energy
density and the potential energy density for $\chi$. The separation of the
kinetic and gradient energy density into components corresponding to each
field is possible because the metric of the field space is diagonal. The
separation of the potential energy density into components corresponding to
each field could be done, because the potential is a sum of two independent
contributions, one from each field. The onset of geometrical destabilization,
corresponding to the moment at which the effective mass (3.3) changes sign to
negative is marked as $N=0$.

We observe that quite soon after geometrical destabilization begins, the
kinetic and potential energy density of $\chi$ starts to grow, and that this
growth is soon followed, but not matched, by gradient energy densities of the
fields. This shows that spatial inhomogeneities have a subdominant role in the
evolution of the Universe. The growth of tachyonic instabilities is very short
and it is quickly terminated by backreaction effects.

Figure 2: Evolution of the slow roll parameter $\epsilon_{1}$ for six models
listed in Table 1. Solid blue lines correspond to the results of lattice
simulations. Dashed orange lines correspond to the approximate analysis
presented in Ref. [40].

In Ref. [40], it was hypothesized that the backreaction was mainly caused by
the growth of amplitudes of perturbations of the field $\chi$, which led to a
geometrical suppression of the curvature term in the effective mass (3.3).
This conclusion was obtained by following equations of motion for homogeneous
fields, with initial conditions with the displacement of the field strength
$\chi$ from zero to a value inferred from stochastic inflation considerations.

Our numerical studies do not confirm this picture on the quantitative level.
In Figure 2 we show the evolution of the slow-roll parameter $\epsilon_{1}$,
superimposing the results of our lattice simulations (solid blue lines) with
the approximate method of [40]. It is clearly visible that the backreaction
effects become important much earlier than the approximate analysis suggests
and this significantly reduces the effects of geometrical destabilization.

Figure 3: Evolution of the standard deviation of $\chi$ and for six models
listed in Table 1. Dashed orange lines correspond to the results for the
average value of $\chi$ obtained from the approximate analysis presented in
Ref. [40].

Similar conclusions can be drawn from the evolution of $\chi$. In Figure 3, we
present standard deviation of this field strength, and compare it with the
results obtained with the approximate method of [40]. Comparing the standard
deviation with the average field value, we note that in all cases the former
is much larger than the latter, which may suggest the formation of patches in
space, with values of $\chi$ differing in sign. Additionally, the approximate
method underestimates the fluctuations growth at the early stages of the
geometrical destabilization and overestimates the range of oscillations of
those fluctuations.

In order to corroborate the formation of patches in the Universe,
corresponding to approximately uniform values of the displaced $\chi$ field,
in Figure 4 we show the distribution of the values of the spectator field
$\chi$ and the histograms of relative probability density of finding $\chi$
close to a particular value for four different instances, corresponding to
times before, during and after geometrical destabilization. We can see that an
initial random distribution is transformed to one that exhibits discernible
patches corresponding to different signs but the same overall amplitude of
$\chi$, and separated by domain walls. We can also determine that the size of
these patches roughly corresponds to the Hubble radius at the onset of
geometrical destabilization.111This effect was first reported based on a much
smaller-scale simulations in [47].

Figure 4: Left column: snapshots of the spatial distribution of the spectator
field $\chi$ before (upper panel), during (both middle panels) and after
(lower panel) geometrical destabilization on a section of the lattice. Right
column: histograms of relative frequency of different values of the spectator
field $\chi$ on a section of the lattice. All results are given for model 1A
listed in Table 1.

In order to analyze further the formation of patches in the Universe,
corresponding to approximately uniform values of the displaced $\chi$ field,
we studied how the distribution of the amplitudes of $\chi$ in nodes of the
lattice evolves in time. The results are presented in Figure 5. We can see
that before the onset of geometrical destabilization the most probable field
value is peaked at zero; later, the evolution of the field bifurcates and the
distribution becomes bimodal, with two sharp peaks corresponding to an almost
constant field value within each patch, and that neighboring patches differ in
the sign of $\chi$, while the overall amplitude of the field is the same.

Figure 5: Time evolution of the distribution of the amplitude of $\chi$ in
nodes of the lattice for six models listed in Table 1. Field values from the
range displayed in the plots are binned and the shade of the bin corresponds
to the proportion of nodes at which the field value correspond to a given bin.

##  4 Discussion and conclusions

When geometrical destabilization of inflation was originally advocated [36],
it was hypothesized that there could be two possible outcomes of this
instability. Either the curvature perturbations could grow uncontrollably,
effectively terminating inflation, or their evolution could be much less
dramatic, with a shift of the spectator field to a new value determined by a
balance between the divergence of the geodesic lines (producing an apparent
force akin to apparent forces in a non-inertial reference frame) and the force
resulting from the field potential.

The first option was adopted in [37], which provided a Bayesian analysis of
the parameter space of the inflationary models under the assumption of a
sudden termination of inflation.

A radically different view was put forth in [38], which assumed that the
spectator field is merely ‘side-kicked’ to a new classical state and that
after geometrical destabilization inflation proceeds along a new classical
trajectory. Inflation along this trajectory can give rise to a number of
interesting phenomena, e.g. a transient instability of the curvature
perturbations related to a negative effective sound speed [48] (see also [49,
50]) and a peculiar shape and amplitude of generated non-Gaussianities [38].
Similarly, it was suggested in [39] that geometrical destabilization of
inflation must be terminated, because the would-be spectator mass evaluated on
the attractor solution is always positive.

The arguments recapitulated above are based on the analysis of a homogeneous
inflationary background, so it is unclear whether they can be applicable once
geometrical destabilization is in place. A step further was made in [40],
showing a negative feedback loop intrinsically built into any model of
geometrical destabilization. What remained elusive at that point was a fully
numerical corroboration of that observation.

Our work fills what we believe was the last gap in the discussion of
geometrical destabilization of inflation. Here we have performed state-of-the-
art lattice simulations to show that that the instability is quickly shut off
by backreaction effects and that they are even stronger than previously
anticipated. Our calculations have conclusively shown that geometrical
destabilization is a short-lasting phenomenon and that its main consequence
consists in the fact that the classical fields leave the configuration which
is no longer stable and move, within a causally connected patch, to a new,
stable one.

#### Acknowledgements

We thank Michał Wieczorek for collaboration at an early stage of this project.
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. T.K. acknowledges hospitality of Rudolf Peierls Centre
for Theoretical Physics at Oxford University, where parts of this work has
been done. K.T. thanks International Physicists’ Tournament at Universidad
Industrial de Santander in Bucaramanga (Colombia) for fruitful discussions and
stimulating atmosphere during the last stages of the completion of this work.

##  Appendix A Numerical discretization scheme for lattice simulations

The method presented here is a slight modification of the method that we used
for studies of dynamics of preheating in $\alpha$-attractor T-models of
inflation in [51, 52]. Only non-canonical kinetic term coupling $e^{2b(\chi)}$
and the potential $V(\phi,\chi)$ had to be modified. For the model under
consideration these functions are:

$\displaystyle e^{2b(\chi)}$
$\displaystyle=1+\frac{2\chi^{2}}{M_{\text{np}}^{2}},$ $\displaystyle
V(\phi,\chi)$ $\displaystyle=V(\phi)+\frac{1}{2}m_{h}^{2}\chi^{2}.$ (A.1)

In the current section we will present the algorithm in its general
formulation.

###  A.1 Lattice formulation of the model

We consider the general model described by the following action

$S=\int
d^{4}x\sqrt{-g}\bigg{[}\frac{{M_{Pl}}^{2}}{2}\mathcal{R}-\frac{1}{2}e^{2b(\chi)}(\partial_{\mu}\phi)(\partial^{\mu}\phi)-\frac{1}{2}(\partial_{\mu}\chi)(\partial^{\mu}\chi)-V(\phi,\chi)\bigg{]}.$
(A.2)

We assume here that the spacetime is spatially homogeneous, isotropic and
flat, thus we approximate the metric tensor field by Friedman-Lemaître-
Robertson-Walker metric expressed in terms of the conformal time $\tau$:

$ds^{2}=a^{2}(-d\tau^{2}+d\mathbf{x}^{2}),$ (A.3)

which implies

$\sqrt{-g}=a^{4}\qquad\textrm{and}\qquad\mathcal{R}=6\frac{a^{\prime\prime}}{a^{3}},$
(A.4)

where the prime denotes the derivative with respect to $\tau$. In our
simulations we neglected the backreaction from metric fluctuations on the
evolution of both the inflaton and the spectator, still keeping tack of the
evolution of the scale factor $a$ of the background metric.

Our aim is to obtain symplectic integrator for model under consideration. We
use the method of lines, i.e. first we discretized the action in space and
then we constructed the symplectic integrator for the theory defined on the
lattice. We simulate some patch of the Universe with finite comoving volume
$U$, thus discrete theory has only finite number of degrees of freedom. Under
the assumption of homogeneity of the Universe it is natural to use periodic
boundary conditions.

We use values of the field strengths of $\chi$ and $\phi$ at the nodes of
cubic regular lattice as degrees of freedom for discrete theory. The spacial
gradients of fields can be approximated using finite difference method. In our
implementation we used first order forward finite difference scheme. This
choice is equivalent to linear interpolation of field strength values at cells
of the lattice.

After discretization in space the action (A.2) can be written as:

$\displaystyle S$ $\displaystyle=$ $\displaystyle\int Ld\tau=$
$\displaystyle=$ $\displaystyle\int\bigg{[}-3a^{\prime
2}U{M_{Pl}}^{2}+\sum_{\vec{x}}\frac{a^{2}}{2}\frac{U}{V_{L}}\Bigg{(}e^{2b(\chi_{\vec{x}})}\bigg{(}(\phi^{\prime}_{\vec{x}})^{2}-G_{\vec{x}}(\phi)\bigg{)}+$
$\displaystyle+$
$\displaystyle\bigg{(}(\chi^{\prime}_{\vec{x}})^{2}-G_{\vec{x}}(\chi)\bigg{)}-a^{2}V(\phi_{\vec{x}},\chi_{\vec{x}})\Bigg{)}\bigg{]}d\tau,$

where $V_{L}$ is the number of cells in the lattice, thus $\frac{U}{V_{L}}$ is
equal to the volume of each cell, $\vec{x}=(x_{1},x_{2},x_{3})^{T}\colon
x_{1},x_{2},x_{3}\in\mathbb{N}$ is a multi-index that numerates lattice points
in three dimensions and

$G_{\vec{x}}(Y)=\frac{1}{2\delta^{2}}\sum_{i=1}^{3}(Y_{\vec{x}+\vec{e_{i}}}-Y_{\vec{x}})^{2}$
(A.6)

is the square of discretization of the spatial gradient with
$\vec{e_{1}}=(1,0,0)^{T}$, $\vec{e_{2}}=(0,1,0)^{T}$ and
$\vec{e_{3}}=(0,0,1)^{T}$ and $\delta$ is comoving length of the edge of a
lattice cell.

After droping common factor $\frac{U}{V_{L}}$ and the Legendre transformation
we obtain the following Hamiltonian:

$H=-\frac{p_{a}^{2}}{12V_{L}{M_{Pl}}^{2}}+\sum_{\vec{x}}a^{4}\bigg{(}\frac{\pi^{2}_{\phi,\vec{x}}}{2a^{6}e^{2b(\chi_{\vec{x}})}}+\frac{\pi^{2}_{\chi,\vec{x}}}{2a^{6}}+e^{2b(\chi_{\vec{x}})}\frac{G_{\vec{x}}(\phi)}{2\delta^{2}a^{2}}+\frac{G_{\vec{x}}(\chi)}{2\delta^{2}a^{2}}+V(\phi_{\vec{x}},\chi_{\vec{x}})\bigg{)},$
(A.7)

where canonical momenta are defined by formulae:

$p_{a}\equiv\frac{\partial L}{\partial
a^{\prime}}=-6a^{\prime}V_{L}{M_{Pl}}^{2},\quad\pi_{\phi,\vec{x}}\equiv\frac{\partial
L}{\partial\phi^{\prime}}=a^{2}\textrm{e}^{2b(\chi_{\vec{x}})}\phi^{\prime}_{\vec{x}}\quad\textrm{and}\quad\pi_{\chi,\vec{x}}\equiv\frac{\partial
L}{\partial\chi^{\prime}}=a^{2}\chi^{\prime}_{\vec{x}}.$ (A.8)

From the point of view of implementation of numerical algorithm it is
convenient to use

$\wp_{a}:=\frac{p_{a}}{V_{L}}=-6a^{\prime}{M_{Pl}}^{2}.$ (A.9)

instead of canonical momentum $p_{a}$.

###  A.2 Time integration scheme

Construction of the time integrator that we used in our simulations is based
on the technique of operator splitting. The hamiltonian $H$ can be divided
into four parts

$H=H_{1}+H_{2}+H_{3}+H_{4}$ (A.10)

in such the way that Hamilton’s equations for each part can be integrated
explicitly. We used following splitting:

$H_{1}\equiv-\frac{p_{a}^{2}}{12V_{L}{M_{Pl}}^{2}},$ (A.11)
$H_{2}\equiv\sum_{\vec{x}}a^{4}\bigg{(}\frac{\pi^{2}_{\phi,\vec{x}}}{2a^{6}e^{2b(\chi_{\vec{x}})}}\bigg{)},$
(A.12)
$H_{3}\equiv\sum_{\vec{x}}a^{4}\bigg{(}\frac{\pi^{2}_{\chi,\vec{x}}}{2a^{6}}\bigg{)}$
(A.13)

and

$H_{4}\equiv\sum_{\vec{x}}a^{4}\bigg{(}e^{2b(\chi_{\vec{x}})}\frac{G_{\vec{x}}(\phi)}{2\delta^{2}a^{2}}+\frac{G_{\vec{x}}(\chi)}{2\delta^{2}a^{2}}+V(\phi_{\vec{x}},\chi_{\vec{x}})\bigg{)}.$
(A.14)

The corresponding flows for time step $\delta\tau$ are as follows:

$\Phi_{1}(h):\bigg{(}a,p_{a},\phi_{\vec{x}},\pi_{\phi,\vec{x}},\chi_{\vec{x}},\pi_{\chi,\vec{x}}\bigg{)}\rightarrow\bigg{(}a+\frac{\partial
H_{1}}{\partial
p_{a}}\delta\tau,p_{a},\phi_{\vec{x}},\pi_{\phi,\vec{x}},\chi_{\vec{x}},\pi_{\chi,\vec{x}}\bigg{)},$
(A.15)
$\Phi_{2}(h):\bigg{(}a,p_{a},\phi_{\vec{x}},\pi_{\phi,\vec{x}},\chi_{\vec{x}},\pi_{\chi,\vec{x}}\bigg{)}\rightarrow\bigg{(}a,p_{a}-\frac{\partial
H_{2}}{\partial a}\delta\tau,\phi_{\vec{x}}+\frac{\partial
H_{2}}{\partial\pi_{\phi,\vec{x}}}\delta\tau,\pi_{\phi,\vec{x}},\chi_{\vec{x}},\pi_{\chi,\vec{x}}-\frac{\partial
H_{2}}{\partial\chi_{\vec{x}}}\delta\tau\bigg{)},$ (A.16)
$\Phi_{3}(h):\bigg{(}a,p_{a},\phi_{\vec{x}},\pi_{\phi,\vec{x}},\chi_{\vec{x}},\pi_{\chi,\vec{x}}\bigg{)}\rightarrow\bigg{(}a,p_{a}-\frac{\partial
H_{3}}{\partial
a}\delta\tau,\phi_{\vec{x}},\pi_{\phi,\vec{x}},\chi_{\vec{x}}+\frac{\partial
H_{3}}{\partial\pi_{\chi,\vec{x}}}\delta\tau,\pi_{\chi,\vec{x}}\bigg{)}$
(A.17)

and

$\Phi_{4}(h):\bigg{(}a,p_{a},\phi_{\vec{x}},\pi_{\phi,\vec{x}},\chi_{\vec{x}},\pi_{\chi,\vec{x}}\bigg{)}\rightarrow\bigg{(}a,p_{a}-\frac{\partial
H_{4}}{\partial a}\delta\tau,\phi_{\vec{x}},\pi_{\phi,\vec{x}}-\frac{\partial
H_{4}}{\partial\phi_{\vec{x}}}\delta\tau,\chi_{\vec{x}},\pi_{\chi,\vec{x}}-\frac{\partial
H_{4}}{\partial\chi_{\vec{x}}}\delta\tau\bigg{)}.$ (A.18)

Using Strang splliting one can easliy construct second order symplectic method
which we use in our numerical simulations:

$\Phi(h)=\Phi_{1}(h/2)\circ\Phi_{2}(h/2)\circ\Phi_{3}(h/2)\circ\Phi_{4}(h)\circ\Phi_{3}(h/2)\circ\Phi_{2}(h/2)\circ\Phi_{1}(h/2).$
(A.19)

It can be written as the following set of explicit upgrades:

$a_{n+1/2}=a_{n}+\frac{\delta\tau}{2}{V_{L}}^{-1}\frac{\partial\tilde{H}_{1}}{\partial\wp_{a}}(\wp_{a,n})$
(A.20)
$\tilde{\wp}_{a,n+1/2}=\wp_{a,n}-\frac{\delta\tau}{2}{V_{L}}^{-1}\frac{\partial
H_{2}}{\partial a}(a_{n+1/2},\pi_{\phi,\vec{x},n},\chi_{\vec{x},n})$ (A.21)
$\tilde{\pi}_{\chi,\vec{x},n+1/2}=\pi_{\chi,\vec{x},n}-\frac{\delta\tau}{2}\frac{\partial
H_{2}}{\partial\chi_{\vec{x}}}(a_{n+1/2},\pi_{\phi,\vec{x},n},\chi_{\vec{x},n})$
(A.22)
$\phi_{\vec{x},n+1/2}=\phi_{\vec{x},n}+\frac{\delta\tau}{2}\frac{\partial
H_{2}}{\partial\pi_{\phi,\vec{x}}}(a_{n+1/2},\pi_{\phi,\vec{x},n},\chi_{\vec{x},n})$
(A.23)
$\tilde{\tilde{\wp}}_{a,n+1/2}=\tilde{\wp}_{a,n+1/2}-\frac{\delta\tau}{2}{V_{L}}^{-1}\frac{\partial
H_{3}}{\partial a}(a_{n+1/2},\tilde{\pi}_{\chi,\vec{x},n+1/2})$ (A.24)
$\chi_{\vec{x},n+1/2}=\chi_{\vec{x},n}+\frac{\delta\tau}{2}\frac{\partial
H_{3}}{\partial\pi_{\chi,\vec{x}}}(a_{n+1/2},\tilde{\pi}_{\phi,\vec{x},n+1/2})$
(A.25)
$\tilde{\tilde{\wp}}_{a,n+1}=\tilde{\tilde{\wp}}_{a,n+1/2}-h{V_{L}}^{-1}\frac{\partial
H_{4}}{\partial a}(a_{n+1/2},\phi_{\vec{x},n+1/2},\chi_{\vec{x},n+1/2})$
(A.26) $\pi_{\phi,\vec{x},n+1}=\pi_{\phi,\vec{x},n}-\delta\tau\frac{\partial
H_{4}}{\partial\phi{\vec{x}}}(a_{n+1/2},\phi_{\vec{x},n+1/2},\chi_{\vec{x},n+1/2})$
(A.27)
$\tilde{\pi}_{\chi,\vec{x},n+1}=\tilde{\pi}_{\chi,\vec{x},n+1/2}-\delta\tau\frac{\partial
H_{4}}{\partial\chi}(a_{n+1/2},\phi_{\vec{x},n+1/2},\chi_{\vec{x},n+1/2})$
(A.28)
$\chi_{\vec{x},n+1}=\chi_{\vec{x},n+1/2}+\frac{\delta\tau}{2}\frac{\partial
H_{3}}{\partial\pi_{\chi,\vec{x}}}(a_{n+1/2},\tilde{\pi}_{\phi,\vec{x},n+1})$
(A.29)
$\tilde{\wp}_{a,n+1}=\tilde{\tilde{\wp}}_{a,n+1}-\frac{\delta\tau}{2}{V_{L}}^{-1}\frac{\partial
H_{3}}{\partial a}(a_{n+1/2},\tilde{\pi}_{\chi,\vec{x},n+1})$ (A.30)
$\phi_{\vec{x},n+1}=\phi_{\vec{x},n+1/2}+\frac{h}{2}\frac{\partial
H_{2}}{\partial\pi_{\phi,\vec{x}}}(a_{n+1/2},\pi_{\phi,\vec{x},n+1},\chi_{\vec{x},n+1})$
(A.31)
$\pi_{\chi,\vec{x},n+1}=\tilde{\pi}_{\chi,\vec{x},n+1}-\frac{\delta\tau}{2}\frac{\partial
H_{2}}{\partial\chi_{\vec{x}}}(a_{n+1/2},\pi_{\phi,\vec{x},n+1},\chi_{\vec{x},n+1})$
(A.32)
$\wp_{a,n+1}=\tilde{\wp}_{a,n+1}-\frac{\delta\tau}{2}{V_{L}}^{-1}\frac{\partial
H_{2}}{\partial a}(a_{n+1/2},\pi_{\phi,\vec{x},n+1},\chi_{\vec{x},n+1})$
(A.33)
$a_{n+1}=a_{n+1/2}+\frac{\delta\tau}{2}{V_{L}}^{-1}\frac{\partial\tilde{H}_{1}}{\partial\wp_{a}}(\wp_{a,n+1})$
(A.34)

where we have defined:

$\tilde{H}_{1}(\wp_{a}):=H_{1}(V_{L}\wp_{a})=-\frac{V_{L}\wp_{a}^{2}}{12{M_{Pl}}^{2}}.$
(A.35)

## References

  * [1] A. A. Starobinsky, Phys. Lett. B 91 (1980), 99-102
  * [2] K. Sato, Mon. Not. Roy. Astron. Soc. 195 (1981), 467-479
  * [3] A. H. Guth, Phys. Rev. D 23 (1981), 347-356
  * [4] A. D. Linde, Phys. Lett. B 108 (1982), 389-393
  * [5] A. Albrecht and P. J. Steinhardt, Phys. Rev. Lett. 48 (1982), 1220-1223
  * [6] A. D. Linde, Phys. Lett. B 129 (1983), 177-181
  * [7] V. F Mukhanov Physical foundations of cosmology, Cambridge University Press 2005
  * [8] A. A. Starobinsky, JETP Lett. 30 (1979), 682-685
  * [9] V. F. Mukhanov and G. V. Chibisov, JETP Lett. 33 (1981), 532-535
  * [10] S. W. Hawking, Phys. Lett. B 115 (1982), 295
  * [11] A. A. Starobinsky, Phys. Lett. B 117 (1982), 175-178
  * [12] A. H. Guth and S. Y. Pi, Phys. Rev. Lett. 49 (1982), 1110-1113
  * [13] J. M. Bardeen, P. J. Steinhardt and M. S. Turner, Phys. Rev. D 28 (1983), 679
  * [14] A. A. Starobinsky, Phys. Lett. B 91 (1980) 99
  * [15] A. A. Starobinsky in: M. A. Markov and P. C. West (Eds.), “Quantum Gravity. Proceedings, 2nd Seminar, Moscow, USSR, October 13-15, 1981,” New York, Usa: Plenum (1984) 548pp
  * [16] A. De Felice and S. Tsujikawa, Living Rev. Rel. 13 (2010) 3 [arXiv:1002.4928 [gr-qc]].
  * [17] N. Aghanim et al. [Planck], Astron. Astrophys. 641 (2020), A6 [erratum: Astron. Astrophys. 652 (2021), C4] [arXiv:1807.06209 [astro-ph.CO]].
  * [18] Y. Akrami et al. [Planck], Astron. Astrophys. 641 (2020), A10 [arXiv:1807.06211 [astro-ph.CO]].
  * [19] J. Martin, C. Ringeval and V. Vennin, Phys. Dark Univ. 5-6 (2014), 75-235 [arXiv:1303.3787 [astro-ph.CO]].
  * [20] J. Martin, C. Ringeval, R. Trotta and V. Vennin, JCAP 03 (2014), 039 [arXiv:1312.3529 [astro-ph.CO]].
  * [21] S. Renaux-Petel, Comptes Rendus Physique 16 (2015), 969-985 [arXiv:1508.06740 [astro-ph.CO]].
  * [22] J. H. Traschen and R. H. Brandenberger, Phys. Rev. D 42 (1990) 2491.
  * [23] L. Kofman, A. D. Linde and A. A. Starobinsky, Phys. Rev. Lett. 73 (1994) 3195 [arXiv:hep-th/9405187].
  * [24] Y. Shtanov, J. H. Traschen and R. H. Brandenberger, Phys. Rev. D 51 (1995) 5438 [arXiv:hep-ph/9407247].
  * [25] L. Kofman, A. D. Linde and A. A. Starobinsky, Phys. Rev. D 56 (1997) 3258 [arXiv:hep-ph/9704452].
  * [26] B. A. Bassett, S. Tsujikawa and D. Wands, Rev. Mod. Phys. 78 (2006) 537 [arXiv:astro-ph/0507632].
  * [27] M. A. Amin, M. P. Hertzberg, D. I. Kaiser and J. Karouby, Int. J. Mod. Phys. D 24 (2014) 1530003 [arXiv:1410.3808 [hep-ph]].
  * [28] D. Battefeld, T. Battefeld, C. Byrnes and D. Langlois, JCAP 08 (2011), 025 [arXiv:1106.1891 [astro-ph.CO]].
  * [29] A. D. Linde, Phys. Lett. B 259 (1991), 38-47
  * [30] A. D. Linde, Phys. Rev. D 49 (1994), 748-754 [arXiv:astro-ph/9307002 [astro-ph]].
  * [31] E. J. Copeland, A. R. Liddle, D. H. Lyth, E. D. Stewart and D. Wands, Phys. Rev. D 49 (1994), 6410-6433
  * [32] D. H. Lyth and A. Riotto, Phys. Rept. 314 (1999), 1-146 [arXiv:hep-ph/9807278 [hep-ph]].
  * [33] S. Clesse, Phys. Rev. D 83 (2011), 063518 [arXiv:1006.4522 [gr-qc]].
  * [34] A. Avgoustidis, S. Cremonini, A. C. Davis, R. H. Ribeiro, K. Turzynski and S. Watson, JCAP 02 (2012), 038 [arXiv:1110.4081 [astro-ph.CO]].
  * [35] J. Martin and V. Vennin, Phys. Rev. D 85 (2012), 043525 [arXiv:1110.2070 [astro-ph.CO]].
  * [36] S. Renaux-Petel and K. Turzyński, Phys. Rev. Lett. 117 (2016) no.14, 141301 [arXiv:1510.01281 [astro-ph.CO]].
  * [37] S. Renaux-Petel, K. Turzyński and V. Vennin, JCAP 1711 (2017) no.11, 006 [arXiv:1706.01835 [astro-ph.CO]].
  * [38] S. Garcia-Saenz, S. Renaux-Petel and J. Ronayne, JCAP 1807 (2018) no.07, 057 [arXiv:1804.11279 [astro-ph.CO]].
  * [39] M. Cicoli, V. Guidetti, F. G. Pedro and G. P. Vacca, JCAP 1812 (2018) no.12, 037 [arXiv:1807.03818 [hep-th]].
  * [40] O. Grocholski, M. Kalinowski, M. Kolanowski, S. Renaux-Petel, K. Turzyński and V. Vennin, JCAP 05 (2019), 008 [arXiv:1901.10468 [astro-ph.CO]].
  * [41] M. Sasaki and E. D. Stewart, Prog. Theor. Phys. 95 (1996), 71-78 [arXiv:astro-ph/9507001 [astro-ph]].
  * [42] V. F. Mukhanov and P. J. Steinhardt, Phys. Lett. B 422 (1998), 52-60 [arXiv:astro-ph/9710038 [astro-ph]].
  * [43] C. Gordon, D. Wands, B. A. Bassett and R. Maartens, Phys. Rev. D 63 (2000), 023506 [arXiv:astro-ph/0009131 [astro-ph]].
  * [44] S. Groot Nibbelink and B. J. W. van Tent, Class. Quant. Grav. 19 (2002), 613-640 [arXiv:hep-ph/0107272 [hep-ph]].
  * [45] C. M. Peterson and M. Tegmark, Phys. Rev. D 87 (2013) no.10, 103507 [arXiv:1111.0927 [astro-ph.CO]].
  * [46] S. Renaux-Petel and K. Turzynski, JCAP 06 (2015), 010 [arXiv:1405.6195 [astro-ph.CO]].
  * [47] M. Wieczorek, “Niestabilność i nieliniowe oddziaływania perturbacji pól skalarnych w modelach inflacji kosmologicznej,” PhD dissertation, University of Warsaw 2019.
  * [48] S. Cremonini, Z. Lalak and K. Turzynski, JCAP 03 (2011), 016 [arXiv:1010.3021 [hep-th]].
  * [49] A. R. Brown, Phys. Rev. Lett. 121 (2018) no.25, 251601 doi:10.1103/PhysRevLett.121.251601 [arXiv:1705.03023 [hep-th]].
  * [50] S. Mizuno and S. Mukohyama, Phys. Rev. D 96 (2017) no.10, 103533 [arXiv:1707.05125 [hep-th]].
  * [51] T. Krajewski, K. Turzyński and M. Wieczorek, Eur. Phys. J. C 79 (2019) no.8, 654 doi:10.1140/epjc/s10052-019-7155-z [arXiv:1801.01786 [astro-ph.CO]].
  * [52] T. Krajewski and K. Turzyński, [arXiv:2204.12909 [astro-ph.CO]].