# Dynamics of Non-minimally Coupled Scalar Fields in the Jordan Frame

Daniel G. Figueroa<EMAIL_ADDRESS>Instituto de Física Corpuscular
(IFIC), Universitat de València-CSIC, E-46980, Valencia, Spain.  Adrien Florio
<EMAIL_ADDRESS>Center for Nuclear Theory, Department of Physics
and Astronomy, Stony Brook University Stony Brook, New York 11794, USA.  Toby
Opferkuch<EMAIL_ADDRESS>Berkeley Center for Theoretical Physics,
University of California, Berkeley, CA 94720  Theoretical Physics Group,
Lawrence Berkeley National Laboratory, Berkeley, CA 94720  Theoretical Physics
Department, CERN, Esplanade des Particules, 1211 Geneva 23, Switzerland  Ben
A. Stefanek<EMAIL_ADDRESS>Physik-Institut, Universität Zürich, CH-8057
Zürich, Switzerland

###### Abstract

The presence of scalar fields with non-minimal gravitational interactions of
the form $\xi|\phi|^{2}R$ may have important implications for the physics of
the early universe. While many studies solve the dynamics of non-minimally
coupled scalars in the Einstein frame, where gravity is simply described by
the Einstein-Hilbert action, we instead propose a procedure to solve the
dynamics directly in the original Jordan frame where the non-minimal couplings
are maintained explicitly. Our algorithm can be applied to scenarios that
include minimally coupled fields and an arbitrary number of non-minimally
coupled scalars, with the expansion of the universe sourced by all fields
present. This includes situations when the dynamics become fully
inhomogeneous, fully non-linear (due to e.g. backreaction or mode rescattering
effects), and/or when the expansion of the universe is dominated by non-
minimally coupled species. As an example, we study geometric preheating with a
non-minimally coupled scalar spectator field when the inflaton oscillates
following the end of inflation. In the future, our technique may be used to
shed light on aspects of the equivalence of the Jordan and Einstein frames at
the quantum level.

##  I Introduction

The dynamics of fields with non-minimal gravitational interactions may have
important implications for the physics of the early universe. In the case of
scalar field $\phi$ (either singlet or charged), one can add to the action a
marginal operator of the form $\xi|\phi|^{2}R$, where $R$ is the Ricci scalar
and $\xi$ is a real coupling constant controlling the strength of the
interaction. The presence of such term is actually required by the
renormalization properties of a scalar field in curved spacetime Parker2009 ;
Birrell1984 , where it is a running parameter that cannot be set to zero at
all energy scales.111Exceptionally, the running vanishes for the conformal
value $\xi=1/6$. In the case of the Standard Model (SM) Higgs field $\Phi$,
the operator $\xi|\Phi|^{2}R$ is actually the only missing operator of
dimension-4 that respects all symmetries of the SM and gravity. The coupling
$\xi$ can therefore be considered as the last unknown parameter of the SM.
However, due to the weakness of the gravitational interaction, current
particle physics experiments provide only extremely weak constraints on this
coupling, $\xi\lesssim 10^{15}$ Atkins:2012yn . It is therefore likely that
only early universe phenomena involving much higher energies than those
accessible to particle colliders can allow us to probe the non-minimal
gravitational interaction of the SM Higgs field, see e.g. Espinosa:2007qp ;
Enqvist:2014bua ; Espinosa:2015qea ; Ema:2016dny ; Enqvist:2016mqj ;
Figueroa:2017slm ; Markkanen:2018pdo ; Horn:2020wif ; Mantziris:2020rzh ;
Lebedev:2021xey .

Other fundamental (yet speculative) scalar fields may also have non-minimal
interactions with gravity. For instance, an early phase of accelerated
expansion in the Universe, known as inflation, is often assumed to be driven
by a scalar field called the inflaton, with an appropriate potential and
initial conditions (for reviews on inflation see e.g. Lyth:1998xn ;
Riotto:2002yw ; Bassett:2005xm ; Linde:2007fr ; Baumann:2009ds ). Indeed, a
scalar field with a non-minimal coupling to gravity can actually serve as a
good inflaton candidate, as any non-minimally coupled scalar theory can be
mapped via a conformal transformation to a minimally coupled theory with an
effective potential that can sustain inflation. Another popular realization of
inflation lies in modified gravity $f(R)$ theories, where $f$ is an arbitrary
function of $R$ (for a review on $f(R)$ theories, see e.g. DeFelice:2010aj ).
If $f^{\prime}(R)\neq 0$ and $f^{\prime\prime}(R)\neq 0$, there always exists
a mapping between the $f(R)$ theory and a scalar-tensor theory with a
propagating scalar degree of freedom non-minimally coupled to gravity with a
scalar potential purely of gravitational origin. As previously mentioned, this
setup can then be mapped onto a minimally coupled theory with an effective
potential suitable for inflation. A paradigmatic example of this is
Starobinsky inflation Starobinsky:1980te , defined by $f(R)=R+\alpha R^{2}$
with $\alpha>0$. After a conformal transformation to obtain a minimally
coupled theory, there is a scalar field – the scalaron – with a potential that
plateaus at large field amplitudes, naturally leading to inflation. Scenarios
where the inflaton has a non-minimal gravitational coupling $\xi|\phi|^{2}R$
lead to inflationary predictions in excellent agreement with current
observational constraints Planck:2018jri ; BICEP:2021xfz . This is independent
of whether the inflaton is of gravitational origin (as in Starobinsky
inflation) or elementary origin as in Higgs-Inflation, where the inflaton is
identified with the SM Higgs Bezrukov:2007ep . It is very interesting that
data from cosmological observations clearly favours plateau-like potentials
that naturally emerge in these scenarios Akrami:2018odb .

Non-minimally coupled inflaton scenarios can also lead to very interesting
phenomenology during the period after inflation. If the inflaton oscillates
around the minimum of its potential following inflation, particle species
coupled with sufficient strength are typically created in energetic bursts.
This non-perturbative process of particle production is known as preheating,
and it often leads to an exponential transfer of energy into particle sectors
(for reviews see Allahverdi:2010xz ; Amin:2014eta ; Lozanov:2019jxc ;
Allahverdi:2020bys ). This can occur in non-minimally coupled inflaton $\chi$
scenarios, where preheating into other degrees of freedom, e.g. scalar fields
$\\{\phi\\}$, can be realized very efficiently when couplings of the form
$g^{2}\chi^{2}\phi^{2}$ or $\xi\phi^{2}R$ are considered. Preheating scenarios
considering a simple monomial inflaton potential and a preheat scalar field
non-minimally coupled to gravity $\xi\phi^{2}R$ were first considered in
Bassett:1997az  and later on in Tsujikawa:1999jh , where an inflaton-preheat
field coupling $g^{2}\chi^{2}\phi^{2}$ was also included. The excitation of
the non-minimally coupled preheat field due to the oscillatory behavior of the
curvature $R$ (dictated by the oscillations of the inflaton field) was coined
as Geometric preheating in Ref. Bassett:1997az , and we keep that terminology
here. Preheating following inflation due to higher order curvature terms
$f(R)=R+\alpha_{n}R^{n}$ has been studied considering geometric preheating
effects in Tsujikawa:1999iv ; Fu:2019qqe . Preheating after Higgs-Inflation
was originally studied in Bezrukov:2008ut ; Garcia-Bellido:2008ycs ;
Figueroa:2009jw , and later on considered in more detail in Ema:2016dny ;
Sfakianakis:2018lzf . Preheating in multi-field inflationary scenarios with
$N$ scalar fields $\\{\phi_{j}\\}$ and couplings $\xi_{j}\phi_{j}^{2}R$ has
also been extensively studied DeCross:2015uza ; DeCross:2016fdz ;
DeCross:2016cbs ; Nguyen:2019kbm . Finally, preheating with a non-minimal
gravitational coupling $f(\phi)R$ with $f$ a general function of $\phi$, has
been also considered, see e.g. Watanabe:2006ku . All of the above
instabilities due to the presence of non-minimal gravitational couplings can
be generically regarded as “gravitational reheating” mechanisms222This
terminology does not apply to scenarios like Higgs-inflation, where only the
inflaton is non-minimally coupled to gravity, while the preheat fields are
directly coupled (non-gravitationally) to the inflaton., as they lead to very
efficient preheating, often exhibiting a violent transfer of energy among
fields. A different type of gravitational reheating mechanism was originally
put forward in Ford:1986sy  (see also Spokoiny:1993kt ). Namely, a massless
scalar field $\phi$ non-minimally coupled to gravity is excited towards the
end of inflation.333In Ford:1986sy  a non-conformal coupling was considered,
i.e. $\xi\neq 1/6$, but assuming a quasi-conformal window $0<|6\xi-1|\ll 1$.
Ref. Damour:1995pd  showed later on that the $(6\xi-1)^{2}$ suppression of the
energy density there is lifted to $\mathcal{O}(1)$ away from the quasi-
conformal window. The inflaton potential is chosen such that there is a
sustained kination dominated era with stiff equation of state $1/3<w\leq 1$
after inflation. As a consequence, the energy stored in $\phi$ (initially
suppressed compared to the inflaton energy) eventually becomes the dominant
energy component of the Universe. This idea is perhaps best exemplified in so-
called Quintessential inflation Peebles:1998qn ; Peloso:1999dm ; Huey:2001ae ;
Majumdar:2001mm ; Dimopoulos:2001ix ; Wetterich:2013jsa ; Wetterich:2014gaa ;
Hossain:2014xha ; Rubio:2017gty . The original gravitational reheating
mechanism, however, was shown in Ref. Figueroa:2018twl  to be inconsistent
with BBN/CMB constraints Caprini:2018mtu ; Clarke:2020bil  due to an excess
amount of gravitational wave production444In these scenarios the gravitational
wave background from inflation develops a large blue-tilt at high frequencies,
see e.g. Giovannini:2009kg ; Boyle:2007zx ; Figueroa:2019paj ;
Opferkuch:2019zbd ; Bernal:2020ywq .. More relevantly, it has been also shown
that a massless spectator field with a non-minimal coupling to gravity does
actually not scale as a radiation degree of freedom during kination domination
(as originally assumed in Ford:1986sy ; Spokoiny:1993kt ; Peebles:1998qn ),
but rather experiences a tachyonic instability due to the change in sign of
$R\propto(1-3w)$ for a stiff equation of state $w>1/3$. If the field is also
self-interacting, its energy grows due to the tachyonic instability until the
self-interaction eventually compensates the tachyonic mass. This was first
considered in Ref. Figueroa:2016dsc , with the SM Higgs as a spectator field
with a non-minimal coupling to gravity. There the universe is reheated into
relativistic SM particles after the Higgs experiences tachyonic growth during
a period of kination domination and then decays. The same mechanism was later
studied in more detail and extended to generic scalar fields non-minimally
coupled to gravity in Opferkuch:2019zbd , see also Dimopoulos:2018wfg ;
Bettoni:2021zhq . In Opferkuch:2019zbd , this mechanism was coined as Ricci
reheating and we stick to that nomenclature here.

All of the above scenarios exemplify the relevance of understanding the
dynamics of scalar fields non-minimally coupled to gravity in the early
universe. The form of the action where the non-minimal coupling to gravity is
maintained explicitly is known as the Jordan frame. In this frame, the
resulting equations of motion are difficult to solve in full generality due to
the non-linear feedback among them. Consequently, most studies rely on a
conformal transformation of the metric that brings the gravitational action to
the canonical Einstein-Hilbert form. This defines the so-called Einstein
frame, where the non-minimal coupling is absent and instead the kinetic terms
and scalar potentials of the matter fields are multiplied by a conformal
factor depending on the non-minimal coupling. After a field redefinition, one
can study the dynamics of the scalar field with Einstein gravity and a
canonical kinetic term, with all metric and field redefinition’s absorbed into
the scalar potential. Most of the studies cited above have worked out the
dynamics of non-minimally coupled scalar fields in the Einstein frame, or in
the linear regime in the Jordan frame, where analytic calculations can be
employed. The two frames are equivalent at the classical level, as long as the
map between them always exists. However, explicit examples exist where the
conformal map does not exist for all field values, such as in the
transformation from a non-minimally coupled theory to a minimally coupled
theory. In this case, the conformal map is given by
$\Omega^{2}=1-\xi(\phi/m_{p})^{2}$ which appears to be non-invertible for
$\phi^{2}=m_{p}^{2}/\xi$ and $\xi>0$ (here $m_{p}\simeq 2.4\cdot 10^{18}$ GeV
is the reduced Planck Mass). Furthermore, it is not known to what extent the
two frames are equivalent at the quantum level, as the conformal factor to
change from the Jordan to the Einstein frame is a local function of the non-
minimally coupled field $\Omega^{2}(\phi(x))$, with $\phi(x)$ often treated as
a quantum field. Some works evaluate the map using the vacuum expectation
value $\langle\phi^{2}\rangle$, but then it is not clear that the Einstein
frame description fully captures the physics of the theory originally written
in the Jordan frame.

These ambiguities can be simply resolved by solving the system directly in the
original Jordan frame. In this paper, we introduce a technique for precisely
this purpose. In particular, we are able to solve the dynamics of an arbitrary
scalar field $\phi$ with a non-minimal coupling to gravity $\xi\phi^{2}R$ in
an expanding background sourced by all fields present, even when the dynamics
become fully inhomogeneous and/or fully non-linear due to backreaction of the
excited species, including when the expansion of the universe is dominated by
the non-minimally coupled species. We can self-consistently evolve the
expansion of the universe while fully capturing field inhomogeneities and non-
linearities in the system, both of which typically develop very rapidly when
there are exponential instabilities like those typically arising in the
presence of a non-minimally coupled scalar field. As a working example, we
study geometric preheating effects involving a real scalar spectator field
non-minimally coupled to gravity, excited via an oscillatory effective mass
from $R$ that is sourced by oscillations of an inflaton with monomial
potential around its minimum.

##  II Continuum Dynamics in the Jordan frame

In this section we derive the equations of motion in the Jordan frame for a
theory with a non-minimally coupled scalar field. We consider a flat
Friedmann-Lemaître-Robertson-Walker (FLRW) background described by

$ds^{2}=-a(\eta)^{2\alpha}d\eta^{2}+a(\eta)^{2}\delta_{ij}dx^{i}dx^{j}\,,$ (1)

with $\eta$ an “$\alpha$-time” variable related to cosmic time by
$dt=a(\eta)^{\alpha}d\eta$. Here $\alpha$ is a (real number) parameter to be
conveniently chosen to suit each particular problem. Given the metric in Eq.
1, the Ricci scalar $R$ can be computed as (see Appendix A)

$R=\frac{6}{a^{2\alpha}}\left[\frac{a^{\prime\prime}}{a}+(1-\alpha)\left(\frac{a^{\prime}}{a}\right)^{2}\right]\,,$
(2)

where primes indicate derivatives with respect to $\eta$. We emphasize that
the Ricci scalar is a spatially homogeneous function, only depending on time,
as expected from consistency with Eq. 1. Let us consider a generic matter
sector ${\\{\varphi_{\rm m}\\}}$ minimally coupled to gravity, together with a
scalar field $\phi$ non-minimally coupled to gravity. Without loss of
generality, the action of this system reads

$\displaystyle\mathcal{S}$ $\displaystyle=\int
d^{4}x\sqrt{-g}\left(\frac{1}{2}m_{p}^{2}R-\frac{1}{2}\xi
R\phi^{2}-\frac{1}{2}g^{\mu\nu}\partial_{\mu}\phi\partial_{\nu}\phi-V(\phi,{\\{\varphi_{\rm
m}\\}})+\mathcal{L}_{\rm m}\right)\,,$ (3)

where $\frac{1}{2}m_{p}^{2}R$ is the standard Hibert-Einstein term,
$\frac{1}{2}\xi R\phi^{2}$ represents a non-minimal gravitational interaction
of $\phi$, and $V(\phi,{\\{\varphi_{\rm m}\\}})$ encompasses both the self
interactions of $\phi$ as well as its non-gravitational interactions with the
minimally-coupled matter sector. The term $\mathcal{L}_{\rm m}$ characterizes
the dynamics of the minimally-coupled fields, including their interactions and
self-interactions (which we do not specify explicitly here since they are
irrelevant for our discussion). Varying the action with respect to $\phi$, we
obtain the following equation of motion for $\phi$

$\displaystyle\Box\phi-\xi R\phi-\frac{\partial V}{\partial\phi}=0\,,$ (4)

where $\Box=g^{\mu\nu}\nabla_{\mu}\nabla_{\nu}$ and $\nabla_{\mu}$ is the
covariant derivative. We see that the non-minimal coupling introduces a term
proportional to $R$ in the equation of motion that acts a time dependent
effective mass for $\phi$. Using the $\alpha$-time metric given in Eq. 1, the
above equation becomes

$\displaystyle\phi^{\prime\prime}+(3-\alpha)\frac{a^{\prime}}{a}\phi^{\prime}-a^{-2(1-\alpha)}\nabla^{2}\phi+a^{2\alpha}\left(\xi
R\phi+\frac{\partial V}{\partial\phi}\right)=0\,.$ (5)

Equivalently, we can think of this gravitational interaction as part of an
effective potential

$\displaystyle V_{\rm eff}(\phi,{\\{\varphi_{\rm m}\\}},R)\equiv
V(\phi,{\\{\varphi_{\rm m}\\}})+\frac{1}{2}\xi R\phi^{2}\,,$ (6)

that includes all together the non-minimal coupling to gravity, the non-
gravitational interactions with the minimally-coupled matter sector, as well
as the self-interactions of $\phi$. For convenience, we can then think of a
Lagrangian for $\phi$ given by

$\displaystyle\mathcal{L}_{\phi}\equiv-\frac{1}{2}g^{\mu\nu}\partial_{\mu}\phi\partial_{\nu}\phi-
V_{\rm eff}(\phi,{\\{\varphi_{\rm m}\\}},R)\,.$ (7)

The Einstein equations are obtained by varying Eq. (3) with respect to
$g^{\mu\nu}$

$\displaystyle
G_{\mu\nu}=R_{\mu\nu}-\frac{1}{2}g_{\mu\nu}R=\frac{1}{m_{p}^{2}}T_{\mu\nu}\,,$
(8)

with

$\displaystyle T_{\mu\nu}$
$\displaystyle=-\frac{2}{\sqrt{-g}}\frac{\delta(\sqrt{-g}\mathcal{L}_{\rm
m})}{\delta
g^{\mu\nu}}-\frac{2}{\sqrt{-g}}\frac{\delta(\sqrt{-g}\mathcal{L}_{\phi})}{\delta
g^{\mu\nu}}\equiv T_{\mu\nu}^{\rm m}+T_{\mu\nu}^{\phi}\,,$ (9)

where $T_{\mu\nu}^{\rm m}$ and $T_{\mu\nu}^{\phi}$ have been defined as the
energy-momentum tensors of the minimally-coupled matter fields and the non-
minimally coupled scalar field $\phi$, respectively. In particular, one finds
(see Appendix B) the energy-momentum tensor of the non-minimally coupled field
to be

$\displaystyle T^{\phi}_{\mu\nu}$
$\displaystyle=\partial_{\mu}\phi\partial_{\nu}\phi-
g_{\mu\nu}\left(\frac{1}{2}g^{\rho\sigma}\partial_{\rho}\phi\partial_{\sigma}\phi+V\right)+\xi(G_{\mu\nu}+g_{\mu\nu}\Box-\nabla_{\mu}\nabla_{\nu})\phi^{2}\,.$
(10)

The trace of $T^{\phi}_{\mu\nu}$, defined as
$T_{\phi}=g^{\mu\nu}T^{\phi}_{\mu\nu}$, takes a simple form and will prove
very useful for simplifying the equations determining the evolution of the
scale factor. In $d+1$ spacetime dimensions, we find

$\displaystyle
T_{\phi}=(1-d)\frac{1}{2}\partial^{\mu}\phi\partial_{\mu}\phi-(d+1)V+\xi
G\phi^{2}+d\,\xi\Box\phi^{2}\,,$ (11)

where $G=g^{\mu\nu}G_{\mu\nu}=(1-d)R/2$ is the trace of the Einstein tensor
with respect to the background metric. Taking $d=3$ and using Eq. 4, we obtain

$\displaystyle
T_{\phi}=\left(6\xi-1\right)\left(\partial^{\mu}\phi\partial_{\mu}\phi+\xi
R\phi^{2}\right)+6\xi\phi V_{,\phi}-4V\,.$ (12)

where $V_{,\phi}=\partial V/\partial\phi$. Notably, if $\xi=1/6$, then for
$V=0$ or $V\propto\phi^{4}$ we find that $T_{\mu\nu}^{\phi}$ is traceless,
i.e. $T_{\phi}=0$, as a consequence of the conformal invariance of
$\mathcal{S}_{\phi}=\int d^{4}x\sqrt{-g}\mathcal{L}_{\phi}$ in these cases.
Given the FLRW metric in Eq. 1, the consistency of the Einstein equations
requires that $T_{\mu\nu}$ takes the form of the energy-momentum tensor of a
perfect fluid $T^{\mu}_{\,\,\,\,\nu}={\rm
diag}\left\\{-\rho(\eta),p(\eta),p(\eta),p(\eta)\right\\}$. We note that while
fields can develop large spatial inhomogeneities, the homogeneous and
isotropic pressure and energy density $p(\eta)$ and $\rho(\eta)$ should be
understood as the result of a volume average over the inhomogeneous local
field expressions. When the averaging volume is sufficiently large compared to
the excitation scales of the fields, this procedure leads to a well-defined
notion of a homogeneous and isotropic pressure and energy density within the
given volume. In this case, taking spatial averages over the off-diagonal
elements of $T_{\mu\nu}$ leads to vanishing results, consistent with
homogeneity and isotropy within the considered volume. Under these conditions,
the Einstein equations reduce to the Friedmann equations in $\alpha$-time

$\displaystyle\mathcal{H}^{2}$
$\displaystyle\equiv\left(\frac{a^{\prime}}{a}\right)^{2}=\frac{a^{2\alpha}}{3m_{p}^{2}}\rho(\eta)\,,$
(13) $\displaystyle\frac{a^{\prime\prime}}{a}$
$\displaystyle=-\frac{a^{2\alpha}}{6m_{p}^{2}}\Big{[}(1-2\alpha)\rho(\eta)+3p(\eta)\Big{]}\,,$
(14)

where we defined $\mathcal{H}=a^{\prime}/a$, which is related to the cosmic
time Hubble rate $H$ as $H=\mathcal{H}/a^{\alpha}$. We define the energy
density and pressure as

$\displaystyle\rho(\eta)$ $\displaystyle=\rho_{\phi}(\eta)+\rho_{\rm
m}(\eta)\equiv a^{-2\alpha}\langle T_{00}^{\phi}\rangle+a^{-2\alpha}\langle
T_{00}^{{\rm m}}\rangle\,,$ (15) $\displaystyle p(\eta)$
$\displaystyle=p_{\phi}(\eta)+p_{\rm
m}(\eta)\equiv\frac{1}{3a^{2}}\delta^{ij}\langle
T_{ij}^{\phi}\rangle+\frac{1}{3a^{2}}\delta^{ij}\langle T_{ij}^{{\rm
m}}\rangle\,,$ (16)

with $\langle\dots\rangle$ denoting volume averages. With these definitions,
the explicit expressions for the energy density and pressure of the non-
minimally coupled field $\rho_{\phi}$ and $p_{\phi}$ are found to be555The
volume averages of the total divergence terms
$\nabla^{2}\phi^{2}=\nabla\cdot\nabla\phi^{2}$ can be converted into surface
integrals that vanish in the case of an infinite volume with well-behaved
fields or in the case of a finite volume with periodic boundary conditions.

$\displaystyle\rho_{\phi}(\eta)$
$\displaystyle=\frac{1}{2a^{2\alpha}}\langle\phi^{\prime
2}\rangle+\frac{1}{2a^{2}}\langle(\nabla\phi)^{2}\rangle+\langle
V(\phi)\rangle+\frac{3\xi}{a^{2\alpha}}\mathcal{H}^{2}\langle\phi^{2}\rangle+\frac{6\xi}{a^{2\alpha}}\mathcal{H}\langle\phi\phi^{\prime}\rangle-\frac{\xi}{a^{2}}\langle\nabla^{2}\phi^{2}\rangle\,,$
(17) $\displaystyle p_{\phi}(\eta)$
$\displaystyle=\frac{(1-4\xi)}{2a^{2\alpha}}\langle\phi^{\prime
2}\rangle-\frac{(1-12\xi)}{6a^{2}}\langle(\nabla\phi)^{2}\rangle-\langle
V(\phi)\rangle+~{}\frac{2\xi}{a^{2\alpha}}\mathcal{H}\langle\phi\phi^{\prime}\rangle-\frac{\xi}{3a^{2}}\langle\nabla^{2}\phi^{2}\rangle$
$\displaystyle+2\xi\langle\phi
V_{,\phi}\rangle+\frac{\xi}{a^{2\alpha}}\left[\mathcal{H}^{2}+12\left(\xi-{1\over
6}\right)\left({a^{\prime\prime}\over
a}+(1-\alpha)\mathcal{H}^{2}\right)\right]\langle\phi^{2}\rangle\,.$ (18)

In principle, one can solve for the scale factor $a(\eta)$ from either Eq. 13
or Eq. 14. However, it is difficult in practice to solve these equations due
to their non-linear dependence on the derivatives of the scale factor. An
alternative approach is to relate the evolution of the scale factor to the
trace of the energy-momentum tensor, which only includes terms involving $R$
and the fields. An expedient way to do this is by computing the trace of Eq.
8, which gives

$\displaystyle R$
$\displaystyle=-\frac{1}{m_{p}^{2}}g^{\mu\nu}\left(T^{\phi}_{\mu\nu}+T^{\rm
m}_{\mu\nu}\right)=-\frac{1}{m_{p}^{2}}\left(T_{\phi}+T_{\rm m}\right)\,.$
(19)

Inserting the expression for $T_{\phi}$ given in Eq. 12, taking the volume
average of both sides, and solving for $R$, we find an expression only in
terms of the fields

$\displaystyle
R=\frac{F(\phi)}{m_{p}^{2}}\Big{[}\left(1-6\xi\right)\langle\partial^{\mu}\phi\partial_{\mu}\phi\rangle+4\langle
V\rangle-6\xi\langle\phi V_{,\phi}\rangle-\langle T_{\rm
m}\rangle\Big{]}\,,\hskip
35.56593ptF(\phi)\equiv\frac{1}{1+\left(6\xi-1\right)\xi\langle\phi^{2}\rangle/m_{p}^{2}}\,.$
(20)

This expression can be directly related to the evolution of the scale factor
using Eq. 2. This leads to the differential equation

$\displaystyle\frac{a^{\prime\prime}}{a}+(1-\alpha)\left(\frac{a^{\prime}}{a}\right)^{2}=\frac{a^{2\alpha}F(\phi)}{6m_{p}^{2}}\Big{[}\left(1-6\xi\right)\langle\partial^{\mu}\phi\partial_{\mu}\phi\rangle+4\langle
V\rangle-6\xi\langle\phi V_{,\phi}\rangle-\langle T_{\rm m}\rangle\Big{]}\,,$
(21)

that together with the equation of motion for $\phi$ in Eq. 5, will allow us
to spell out a simple and concise numerical scheme to evolve this system. To
start, it is convenient write the equations in terms of natural variables, by
rescaling fields and coordinates as

$\displaystyle\tilde{\phi}=\frac{1}{f_{*}}\phi\,,~{}~{}\differential\tilde{\eta}=\omega_{*}\differential\eta\,,~{}~{}\differential\tilde{x}_{i}=\omega_{*}\differential
x_{i}\,.$ (22)

with $f_{*}$ some typical field amplitude and $\omega_{*}$ a characteristic
(inverse) time scale of the problem to be studied. The choice of $f_{*}$ and
$\omega_{*}$ depends entirely on the scenario at hand (we will provide an
explicit example in Section IV. We also need introduce an appropriate
rescaling of the matter sector (see Ref. Figueroa:2020rrl  for examples). If
the matter sector simply comprises of a set of scalar fields $\\{\varphi_{\rm
m}\\}$, these are normalized as in Eq. 22. We note that rescaling the
coordinates by $\omega_{*}$ naturally induce the following rescaling in $R$

$\tilde{R}=\omega_{*}^{-2}R\ .$ (23)

It is also natural to introduce rescaled energy densities and pressure

$\displaystyle\widetilde{V}=\frac{1}{f_{*}^{2}\omega_{*}^{2}}V\,,\hskip
42.67912pt{\tilde{\rho}}=\frac{1}{f_{*}^{2}\omega_{*}^{2}}\rho\,,\hskip
42.67912pt{\tilde{p}}=\frac{1}{f_{*}^{2}\omega_{*}^{2}}p\,.$ (24)

Next, we reduce the order of the equation of motion for $\tilde{\phi}$ by
introducing a conjugate momentum variable as

$\displaystyle\tilde{\pi}_{\phi}=a^{3-\alpha}\tilde{\phi}^{\prime}\,.$ (25)

The matter sector is treated in a similar way, with the rescaling of the
conjugate momenta variables depending on the spin of the species, see
Figueroa:2020rrl . If the matter sector is comprised of scalar fields, we
simply introduce a set of conjugate momenta $\\{\tilde{\pi}_{\varphi_{\rm
m}}\\}$, analogously to Eq. (25). In the new variables, the evolution of the
non-minimally coupled scalar field is governed by a system of coupled first-
order differential equations, in terms of a kernel functional
$\tilde{\mathcal{K}}_{\phi}$, as follows

$\displaystyle\left\\{\begin{array}[]{l}\tilde{\phi}^{\prime}=a^{\alpha-3}\tilde{\pi}_{\phi}\,,\vspace*{0.2cm}\\\
\tilde{\pi}^{\prime}_{\phi}=\tilde{\mathcal{K}}_{\phi}[a,\tilde{\phi},\\{\tilde{\varphi}_{\rm
m}\\},\tilde{R}]\,,~{}~{}~{}~{}~{}~{}~{}{\rm
with}~{}~{}~{}~{}~{}~{}~{}\mathcal{K}_{\phi}[a,\tilde{\phi},\\{\tilde{\varphi}_{\rm
m}\\},\tilde{R}]\equiv
a^{1+\alpha}\,\tilde{\nabla}^{2}\tilde{\phi}-a^{3+\alpha}\left(\xi\tilde{R}\tilde{\phi}+\frac{\partial\widetilde{V}}{\partial\tilde{\phi}}\right)\,.\end{array}\right.$
(28)

Similarly, to evolve the scale factor we use Eq. 21 as derived from the trace
of the energy-momentum tensor. Defining the conjugate momentum of $a(\eta)$ as

$\displaystyle\pi_{a}=a^{1-\alpha}a^{\prime}\,,$ (29)

we arrive to a system of coupled first-order differential equations depending
on another kernel functional,

$\displaystyle\left\\{\begin{array}[]{l}a^{\prime}=a^{\alpha-1}\tilde{\pi}_{a}\,,\vspace*{0.2cm}\\\
\tilde{\pi}_{a}^{\prime}=\tilde{\mathcal{K}}_{a}[a,\tilde{R}]\,,~{}~{}~{}~{}~{}~{}~{}{\rm
with}~{}~{}~{}~{}~{}~{}~{}\tilde{\mathcal{K}}_{a}[a,\tilde{R}]\equiv\frac{a^{{2+\alpha}}}{6}\tilde{R}\
.\end{array}\right.$ (32)

To close the system, an expression for $\tilde{R}$ is needed in both kernels
$\mathcal{K}_{\phi},~{}\mathcal{K}_{a}$. Using Eq. 20, we can write

$\displaystyle\tilde{R}$
$\displaystyle=\frac{f_{*}^{2}}{m_{p}^{2}}\left[\frac{2\left(1-6\xi\right)\big{(}\tilde{E}_{G}^{\tilde{\phi}}-\tilde{E}_{K}^{\tilde{\phi}}\big{)}+4\langle\tilde{V}\rangle-6\xi\langle\tilde{\phi}\,\tilde{V}_{,\tilde{\phi}}\rangle+({\tilde{\rho}}_{\rm
m}-3{\tilde{p}}_{\rm
m})}{1+\left(6\xi-1\right)\xi\langle\tilde{\phi}^{2}\rangle
f_{*}^{2}/m_{p}^{2}}\right]\,,$ (33)

where we have used $\langle T_{\rm m}\rangle=3p_{\rm m}-\rho_{\rm m}$ and
introduced the volume-averaged kinetic $\tilde{E}_{K}^{\tilde{\phi}}$ and
gradient $\tilde{E}_{G}^{\tilde{\phi}}$ energy densities

$\displaystyle\tilde{E}_{K}^{\tilde{\phi}}=\frac{1}{2a^{6}}\langle\tilde{\pi}_{\tilde{\phi}}^{2}\rangle\,,\hskip
42.67912pt\tilde{E}_{G}^{\tilde{\phi}}=\frac{1}{2a^{2}}\sum_{i}\langle\tilde{\partial}_{i}\tilde{\phi}\partial_{i}\tilde{\phi}\rangle\
.$ (34)

In summary, Eqs. 28 and 32, together with the expression for $\tilde{R}$ in
Eq. 33 (plus the equations of motion of the unspecified matter sector),
represent a set of equations that completely characterizes the dynamics of a
system with a scalar field non-minimally coupled to gravity in the Jordan
frame. Generalization to multiple non-minimally coupled scalars is obtained
straight forwardly by summing over the terms with non-minimal coupling in Eq.
33.

##  III Lattice formulation

In order to evolve our system of equations Eqs. 28, 32 and 33 in a way that
fully captures the spatial dependence of the fields, we need to choose a time
evolution scheme and to introduce a spatial discretization prescription. We
use a lattice with $N$ sites per dimension with periodic boundary conditions.
We will consider the lattice sites to represent comoving coordinates. If the
(comoving) length of the grid is $L$, the resulting (comoving) lattice spacing
between sites is $\delta x=L/N$. We work with finite differences and use the
following notation for the forward and backward derivatives

$\displaystyle\nabla_{i}^{\pm}f({\bf n})=\frac{\pm f({\bf n})\mp f({\bf
n}\pm\hat{i})}{\delta x}\,,$ (35)

where $f$ is an arbitrary scalar function defined on the lattice sites ${\bf
n}=(n_{1},n_{2},n_{3})$, and $\hat{i}$ represents a displacement vector of one
unit in the $i$-th direction. We discretize the gradient terms using forward
differences and the Laplacian using a symmetric discretization

$\displaystyle\sum_{i}\langle\partial_{i}\phi\partial_{i}\phi\rangle$
$\displaystyle~{}~{}\longrightarrow~{}~{}\sum_{i}\langle\nabla_{i}^{+}\phi\nabla_{i}^{+}\phi\rangle\,,$
(36) $\displaystyle\vec{\nabla}^{2}\phi$
$\displaystyle~{}~{}\longrightarrow~{}~{}\sum_{i}\nabla^{-}_{i}\nabla_{i}^{+}\phi\,.$
(37)

We are now in a position to define the evolution equations by introducing the
discrete kernels

$\displaystyle\tilde{\mathcal{K}}_{\phi}\left[a,\tilde{\phi},\\{\tilde{\varphi}_{\rm
m}\\},\tilde{R}\right]$
$\displaystyle=a^{1+\alpha}\tilde{\nabla}^{-}_{i}\tilde{\nabla}_{i}^{+}\tilde{\phi}-a^{3+\alpha}\left(\xi\tilde{R}\tilde{\phi}+\frac{\partial\widetilde{V}}{\partial\tilde{\phi}}\right)\,,$
(38) $\displaystyle\tilde{\mathcal{K}}_{a}\left[a,\tilde{R}\right]$
$\displaystyle=\frac{a^{{2+\alpha}}}{6}\tilde{R}\,,$ (39)

which we have already written in terms of natural field and spacetime
variables, c.f. Eq. 22. We have also introduced dimensionless discrete
derivatives $\tilde{\nabla}$ given by Eq. (35) in terms of the dimensionless
lattice spacing $\delta\tilde{x}=\tilde{L}/N=\omega_{*}\delta x$, with
$\tilde{L}=\omega_{*}L$.

At this point, it is important to realize that
$\tilde{R}=\tilde{R}[\tilde{\phi},\tilde{\pi}_{\phi},\\{\tilde{\varphi}_{\rm
m}\\},\\{\tilde{\pi}_{\varphi_{\rm m}}\\}]$ depends on all fields and
conjugate momentum variables, and hence the kernel for the non-minimally
coupled field $\tilde{\phi}$ depends on its own conjugate momentum. Because of
this, preferred symplectic algorithms such as staggered Leapfrog, velocity- or
position-Verlet, cannot be used (see Ref. Figueroa:2020rrl  for a discussion
on this). We can instead use Runge-Kutta (RK) methods, in particular explicit
RK algorithms. We have adapted the well known mid-point method to our set of
equations, corresponding to a second order RK method. To account for
situations where a high time-accuracy may be required, we have also
implemented a particularly interesting family of explicit low-storage RK
methods of higher order following Refs. WILLIAMSON198048 ; Bazavov:2021pik .
These present multiple advantages: they are easy to implement, the memory cost
does not increase when increasing the accuracy order, and in some cases an
adaptive time-step scheme is allowed. The interested reader can find more
information and an explicit description of all these RK algorithms applied to
our system of equations in Appendix C.

One last important point is to have a discrete version of the Hubble
constraint given in Eq. 13. Verifying that this constraint is preserved by our
numerical evolution scheme provides an important check of the method. In terms
of rescaled variables, it reads

$\displaystyle\left(\frac{a^{\prime}}{a}\right)^{2}=\frac{a^{2\alpha}f_{*}^{2}}{3m_{p}^{2}}\bigg{[}{\tilde{\rho}}_{{\rm
m}}+\tilde{E}_{K}^{\tilde{\phi}}+\tilde{E}_{G}^{\tilde{\phi}}+\langle\widetilde{V}\rangle+\frac{3\xi}{a^{2\alpha}}\left(\frac{a^{\prime}}{a}\right)^{2}\langle\tilde{\phi}^{2}\rangle+\frac{6\xi}{a^{\alpha+3}}\left(\frac{a^{\prime}}{a}\right)\langle\tilde{\phi}\,\tilde{\pi}_{\tilde{\phi}}\rangle\bigg{]}\,,$
(40)

where we have dropped the $\langle\nabla^{2}\phi^{2}\rangle$ term of Eq. 17
because it is a total derivative whose volume average vanishes due to the
periodic boundary conditions of the lattice. We now have all the tools to
evolve our system of equations on the lattice. In the next section, we present
an explicit example in the context of geometric preheating. Lastly, note that
all numerical algorithms presented above have been implemented in the package
$\mathcal{C}\mathrm{osmo}\mathcal{L}\mathrm{attice}$ Figueroa:2020rrl ;
Figueroa:2021yhd , which can perform user-friendly and versatile field theory
simulations. These new algorithms will be made publicly available in a future
update of $\mathcal{C}\mathrm{osmo}\mathcal{L}\mathrm{attice}$.

##  IV Example: Geometric Preheating

We now study an example of geometric preheating directly in the Jordan frame,
using the formalism developed in the previous sections. By geometric
preheating, we refer to the excitation of a light spectator field $\phi$ non-
minimally coupled to gravity. This occurs due to the oscillatory behavior of
the spacetime curvature $R$ that follows after inflation, when a homogeneous
inflaton field oscillates around the minimum of its potential Bassett:1997az
illustrated in Fig. 1. The fact that $R$ becomes oscillatory can be seen from
the traced Einstein equations, assuming that the homogeneous inflaton field
$\chi$ initially dominates the energy density of the universe, such that
$T=T_{\phi}+T_{\rm inf}\approx T_{\rm
inf}=-\partial^{\mu}\chi\partial_{\mu}\chi-4V_{\rm inf}(\chi)$. This leads to

Figure 1: Time evolution of the effective potential of the non-minimally
coupled field $\phi$ as a function of the Ricci scalar. Oscillations between
an unbounded and bounded potential occur once the Ricci scalar begins
oscillating at around one e-fold.

$\displaystyle
R=-\frac{T}{m_{p}^{2}}\approx\frac{1}{m_{p}^{2}}\left[\partial^{\mu}\chi\partial_{\mu}\chi+4V_{\rm
inf}(\chi)\right]=\frac{1}{m_{p}^{2}}\left[4V_{\rm
inf}(\chi)-\frac{1}{a^{2\alpha}}\chi^{\prime 2}\right]\,.$ (41)

One illustrative example is the case where $V_{\rm inf}(\chi)={1\over
2}m^{2}\chi^{2}$, in which case $R$ can be approximated in cosmic time
($\alpha=0$) as

$\displaystyle R=\frac{1}{m_{p}^{2}}\left[4V_{\rm
inf}(\chi)-\dot{\chi}^{2}\right]\approx\frac{R_{0}}{4}\left(\frac{a_{0}}{a}\right)^{3}\left[1+3\cos(2mt)\right]\,.$
(42)

In this expression, it is manifest that $R$ oscillates between positive and
negative values due to the harmonic oscillations of $\chi$. In general, for an
inflaton potential with a minimum around the origin and an arbitrary power law
behavior $V_{\rm inf}(\chi)\propto|\chi|^{p}$ ($p>1$), (or even for a linear
combination of various power laws) the oscillations of the inflaton will not
be harmonic. This does not change the fact that $R$, and hence the effective
mass squared of the spectator field $m_{\phi,{\rm eff}}^{2}=\xi R$, will still
alternate periodically between positive and negative values. As a consequence
of the periodic tachyonic stages ($m_{\phi,{\rm eff}}^{2}<0$), initial quantum
vacuum fluctuations of the spectator field $\phi$ can be exponentially
amplified if the strength of its non-minimal coupling $\xi$ is large enough.
The amplification may persist until the effective tachyonic mass of $\phi$ is
fully screened by its own self-interactions, or until the energy of $\phi$
grows to the same order as the energy available in the system. In either case,
a detailed lattice study is required due to the non-linearity of the system.
We present first in Section IV.1 a linear analysis of the initial instability
of the mode functions of $\phi$, then in Section IV.2 we present an analysis
of the evolution of the system once the dynamics become non-linear.

###  IV.1 Initial Conditions via Linear Analysis

Our procedure will consist of computing the power spectrum of the $\phi$
fluctuations induced during inflation and the subsequent transition period via
a linear analysis, which we then use as the initial condition for the lattice
evolution before the dynamics enter the non-linear regime. To proceed, we
consider a theory involving an inflaton field $\chi$ and a light spectator
field $\phi$ with a non-minimal coupling to gravity, similar to Refs.
Bassett:1997az ; Bertolami:2010ke ; Fu:2019qqe ; Tsujikawa:1999iv

$\displaystyle\mathcal{S}=\int\mathop{\mathrm{\mathstrut
d}}\nolimits^{4}\\!x\sqrt{-g}\left[\frac{m_{p}^{2}}{2}R+\mathcal{L}_{\phi}+\mathcal{L}_{\text{inf}}\right]\,,$
(43)

with

$\displaystyle\mathcal{L}_{\phi}$
$\displaystyle=-\frac{1}{2}g^{\mu\nu}\partial_{\mu}\phi\partial_{\nu}\phi-\frac{1}{2}\xi
R\phi^{2}-V(\phi)\,,$ (44) $\displaystyle\mathcal{L}_{\text{inf}}$
$\displaystyle=-\frac{1}{2}g^{\mu\nu}\partial_{\mu}\chi\partial_{\nu}\chi-
V_{\text{inf}}(\chi)\,.$ (45)

In this theory, the inflaton $\chi$ and the spectator field $\phi$ interact
only gravitationally through the non-minimal coupling $\xi R\phi^{2}$. During
slow-roll inflation, we have a quasi de-Sitter phase where $R\approx 12H^{2}$
and $H\approx$ constant (its time derivative is slow-roll suppressed). This
means that for $\xi>0$, the spectator field has a heavy effective mass
$m_{\phi,{\rm eff}}^{2}\approx 12\xi H^{2}$ during inflation. We assume that
this effective mass dominates over $V^{\prime\prime}(\phi)$, such that the
potential $V(\phi)$ can be neglected during inflation. This, combined with the
fact that the non-minimally coupled spectator field $\phi$ is energetically
subdominant during inflation, justifies the use of a linear analysis. It will
be convenient to work in conformal time ($\alpha=1$) where the metric is
conformally flat and quantization proceeds as in Minkowski space. In that
case, we can write the action for $\phi$ in terms of the canonically
normalized field $\varphi=a\phi$

$\mathcal{S}_{\varphi}=\frac{1}{2}\int d\tau
d^{3}x\left[(\varphi^{\prime})^{2}-(\nabla\varphi)^{2}-a^{2}\left(\xi-\frac{1}{6}\right)R\varphi^{2}\right]\,.$
(46)

We then canonically quantize $\varphi$ as

$\hat{\varphi}(x,\tau)=\int\frac{d^{3}k}{(2\pi)^{3}}\left[\varphi_{k}(\tau)\hat{a}_{\bf
k}e^{i{\bf k\cdot x}}+\varphi^{*}_{k}(\tau)\hat{a}_{\bf k}^{\dagger}e^{-i\bf
k\cdot x}\right]\,,$ (47)

where $[\hat{a}_{\bf k},\hat{a}_{\bf
k^{\prime}}^{\dagger}]=(2\pi)^{3}\delta({\bf k-k^{\prime}})$ and the modes are
normalized such that
$\varphi_{k}\varphi_{k}^{\prime*}-\varphi^{\prime}_{k}\varphi_{k}^{*}=i$. The
mode functions $\varphi_{k}(\tau)$ obey the equation of motion given by Eq. 46
in momentum space

$\varphi_{k}^{\prime\prime}(\tau)+\left[k^{2}+a^{2}\left(\xi-\frac{1}{6}\right)R\right]\varphi_{k}(\tau)=0\,.$
(48)

We assume that the evolution of each mode starts far inside the Hubble radius,
namely $-k\tau\gg 1$, where the curvature is negligible. In that case, the
modes should approach the _Bunch-Davies_ vacuum

$\displaystyle\varphi_{k}(k\tau\to-\infty)$
$\displaystyle=\frac{1}{\sqrt{2k}}e^{-ik\tau}\,.$ (49)

We are interested in the power spectrum $\Delta_{\varphi}(k,\tau)$, which in
terms of the two-point function is defined as

$\displaystyle\langle\varphi^{2}\rangle=\langle
0|\hat{\varphi}(\tau,x)\hat{\varphi}(\tau,x)|0\rangle=\int\frac{dk}{k}\Delta_{\varphi}(k,\tau)\,.$
(50)

The power spectrum of the original field $\phi$ is then related as

Figure 2: Resulting power spectrum of the non-minimally coupled spectator
field from the linear analysis performed in Section IV.1. Colored lines
illustrate the spectrum at the indicated number of e-folds post inflation.

$\displaystyle\Delta_{\phi}(k,\tau)=\frac{1}{a^{2}}\Delta_{\varphi}(k,\tau)=\frac{k^{3}}{2\pi^{2}a^{2}}|\varphi_{k}(\tau)|^{2}\equiv\frac{k^{3}}{2\pi^{2}}\mathcal{P}_{\phi}(k,\tau)\,.$
(51)

For our numerical results (which involve integrating over many e-folds of
inflaton), it proves easier to solve Eq. (48) in cosmic time ($\alpha=0$),
where it reads

$\displaystyle\ddot{\varphi}_{k}+H\dot{\varphi}_{k}+\left[\frac{k^{2}}{a^{2}}+\left(\xi-\frac{1}{6}\right)R\right]\varphi_{k}=0\,,$
(52)

with the Bunch-Davies initial condition now expressed as

$\displaystyle\varphi_{k}(k/(aH)\gg
1)\approx\frac{1}{\sqrt{2k}}e^{\frac{ik}{aH}}\,,$ (53)

where we have used $\tau\approx-1/(aH)$, given that $H$ changes very slowly.
As previously mentioned, the non-minimally coupled spectator field $\phi$ is
initially energetically subdominant by assumption, so we neglect its
contribution to the background evolution for the linear analysis. In this
case, all modes evolve independently and the background evolution is
completely determined by the homogeneous energy components of the inflaton

$H^{2}\equiv\left(\frac{\dot{a}}{a}\right)^{2}=\frac{1}{3m_{p}^{2}}\left(\frac{1}{2}\dot{\chi}^{2}+V_{\rm
inf}(\chi)\right)\,,$ (54)

with the evolution of the homogenous inflaton governed by the standard Klein-
Gordon equation of motion in cosmic time

$\ddot{\chi}+3H\dot{\chi}+\frac{\partial V_{\rm inf}}{\partial\chi}=0\,.$ (55)

In our numerical analysis we consider an observationally viable inflationary
model inspired by $\alpha$-attractors Kallosh:2013hoa , with inflaton
potential parametrised as Planck:2018jri

$\displaystyle V_{\text{inf}}(\chi)$
$\displaystyle=\Lambda^{4}\tanh^{p}\left(\frac{c\chi}{m_{p}}\right)\,,\qquad\text{with}\qquad
p=4,6\,,$ (56)

which flattens out for $|\chi|\gg m_{p}/c$ (where $c$ is a dimensionless
parameter), and takes a power-law form $V\propto\chi^{p}$ for $|\chi|\ll
m_{p}/c$. We take $c=0.1$ which reproduces the observed value of the scalar
perturbations at CMB scales for $V_{\inf}(\chi_{\rm
CMB})=($1.6\text{\times}{10}^{16}\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$)^{4}$
and saturates the upper bound on the scale of inflation Planck:2018jri
corresponding to $\Lambda=1.79\times 10^{16}$ GeV.

We solve Eq. 52 numerically by discretizing $k$ on a grid of 512 log-spaced
modes. We begin evolving each mode considering the Bunch-Davies initial
condition when $k/(aH)=\beta$, with $\beta\gg 1$ a penetration factor. Larger
values of $\beta$ better approximate the Bunch-Davies initial condition, but
also increase simulation time, so as a compromise we choose $\beta=10^{3}$. At
the end of inflation, we would like to have a superhorizon power spectrum of
simulated modes spanning at least three orders of magnitude in $k$-space.
Since the lowest $k$ mode starts a factor $\beta$ inside the horizon, we
require approximately $\Delta N\simeq\log(10^{3}\beta)\approx 14$ e-folds of
simulated inflation for all modes of interest to exit the horizon. We
therefore choose the initial conditions of the homogeneous inflaton field
$\chi_{i}$ such that we obtain 14 e-folds of inflation (this corresponds to
$\chi_{i}=9.23(10.87)m_{p}$ for $p=4(6)$, respectively). Following this
procedure, we numerically integrate Eqs. (52), (54) and (55) for 14 e-folds of
inflation and through the transition to the post-inflationary stage. We show
the resulting power spectrum in Fig. 2 for $\xi=50$.

For comparison, our numerical results for $\phi$ can be compared to the
predicted power spectrum during inflation in pure de-Sitter space, which was
computed analytically in Ref. Opferkuch:2019zbd  as

$H_{*}^{-2}\,\Delta_{\phi}(z)=\frac{z^{3}}{8\pi}|H_{i\mu}^{(1)}(z)|^{2}e^{-\pi\mu}\approx\frac{1}{4\pi^{2}}\begin{cases}\frac{z^{3}}{\mu}\,,&z\ll
1\hskip 11.38109pt({\rm superhorizon})\vspace*{1mm}\\\ z^{2}\,,&z\gg 1\hskip
11.38109pt({\rm subhorizon})\\\ \end{cases}\,.$ (57)

Here $z=k/(aH_{*})$, $\mu^{2}=12(\xi-3/16)$, and $H_{*}$ is taken to be the
Hubble rate at the end of inflation. The approximate equality holds for
$\xi>3/16$. According to this expression, we see that the superhorizon
fluctuations during inflation follow a $k^{3}$ power law. This is expected as
after Hubble radius exit the modes are damped because of the heavy effective
mass induced by the non-minimal coupling. On the other hand, the subhorizon
modes deep inside the Hubble radius remain in the Bunch-Davies vacuum $\propto
k^{2}$, indicating that they are not excited. The transition between power
laws occurs around $k/(aH_{*})\approx\mu$, where we have
$\mu\approx\sqrt{12\xi}$ for $\xi\gg 1$. We see that this analytic
approximation explains well the behavior of the power spectrum shown in Fig. 2
at the end of inflation.

After inflation ends, the inflaton begins oscillating around the minimum of
its potential, which also induces oscillations in $R$, as shown in Fig. 1.
This allows for tachyonic growth of the non-minimally coupled field during the
periods when $R<0$, leading to the emergence of a peak in its power spectrum
after inflation. That fact that the peak is stationary when the power spectrum
is plotted in terms of $k/(aH)$ can be obtained from inspecting Eq. 52, as the
tachyon is regularized (on a per mode basis) when $k/(aH)\approx\sqrt{|\xi
R/H^{2}|}\approx\mathcal{O}(1)\sqrt{\xi}\lesssim\sqrt{6\xi}$, where we used
that tachyonic growth happens for $-6\leq R/H^{2}<0$ and assumed $\xi\gg 1$.
Correspondingly, the peak in the power spectrum, which clearly emerges about
$\sim$1 e-fold after the end of inflation, is observed at
$k/(aH)\approx\mathcal{O}(1)\sqrt{\xi}$. It is at this moment when we
introduce the power spectrum from our linear analysis, in order to initialize
the non-minimally coupled field in the lattice simulation. The shape of the
peak of the power spectrum determines the range of comoving momenta we need to
consider on the lattice. In terms of comoving momenta, the peak shifts to
smaller values of $k$ while its amplitude grows during the linear regime.
Hence, the most important scales to capture in the lattice are those spanned
by the peak itself and its infrared tail to some extent, so there is room for
the peak to shift further to the infrared as we simulate the dynamics in the
lattice.

###  IV.2 Non-linear Lattice Analysis

Figure 3: Power spectra of the non-minimally coupled spectator field $\phi$ at
the end of inflation defined by $N=N_{*}$ and $N-N_{*}$ e-folds thereafter.
The solid/dashed lines are the results of the lattice/linear simulations.
Left: Power spectrum for the case of $\xi=50$ where the spectator field has no
potential, $V=0$. Initially the spectrum is well described by the linear
analysis until the growth becomes large enough that backreaction occurs at
approximately $N-N_{*}=1.75$. Right: Power spectrum for $\xi=100$ where the
spectator field has a small quartic coupling $\lambda=${10}^{-5}\text{\,}$$.
Deviations from the linear analysis begin almost immediately where we see the
non-linear effect of the quartic term transferring power to higher $k$-modes.

After a clear peak in the power spectrum emerges in the linear analysis, but
still before any interactions of $\phi$ or its backreaction onto the
background dynamics are relevant, we move to solving the system on the
lattice. In particular, we treat $\phi$ as a classical field whose initial
fluctuations are drawn randomly from a Gaussian distribution with power
spectrum given by $\Delta_{\phi}$ that we computed in the linear analysis up
to one e-fold after the end of inflation. We now fully include all
interactions as well as a self-consistent evolution of the expanding
background including contributions by all fields present. The relevant
equations of motions to solve are

$\displaystyle\phi^{\prime\prime}+(3-\alpha)\frac{a^{\prime}}{a}\phi^{\prime}-\frac{\nabla^{2}\phi}{a^{2(1-\alpha)}}+a^{2\alpha}\xi
R\phi$ $\displaystyle=-a^{2\alpha}\frac{\partial V}{\partial\phi}\,,$ (58)
$\displaystyle\chi^{\prime\prime}+(3-\alpha)\frac{a^{\prime}}{a}\chi^{\prime}-\frac{\nabla^{2}\chi}{a^{2(1-\alpha)}}$
$\displaystyle=-a^{2\alpha}\frac{\partial V_{\text{inf}}}{\partial\chi}\,,$
(59)
$\displaystyle\frac{a^{\prime\prime}}{a}+(1-\alpha)\left(\frac{a^{\prime}}{a}\right)^{2}$
$\displaystyle=\frac{a^{2\alpha}}{6}R\,,$ (60)

where Eq. 60 determines the background evolution as discussed in Section II
and $R$ is defined as

$\displaystyle
R=\frac{F(\phi)}{m_{p}^{2}}\left[\left(6\xi-1\right)\left(\frac{1}{a^{2\alpha}}\langle\phi^{\prime
2}\rangle-\frac{1}{a^{2}}\langle(\nabla\phi)^{2}\rangle\right)-6\xi\langle\phi
V_{,\phi}\rangle+4\langle V+V_{\rm
inf}\rangle-\frac{1}{a^{2\alpha}}\langle\chi^{\prime
2}\rangle+\frac{1}{a^{2}}\langle(\nabla\chi)^{2}\rangle\right]\,,$ (61)

with $F(\phi)$ given in Eq. (20). We adapt the continuum equations to the
lattice following Section III and Appendix C. All simulations are run on
lattices of size $N=240$ points per spatial dimension and evolved with the
RK2MP method described in detail in Appendix C. In the lattice, we use natural
variables as defined in Eq. 22, with $f_{*}=m_{p}$ and $\omega_{*}=H_{i}$,
where $H_{i}$ is the Hubble rate at the start of the linear analysis $\simeq
14$ e-folds before the end of inflation.

Figure 4: Left: Evolution of the Ricci scalar normalized to Hubble squared
$H(t)^{2}$. The dashed-black line is the evolution of the linear free-field
analysis which is used as an input to the lattice simulation. The colored
lines $\xi=\\{10,50,100\\}$, $\\{$yellow, blue, green$\\}$ are the lattice
simulation results which illustrate a strong deviation once the energy density
of the non-minimally coupled field is comparable to that of the inflaton.
Shaded in red is the region where the Ricci scalar is negative, driving
tachyonic growth of the spectator field. Right: Expectation value
$\langle\phi^{2}\rangle$, see Eq. 50, of the non-minimally coupled spectator
field for the same values of $\xi$ as before.

Figure 5: Evolution of the total equation of state as a function of e-folds
post inflation. Here we define $w\equiv p(\eta)/\rho(\eta)$ using the non-
minimally coupled fields contributions from Eq. 17 and Eq. 18 plus the
standard minimally coupled contributions from the inflaton. We show the
results for $p=4$ ($p=6$) on the left (right), respectively. Figure 6: Time
evolution of the energy density of both the inflaton (red) and non-minimally
coupled spectator field (black). Parameter values for the inflationary
potential are $c=0.1$,
$\Lambda\simeq$1.79\text{\times}{10}^{16}\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$$,
while the energy density is normalized as in Eq. 24. We show the absolute
value of the energy density but change the line style to dashed to illustrate
when the energy density turns negative. Top row: Three values of the non-
minimal coupling $\xi$ for the inflationary potential with $p=4$. In the top-
right panel we also show the effect of a Hubble scale mass (thin purple) and a
small quartic coupling of the spectator field (thin gray). Bottom row: Same
three values of $\xi$ but for $p=6$, where the energy density of the inflaton
redshifts faster than radiation.

Let us first consider the case where $V(\phi)=0$. The resulting evolution of
the Ricci scalar in the $p=4$ case for the inflaton potential is illustrated
by the dashed-black line in the left-hand panel of Fig. 4. Here we see that
shortly after the end of inflation the value of $R/H^{2}$ tends to negative
values before oscillating in the range $-6\leq R/H^{2}\leq 12$. While $R$ is
negative, the non-minimally coupled field has a tachyonic effective mass and
can experience exponential growth. This can easily be seen in the right-hand
panel of Fig. 4 where we see the expectation value $\langle\phi^{2}\rangle$
growing exponentially when $R$ takes negative values, as expected. For
comparison purposes, the dashed lines show the results of the linear analysis,
where the field grows unbounded when the tachyonic growth is strong enough to
overcome the expansion of the universe (as in the $\xi=50,100$ cases) since we
neglect the contribution of $\phi$ to the background evolution. This is not
the case in the lattice analysis, where for these large values of $\xi$ the
backreaction of $\phi$ on the background evolution asymptotically drives
$\xi\langle\phi^{2}\rangle/m_{p}^{2}$ to a constant value below unity, and $R$
ends oscillating around zero with a damped amplitude, as shown by the solid
lines in Fig. 4. More specifically, at the onset of backreaction the kinetic
term $(6\xi-1)\langle\phi^{\prime 2}\rangle$ drives initially $R$ to a large
positive value, as represented by the green (blue) spikes for $\xi=100$ ($50$)
in the left panel of Fig. 4. This results in a large positive effective mass
squared $\xi R$, that induces a restoring force for $\phi$, opposing its
growth. The field velocity $\phi^{\prime}$ is then suppressed causing $R$ to
start a rapid descent down to a small negative value, after which it begins
oscillating around zero. After the spike, $R$ oscillates with a small damped
amplitude, so the successive tachyonic mass stages cannot overtake anymore the
expansion of the universe, and $\xi\langle\phi^{2}\rangle/m_{p}^{2}$
approaches asymptotically a constant value. If we define the equation of state
in terms of the total pressure and energy densities $w=p(\eta)/\rho(\eta)$,
then $w$ can be written in terms of $R$ by combining Eqs. 2, 13 and 14

$w=\frac{p(\eta)}{\rho(\eta)}=\frac{1}{3}\left(1-\frac{R}{3H^{2}}\right)\,,$
(62)

where $H=\mathcal{H}/a^{\alpha}$. We see that the large positive spikes where
$R/H^{2}>12$ when the non-minimally coupled field backreacts correspond to
periods where the equation of state spikes below $w=-1$, as shown in the left-
hand figure of Fig. 5. Since $\rho$ is always constrained to be positive
definite by Eq. 13, it can be seen from Eq. 62 that $R/H^{2}>12$ corresponds
to a large negative pressure, namely $p<-\rho$, which violates the classical
energy conditions due to the non-minimal interaction of $\phi$ with the
gravitational field. In particular, the dominant contribution to the pressure
during these spikes comes from the $(1-4\xi)\langle\phi^{\prime 2}\rangle$
term in Eq. 18, which is always negative for $\xi>1/4$.

Finally, in Fig. 6 we show the evolution of the energy density of both the
inflaton $\chi$ (red) and the non-minimally coupled spectator field $\phi$
(black). In the top (bottom) row we consider the hypertangent inflaton
potential with $p=4$ ($p=6$), while in the three columns we again consider the
three benchmark values of $\xi=\\{10,50,100\\}$. In the $p=4$ case, the
inflaton energy density drops like radiation since the potential is quartic
around the minimum, while in the $p=6$ case the inflaton energy density decays
faster than radiation. Though the total energy density is always positive, it
is well known that the energy density of the non-minimally coupled field as
defined in Eq. 17 is not positive definite and we indicate when it becomes
negative using dashed lines. In the cases where $V(\phi)=0$, the energy
density of the non-minimally coupled field scales as radiation at late times,
as can be seen in Fig. 5. In the upper right-hand panel of Fig. 6 we also show
the behavior when the non-minimally coupled field has a non-zero potential
$V(\phi)$. We consider the cases of $V(\phi)=m^{2}\phi^{2}/2$ with mass
$m/H_{i}=0.02$ (thin purple line) as well as $V(\phi)=\lambda\phi^{4}/4$ with
a quartic $\lambda=${10}^{-5}\text{\,}$$ (thin gray line). The case
$V(\phi)=m^{2}\phi^{2}/2$ exhibits similar behavior to the $V(\phi)=0$ case
until around 2 e-folds where bare mass term begins to dominate over the
effective mass induced by $R$ and the energy density of the non-minimally
coupled field begins to scale like matter. This allows the energy density of
the non-minimally coupled field to dominate over that of the inflaton,
providing an efficient reheating mechanism. The effect of the quartic is much
more drastic because it acts to regulate the tachyonic growth occurring when
$R$ is negative, preventing the non-minimally coupled field from reaching
large field values. After a transition period where the energy density dilutes
faster than radiation due to the $\xi$ dependent terms in Eq. 17, the energy
becomes dominated by the field oscillating in its quartic potential which
leads to radiation scaling of the energy density. This case does not result in
the non-minimally coupled field fully reheating the universe unless the
inflaton energy drops faster than radiation. This is precisely what occurs in
the $p=6$ case, shown in the bottom row of Fig. 6. In this case, the energy
density of the non-minimally coupled field can quickly come to dominate over
that of the inflaton.

##  V Summary and Conclusions

The presence of at least one fundamental scalar field in the SM raises the
question of the role non-minimal couplings to gravity may play in the
evolution of the early Universe. Any scalar $\phi$ in curved spacetime, be it
the SM Higgs or otherwise, inevitably acquires a non-minimal coupling to
gravity of the form $\xi|\phi|^{2}R$ through renormalization group evolution.
Typically, the dynamics of non-minimally coupled scalars are not studied
directly in the original Jordan frame, but rather in the Einstein frame via a
conformal transformation of the metric that brings the action to the canonical
Einstein-Hilbert form. This approach allows for a more intuitive
interpretation of the dynamics, but is not entirely satisfactory because the
equivalence of the two frames at the quantum level is currently unclear.

In this work, we have developed an approach to solve the dynamics of non-
minimally coupled scalars in an expanding universe directly in the original
Jordan frame, where the non-minimal couplings are maintained explicitly. In
the Jordan frame, the equations of motion describing the background evolution
are typically non-linear in the derivatives of the scale factor, making them
difficult to solve in practice. We tackle this problem by considering the
trace of the energy-momentum tensor, a simpler object that can be related to
the background evolution. This admits a simple system of coupled first-order
differential equations for the background evolution that can be
straightforwardly numerically integrated. In Section III, we demonstrate how
this method can be implemented in the
$\mathcal{C}\mathrm{osmo}\mathcal{L}\mathrm{attice}$ Figueroa:2020rrl ;
Figueroa:2021yhd  package by specifying the discrete evolution kernels. There,
we see that the resulting kernel for the non-minimally coupled field evolution
depends on its own conjugate momentum, preventing the use of symplectic
algorithms typically employed. We have therefore implemented explicit “low-
storage” Runge-Kutta methods that allow for high-order methods while keeping
memory usage constant and permitting adaptive time-steps, see appendix C for
further details.

To demonstrate the viability of our method, we study geometric preheating as
an illustrative example in Section IV. The model involves a real spectator
scalar field $\phi$ that is excited through its non-minimal coupling to
gravity when the inflaton oscillates around the minimum of its potential
following the end of inflation. The oscillations of the inflaton source
oscillations in $R$, inducing a time-dependent effective mass for $\phi$.
While the effective squared mass is negative, tachyonic growth of the non-
minimally coupled field can occur if the value of $\xi$ is large enough to
overcome the friction due to the expansion of the universe, as shown in Fig.
6. In this case, the growth of the non-minimally coupled field is highly
efficient and the energy density of $\phi$ reaches an $\mathcal{O}(1)$ value
of the total energy density within an $\mathcal{O}(1)$ number of e-folds,
representing an extremely efficient preheating mechanism. We find that if
there is no explicit scale in the potential of the non-minimally coupled
field, then its energy density scales as radiation at late times. Therefore,
in the cases where $\xi$ is large enough, whether the inflaton or non-
minimally coupled field dominates the energy density at late times depends on
the choice of potential for both fields.

To conclude, we have introduced a robust method to solve the dynamics of non-
minimally coupled scalar fields directly in the original Jordan frame with the
expansion sourced by all fields present, even when the dynamics becomes fully
inhomogeneous and/or non-linear due to the backreaction of the excited
species. This will provide an important tool to study the equivalence of the
Einstein and Jordan frame at the quantum level, in particular when the initial
conditions of the fields are set by quantum vacuum fluctuations in both
frames. All numerical algorithms presented in this paper will be made publicly
available in a future update to the package
$\mathcal{C}\mathrm{osmo}\mathcal{L}\mathrm{attice}$, which will allow anyone
to perform user-friendly and versatile cosmological simulations involving non-
minimally coupled scalar fields.

## Acknowledgments

TO would like to thank Valerie Domcke, Nadav Outmezguine and Marko Simonovic
for insightful discussions. DGF (ORCID 0000-0002-4005-8915) is supported by a
Ramón y Cajal contract with Ref. RYC-2017-23493. This work is also supported
by project PROMETEO/2021/083 from Generalitat Valenciana, and by project
PID2020-113644GB-I00 from Ministerio de Ciencia e Innovación.

##  Appendix A Curvature in $\alpha$-time

Given a flat “$\alpha$-time” FLRW background described by the metric

$ds^{2}=-a(\eta)^{2\alpha}d\eta^{2}+a(\eta)^{2}\delta_{ij}dx^{i}dx^{j}\,,$
(63)

the Christoffel symbols can be computed from

$\Gamma_{\mu\nu}^{\sigma}=\frac{1}{2}g^{\sigma\rho}\left(\partial_{\nu}g_{\mu\rho}+\partial_{\mu}g_{\nu\rho}-\partial_{\rho}g_{\mu\nu}\right)\,.$
(64)

The non-vanishing Christoffel symbols are found to be

$\Gamma_{00}^{0}=\alpha\left(\frac{a^{\prime}}{a}\right)\,,\hskip
42.67912pt\Gamma_{ij}^{0}=\delta_{ij}\left(\frac{a^{\prime}}{a}\right)\,a^{2(1-\alpha)}\,,\hskip
42.67912pt\Gamma_{0j}^{i}=\Gamma_{j0}^{i}=\delta^{i}_{j}\left(\frac{a^{\prime}}{a}\right)\,,$
(65)

where primes indicate derivatives with respect to $\eta$. From this, the
components of the Ricci tensor

$R_{\mu\nu}=\partial_{\rho}\Gamma^{\rho}_{\mu\nu}-\partial_{\nu}\Gamma^{\rho}_{\mu\rho}+\Gamma_{\mu\nu}^{\rho}\Gamma^{\sigma}_{\rho\sigma}-\Gamma_{\mu\sigma}^{\rho}\Gamma_{\nu\rho}^{\sigma}\,,$
(66)

are found to be

$\displaystyle R_{00}$
$\displaystyle=\frac{3}{a^{2\alpha}}\left[\frac{a^{\prime\prime}}{a}-\alpha\left(\frac{a^{\prime}}{a}\right)^{2}\right]g_{00}\,,\hskip
56.9055ptR_{ij}=\frac{1}{a^{2\alpha}}\left[\frac{a^{\prime\prime}}{a}+(2-\alpha)\left(\frac{a^{\prime}}{a}\right)^{2}\right]g_{ij}\,.$
(67)

From this, the Ricci scalar $R=g^{\mu\nu}R_{\mu\nu}$ can be computed directly.
We find

$R=\frac{6}{a^{2\alpha}}\left[\frac{a^{\prime\prime}}{a}+(1-\alpha)\left(\frac{a^{\prime}}{a}\right)^{2}\right]\,.$
(68)

##  Appendix B Energy-momentum tensor of a non-minimally coupled scalar

The energy momentum tensor for the non-minimally coupled $\phi$ is defined in
Eq. (9) as

$T^{\phi}_{\mu\nu}=-\frac{2}{\sqrt{-g}}\frac{\delta(\sqrt{-g}\mathcal{L_{\phi}})}{\delta
g^{\mu\nu}}=-2\frac{\delta\mathcal{L_{\phi}}}{\delta
g^{\mu\nu}}+g_{\mu\nu}\mathcal{L_{\phi}}\,.$ (69)

The remaining variation reads

$\begin{split}\frac{\delta\mathcal{L_{\phi}}}{\delta
g^{\mu\nu}}&=-\frac{1}{2}\partial_{\mu}\phi\partial_{\nu}\phi-\frac{1}{2}\xi\frac{\delta(\phi^{2}R)}{\delta
g^{\mu\nu}}\,.\end{split}$ (70)

A useful identity one can show by explicit computation from the expression of
$R$ in terms of the metric is the following

$\frac{\delta(fR)}{\delta
g^{\mu\nu}}=fR_{\mu\nu}+(g_{\mu\nu}\nabla^{\sigma}\nabla_{\sigma}-\nabla_{\mu}\nabla_{\nu})f\,,$
(71)

for any scalar function $f$ and $\nabla_{\nu}$ is the covariant derivative
associated to $g_{\mu\nu}$. Applying this to Eq. (70), we obtain

$\displaystyle\frac{\delta\mathcal{L_{\phi}}}{\delta g^{\mu\nu}}$
$\displaystyle=-\frac{1}{2}\partial_{\mu}\phi\partial_{\nu}\phi-\frac{1}{2}\xi\phi^{2}R_{\mu\nu}-\frac{1}{2}\xi(g_{\mu\nu}\nabla^{\sigma}\nabla_{\sigma}-\nabla_{\mu}\nabla_{\nu})\phi^{2}\,,$

and putting it altogether we get

$\displaystyle T^{\phi}_{\mu\nu}$
$\displaystyle=\partial_{\mu}\phi\partial_{\nu}\phi-
g_{\mu\nu}\left(\frac{1}{2}g^{\rho\sigma}\partial_{\rho}\phi\partial_{\sigma}\phi+V(\phi)\right)+\xi\left(R_{\mu\nu}-\frac{1}{2}Rg_{\mu\nu}+g_{\mu\nu}\Box-\nabla_{\mu}\nabla_{\nu}\right)\phi^{2}\,.$
(72)

where we have defined
$\Box\equiv\nabla^{\sigma}\nabla_{\sigma}=g^{\rho\sigma}\nabla_{\rho}\nabla_{\sigma}$.
We would like to compute the trace of this energy-momentum tensor for $\phi$.
For a moment, let us write down the result working in $d+1$ dimensions. Then,
we find

$\displaystyle
T_{\phi}=(1-d)\frac{1}{2}\partial^{\mu}\phi\partial_{\mu}\phi-(d+1)V+\xi
G\phi^{2}+d\,\xi\Box\phi^{2}\,,$ (73)

where $G=g^{\mu\nu}G_{\mu\nu}=(1-d)R/2$ is the trace of the Einstein tensor
with respect to the metric. This expression can be further simplified using
$\Box\phi^{2}=2\phi\Box\phi+2\partial^{\mu}\phi\partial_{\mu}\phi$ and the
equation of motion Eq. 5 for $\phi$ which gives $\phi\Box\phi=\xi
R\phi^{2}+\phi V_{,\phi}$. This leads to the following expression

$\displaystyle g^{\mu\nu}T^{\phi}_{\mu\nu}$
$\displaystyle=2d\left[\xi-\frac{(d-1)}{4d}\right]\left(\partial^{\mu}\phi\partial_{\mu}\phi+\xi
R\phi^{2}\right)-V\left[(d+1)-2\xi d\,\phi\frac{V_{,\phi}}{V}\right]\,.$ (74)

Note that the coefficient of $\left(\partial^{\mu}\phi\partial_{\mu}\phi+\xi
R\phi^{2}\right)$ vanishes for $\xi=(d-1)/(4d)$, which is indeed the conformal
value of $\xi$ in $d+1$ dimensions. Now setting $d=3$, we get

$\displaystyle g^{\mu\nu}T^{\phi}_{\mu\nu}$
$\displaystyle=\left(6\xi-1\right)\left(\partial^{\mu}\phi\partial_{\mu}\phi+\xi
R\phi^{2}\right)-\left(4V-6\xi\phi V_{,\phi}\right)\,.$ (75)

which indeed leads to Eq. 12 for $T_{\phi}=g^{\mu\nu}T^{\phi}_{\mu\nu}$ quoted
in the main text. Let us comment on some specific cases of Eq. (75) which
reproduce known results. Consider the conformal value for $\xi$ in $d+1=4$
dimensions, $\xi=1/6$. Then, we find

$g^{\mu\nu}T^{\phi}_{\mu\nu}=-\left(4V-\phi V_{,\phi}\right)\,,$ (76)

which vanishes identically in the scaleless cases of $V=0$ or
$V\propto\phi^{4}$, as expected. For a quadratic potential
$V=m^{2}\phi^{2}/2$, we get

$g^{\mu\nu}T^{\phi}_{\mu\nu}=-m^{2}\phi^{2}\,,$ (77)

which reproduces the result of Ref. Parker:2009uva .

##  Appendix C Time evolution and low-storage RK methods

In this appendix, we present the implemented ‘low-storage’ Runge-Kutta (RK)
methods . We also write down an explicit algorithm to evolve our system of
equations using these methods. We begin by recalling the reader some facts
about RK methods, following Ref. Figueroa:2020rrl . Consider a vector
$\vec{x}(t)$ of $M$-variables $\vec{x}(t)=(x_{1}(t),\dots,x_{M}(t))^{T}$ and a
system of first order differential equations of the type

$\displaystyle\dot{\vec{x}}(t)=\mathcal{K}\left[\vec{x}(t)\right]\,.$ (78)

Then, a RK method of order $s$ is characterized by a one-step iteration of the
type

$\displaystyle\vec{x}_{n+1}=\vec{x}_{n}+\sum_{i=1}^{s}c_{i}k^{(i)}\,,$ (79)

with

$\displaystyle x^{(i)}$
$\displaystyle=\vec{x}_{n}+\sum_{j=1}^{s}b_{ij}k^{(j)}\,,$ (80) $\displaystyle
k^{(i)}$ $\displaystyle=\mathcal{K}\left[\vec{x}^{(i)}\right]\,,$ (81)
$\displaystyle\sum_{i=1}^{s}c_{i}$ $\displaystyle=1\,,\ \ \,c_{i}<1\,.$ (82)

This iteration effectively split the time interval $\delta t$ into $s$
subintervals $\delta t=\sum_{i=1}^{s}c_{i}\delta t$. Note also that, after
having introduced conjugate momenta, this is precisely the type of equations
we are dealing with. These methods are often represented in terms of Butcher
tableaux

$\displaystyle\begin{array}[]{r|cccc|l}&b_{11}&b_{12}&\cdots&b_{1s}&\\\
&b_{21}&b_{22}&\cdots&b_{2s}&\\\ &\ddots&\ddots&\cdots&\ddots&\\\
&b_{s1}&b_{s2}&\cdots&b_{ss}&\\\
\hline\cr&c_{1}&c_{2}&\cdots&c_{s}&\end{array}\,.$ (88)

Explicit RK methods have the property $b_{ij}=0$ for all $j\geq i$. Well known
methods of order two are the modified Euler-method (RK2ME) and the midpoint
method (RK2MP). For comparison we also show the Butcher tableau of the widely
used RK4 algorithm below:

$\displaystyle{\tt RK2ME}:~{}\begin{array}[]{r|cc|l}&0&0&\\\ &1&0&\\\
\hline\cr&1/2&1/2&\end{array}\qquad\qquad{\tt
RK2MP}:~{}\begin{array}[]{r|cc|l}&0&0&\\\ &1/2&0&\\\
\hline\cr&0&1&\end{array}\qquad\qquad{\tt
RK4}:~{}\begin{array}[]{r|cccc|l}&0&0&0&0&\\\ &1/2&0&0&0&\\\ &0&1/2&0&0&\\\
&0&0&1&0&\\\
\hline\cr&1/6&1/3&1/3&1/6&\end{array}.~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
(100)

In cases where the limiting factor is memory, such as when solving a system of
partial differential equations on large lattices, the memory cost of using
higher-order RK methods can become prohibitive. Indeed, generically, one needs
to solve (almost) all of the $k^{(i)}$ coefficients. For a method with
$s$-stages, the required additional memory is analogous to simulating $s$ new
fields per field and momentum.

Interestingly, there exists a subclass of RK methods which eludes this memory
requirement; they are referred to as ‘low-storage’ RK methods WILLIAMSON198048
(see also Ref. Bazavov:2021pik  for a recent application in lattice QCD). This
method hinges on rewriting of Eqs. 79, 80 and 81 as

$\displaystyle\vec{x}_{n+1}=\vec{y}^{(s)}\,,$ (101)

with

$\displaystyle\vec{y}^{(0)}$ $\displaystyle=\vec{x}_{n}\,,$ (102)
$\displaystyle\vec{y}^{(i)}$
$\displaystyle=\vec{y}^{(i-1)}+B_{i}\Delta\vec{y}^{(i)}\,,$ (103)
$\displaystyle\Delta\vec{y}^{(i)}$
$\displaystyle=A_{i}\Delta\vec{y}^{(i-1)}+\delta
t\,,\mathcal{K}\left[\vec{y}^{(i-1)}\right]\,,$ (104)

with the further requirement $A_{1}=0$. Note that all second-order, and some
third-order RK methods can be put in this form; we refer the interested reader
to Refs. WILLIAMSON198048 ; Bazavov:2021pik  and those therein for more
information. It is easy to see that the second order methods introduced above
can be recast in this form using the following coefficients:

$\displaystyle\texttt{RK2ME}:~{}\begin{array}[]{r|cc|l}&A_{i}&B_{i}&\\\
\hline\cr&0&1&\\\
&~{}-1{}&1/2&\end{array}\qquad\qquad\qquad\qquad\texttt{RK2MP}:~{}\begin{array}[]{r|cc|l}&A_{i}&B_{i}&\\\
\hline\cr&0&1/2&\\\ &-1/2&1&\end{array}$ (111)

We have also implemented the following third order method from Ref.
Carpenter1994Thirdorder2R  which is argued to have desirable stability
properties. The coefficients below are the rational form of the ones presented
in Section 3. of this reference, for $c_{3}=(2+10^{1/3})/6$:

$\displaystyle\texttt{RK3\\_4}:~{}\begin{array}[]{r|cc|l}&A_{i}&B_{i}&\\\
\hline\cr&0&0.06688758201974097&\\\ &-0.7825460361923583&2.876554598956719&\\\
&-2.042914325731225&0.5534657361343982&\\\
&-1.799337253940777&0.3912730180961791&\end{array}$ (117)

Finally, fourth order $2N$-storage schemes also exist, we refer the interested
reader to Refs. Carpenter1994Fourthorder2R ; Bazavov:2021pik  for examples.

Before writing down explicitly the $2N$-storage method applied to our problem,
we note that the scheme RK3_4 has the additional property that the third
iteration $\vec{y}^{(3)}$ is already at second order accuracy in $\delta t$.
At the extra memory cost of saving the previous solution $\vec{x}_{n}$ in case
the update fails, one can then easily turn it into an adaptive time-step RK
scheme. As reviewed in Refs. Carpenter1994Thirdorder2R ; Bazavov:2021pik ,
this can be achieved by estimating the distance $\Delta$ in some norm between
$\vec{y}^{(3)}$ and $\vec{y}^{(4)}$. If this distance is smaller than some
requested tolerance $\epsilon$, the update is accepted; if not it is rejected
and the step is repeated. In both situation, the time step is updated to

$\displaystyle\delta
t_{\text{new}}=0.95\cdot\left(\frac{\epsilon}{\Delta}\right)^{1/3}\delta
t_{\text{old}}\,.$ (118)

This updated always decrease $\delta t$ when the time step needs to be
repeated and almost always increases it when the error is below the requested
tolerance. The factor $0.95$ and the power $1/3$ are empirical determined
based on performance. A practical way to define $\Delta$ is to compute the
Euclidean distance between the solutions

$\displaystyle\Delta=|\vec{y}^{(3)}-\vec{y}^{(4)}|=B_{4}|\Delta\vec{y}^{(4)}|$
(119)

with $|\vec{y}|=\sqrt{\sum_{i=1}^{L}y_{l}^{2}}$ for the $L$-component vector
$\vec{y}=(y_{1},\dots,y_{L})^{T}$.

We are now in a position to present a concrete algorithm to evolve the
equations presented in Section III. For every field, momentum and scale
factor, we introduce associated auxiliary variables:
$\Delta\tilde{\phi},\Delta\tilde{\pi}_{\tilde{\phi}},\\{\Delta\tilde{\varphi}_{\rm
m}\\},\\{\Delta\tilde{\pi}_{\tilde{\varphi}_{\rm m}}\\},\Delta
a,\Delta\pi_{a}$. We can then implement a generic $s$-stage $2N$-storage RK
method as follows

$\displaystyle\left.\begin{array}[]{l}\tilde{\pi}_{\tilde{\phi}}^{(0)}\equiv\tilde{\pi}_{\tilde{\phi}}({\bf
n},n_{0})\vspace*{0.15cm}\\\ \tilde{\phi}^{(0)}\equiv\tilde{\phi}({\bf
n},n_{0})\vspace*{0.15cm}\\\
\pi_{a}^{(0)}\equiv\pi_{a}(n_{0})\vspace*{0.15cm}\\\ a^{(0)}\equiv
a(n_{0})\vspace*{0.15cm}\\\ \left\\{\tilde{\pi}_{\tilde{\varphi}_{\rm
m}}^{(0)}\right\\}\equiv\left\\{\tilde{\pi}_{\tilde{\varphi}_{\rm m}}({\bf
n},n_{0})\right\\}\vspace*{0.15cm}\\\ \left\\{\tilde{\varphi}_{\rm
m}^{(0)}\right\\}\equiv\left\\{\tilde{\varphi}_{\rm m}({\bf
n},n_{0})\right\\}\vspace*{0.15cm}\\\
\tilde{R}^{(0)}\equiv\tilde{R}(n_{0})\end{array}\right\\}~{}~{}~{}\Longrightarrow\left\\{\begin{array}[]{l}\Delta\tilde{\pi}_{\tilde{\phi}}^{(p)}=A_{p}\Delta\tilde{\pi}_{\tilde{\phi}}^{(p-1)}+\delta\tilde{\eta}\mathcal{K}_{\tilde{\phi}}\left[\tilde{\phi}^{(p-1)},\left\\{\tilde{\varphi}_{\rm
m}^{(p-1)}\right\\},\tilde{R}^{(p-1)}\right]\vspace*{0.15cm}\\\
\Delta\tilde{\phi}^{(p)}=A_{p}\Delta\tilde{\phi}^{(p-1)}+\delta\tilde{\eta}a^{(1)\,3-\alpha}\tilde{\pi}_{\tilde{\phi}}^{(p-1)}\vspace*{0.15cm}\\\
\Delta\pi_{a}^{(p)}=A_{p}\Delta\pi_{a}^{(p-1)}+\delta\tilde{\eta}\mathcal{K}_{a}\left[\tilde{R}^{(p-1)}\right]\vspace*{0.15cm}\\\
\Delta a^{(p)}=A_{p}\Delta
a^{(p-1)}+\delta\tilde{\eta}a^{(1)\alpha-1}\pi_{a}^{(p-1)}\vspace*{0.15cm}\\\
\\{\Delta\tilde{\pi}_{\tilde{\varphi}_{\rm
m}}^{(p)}\\}=\left\\{A_{p}\Delta\tilde{\pi}_{\tilde{\varphi}_{\rm
m}}^{(p-1)}+\delta\tilde{\eta}\mathcal{K}_{\tilde{\varphi}_{\rm
m}}\left[\tilde{\phi}^{(p-1)},\left\\{\tilde{\varphi}_{\rm
m}^{(p-1)}\right\\}\right]\right\\}\vspace*{0.15cm}\\\
\\{\Delta\tilde{\varphi}_{\rm
m}^{(p)}\\}=\left\\{A_{p}\Delta\tilde{\varphi}_{\rm
m}^{(p-1)}+\delta\tilde{\eta}\mathcal{D}_{\tilde{\varphi}_{\rm
m}}\left[\tilde{\phi}^{(p-1)},\left\\{\tilde{\varphi}_{\rm
m}^{(p-1)}\right\\}\right]\right\\}\vspace*{0.15cm}\\\
\tilde{\pi}_{\tilde{\phi}}^{(p)}=\tilde{\pi}_{\tilde{\phi}}^{(p-1)}+B_{p}\Delta\tilde{\pi}_{\tilde{\phi}}^{(p)}\vspace*{0.15cm}\\\
\tilde{\phi}^{(p)}=\tilde{\phi}^{(p-1)}+B_{p}\Delta\tilde{\phi}^{(p)}\vspace*{0.15cm}\\\
\pi_{a}^{(p)}=\pi_{a}^{(p-1)}+B_{p}\Delta\pi_{a}^{(p)}\vspace*{0.15cm}\\\
a^{(p)}=a^{(p-1)}+B_{p}\Delta a^{(p)}\vspace*{0.15cm}\\\
\left\\{\tilde{\pi}_{\tilde{\varphi}_{\rm
m}}^{(p)}\right\\}=\left\\{\tilde{\pi}_{\tilde{\varphi}_{\rm
m}}^{(p-1)}+B_{p}\Delta\tilde{\pi}_{\tilde{\varphi}_{\rm
m}}^{(p)}\right\\}\vspace*{0.15cm}\\\ \left\\{\tilde{\varphi}_{\rm
m}^{(p)}\right\\}=\left\\{\tilde{\varphi}_{\rm
m}^{(p-1)}+B_{p}\Delta\tilde{\varphi}_{\rm
m}^{(p)}\right\\}\vspace*{0.15cm}\\\
\tilde{R}^{(p)}=\tilde{R}\left[\tilde{\phi}^{(p)},\tilde{\pi}_{\tilde{\phi}}^{(p)},\\{\tilde{\varphi}_{\rm
m}^{(p)}\\},\\{\tilde{\pi}_{\tilde{\varphi}_{\rm
m}}^{(p)}\\}\right]\end{array}\right\\}_{p=1,\dots,s}$ (140)
$\displaystyle\Longrightarrow\left\\{\begin{array}[]{l}\tilde{\pi}_{\tilde{\phi}}({\bf
n},n_{0}+1)=\tilde{\pi}_{\tilde{\phi}}^{(s)}\vspace*{0.15cm}\\\
\tilde{\phi}({\bf n},n_{0}+1)=\tilde{\phi}^{(s)}\vspace*{0.15cm}\\\
\pi_{a}(n_{0}+1)=\pi_{a}^{(s)}\vspace*{0.15cm}\\\
a(n_{0}+1)=a^{(s)}\vspace*{0.15cm}\\\
\left\\{\tilde{\pi}_{\tilde{\varphi}_{\rm m}}({\bf
n},n_{0}+1)\right\\}=\left\\{\tilde{\pi}_{\tilde{\varphi}_{\rm
m}}^{(s)}\right\\}\vspace*{0.15cm}\\\ \left\\{\tilde{\varphi}_{\rm m}({\bf
n},n_{0}+1)\right\\}=\left\\{\tilde{\varphi}_{\rm
m}^{(s)}\right\\}\vspace*{0.15cm}\\\
\tilde{R}(n_{0}+1)=\tilde{R}^{(s)}\end{array}\right.\hskip 284.52756pt$ (148)

The final piece involves implementing Eq. 40, this can be explicitly checked
at every time step and provides a robust way to check the stability of the
algorithm. To implement the adaptive time step, one proceed as explained
above. In particular, the error $\Delta$ is computed as a sum over all the
fields and all lattice points of the type

$\displaystyle\Delta=B_{4}\sum_{x\in\Lambda}\sum_{f\in\left\\{\tilde{\pi}_{\tilde{\phi}},\tilde{\phi},\tilde{a},\tilde{\pi}_{\tilde{a}},\left\\{\tilde{\varphi}_{\rm
m}^{(p)}\right\\},\left\\{\tilde{\pi}_{\tilde{\varphi}_{\rm
m}}^{(p)}\right\\}\right\\}}|\Delta f(x)|\ .$ (150)

## References

  * (1) L. Parker and D. Toms, Quantum Field Theory in Curved Spacetime: Quantized Fields and Gravity.  Cambridge Monographs on Mathematical Physics, 2009.
  * (2) N. Birrell and P. Davies, Quantum Fields in Curved Space.  Cambridge Monographs on Mathematical Physics, 1984.
  * (3) M. Atkins and X. Calmet, Bounds on the Nonminimal Coupling of the Higgs Boson to Gravity, Phys. Rev. Lett. 110 (2013), no. 5 051301, [1211.0281].
  * (4) J. R. Espinosa, G. F. Giudice, and A. Riotto, Cosmological implications of the Higgs mass measurement, JCAP 05 (2008) 002, [0710.2484].
  * (5) K. Enqvist, T. Meriniemi, and S. Nurmi, Higgs Dynamics during Inflation, JCAP 07 (2014) 025, [1404.3699].
  * (6) J. R. Espinosa, G. F. Giudice, E. Morgante, A. Riotto, L. Senatore, A. Strumia, and N. Tetradis, The cosmological Higgstory of the vacuum instability, JHEP 09 (2015) 174, [1505.04825].
  * (7) Y. Ema, R. Jinno, K. Mukaida, and K. Nakayama, Violent Preheating in Inflation with Nonminimal Coupling, JCAP 02 (2017) 045, [1609.05209].
  * (8) K. Enqvist, M. Karciauskas, O. Lebedev, S. Rusak, and M. Zatta, Postinflationary vacuum instability and Higgs-inflaton couplings, JCAP 11 (2016) 025, [1608.08848].
  * (9) D. G. Figueroa, A. Rajantie, and F. Torrenti, Higgs field-curvature coupling and postinflationary vacuum instability, Phys. Rev. D 98 (2018), no. 2 023532, [1709.00398].
  * (10) T. Markkanen, A. Rajantie, and S. Stopyra, Cosmological Aspects of Higgs Vacuum Metastability, Front. Astron. Space Sci. 5 (2018) 40, [1809.06923].
  * (11) B. Horn, The Higgs field and early universe cosmology: a (brief) review, MDPI Physics 2 (2020), no. 3 503–520, [2007.10377].
  * (12) A. Mantziris, T. Markkanen, and A. Rajantie, Vacuum decay constraints on the Higgs curvature coupling from inflation, JCAP 03 (2021) 077, [2011.03763].
  * (13) O. Lebedev, The Higgs portal to cosmology, Prog. Part. Nucl. Phys. 120 (2021) 103881, [2104.03342].
  * (14) D. H. Lyth and A. Riotto, Particle physics models of inflation and the cosmological density perturbation, Phys. Rept. 314 (1999) 1–146, [hep-ph/9807278].
  * (15) A. Riotto, Inflation and the theory of cosmological perturbations, ICTP Lect. Notes Ser. 14 (2003) 317–413, [hep-ph/0210162].
  * (16) B. A. Bassett, S. Tsujikawa, and D. Wands, Inflation dynamics and reheating, Rev. Mod. Phys. 78 (2006) 537–589, [astro-ph/0507632].
  * (17) A. D. Linde, Inflationary Cosmology, Lect. Notes Phys. 738 (2008) 1–54, [0705.0164].
  * (18) D. Baumann, Inflation, in Theoretical Advanced Study Institute in Elementary Particle Physics: Physics of the Large and the Small, pp. 523–686, 2011.  0907.5424.
  * (19) A. De Felice and S. Tsujikawa, f(R) theories, Living Rev. Rel. 13 (2010) 3, [1002.4928].
  * (20) A. A. Starobinsky, A New Type of Isotropic Cosmological Models Without Singularity, Phys. Lett. B 91 (1980) 99–102.
  * (21) Planck Collaboration, Y. Akrami et al., Planck 2018 results. X. Constraints on inflation, Astron. Astrophys. 641 (2020) A10, [1807.06211].
  * (22) BICEP, Keck Collaboration, P. A. R. Ade et al., Improved Constraints on Primordial Gravitational Waves using Planck, WMAP, and BICEP/Keck Observations through the 2018 Observing Season, Phys. Rev. Lett. 127 (2021), no. 15 151301, [2110.00483].
  * (23) F. L. Bezrukov and M. Shaposhnikov, The Standard Model Higgs boson as the inflaton, Phys. Lett. B 659 (2008) 703–706, [0710.3755].
  * (24) Planck Collaboration, Y. Akrami et al., Planck 2018 results. X. Constraints on inflation, Astron. Astrophys. 641 (2020) A10, [1807.06211].
  * (25) R. Allahverdi, R. Brandenberger, F.-Y. Cyr-Racine, and A. Mazumdar, Reheating in Inflationary Cosmology: Theory and Applications, Ann. Rev. Nucl. Part. Sci. 60 (2010) 27–51, [1001.2600].
  * (26) 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].
  * (27) K. D. Lozanov, Lectures on Reheating after Inflation, 1907.04402.
  * (28) R. Allahverdi et al., The First Three Seconds: a Review of Possible Expansion Histories of the Early Universe, 2006.16182.
  * (29) B. A. Bassett and S. Liberati, Geometric reheating after inflation, Phys. Rev. D 58 (1998) 021302, [hep-ph/9709417]. [Erratum: Phys.Rev.D 60, 049902 (1999)].
  * (30) S. Tsujikawa, K.-i. Maeda, and T. Torii, Resonant particle production with nonminimally coupled scalar fields in preheating after inflation, Phys. Rev. D 60 (1999) 063515, [hep-ph/9901306].
  * (31) S. Tsujikawa, K.-i. Maeda, and T. Torii, Preheating with nonminimally coupled scalar fields in higher curvature inflation models, Phys. Rev. D 60 (1999) 123505, [hep-ph/9906501].
  * (32) C. Fu, P. Wu, and H. Yu, Nonlinear preheating with nonminimally coupled scalar fields in the Starobinsky model, Phys. Rev. D 99 (2019), no. 12 123526, [1906.00557].
  * (33) F. Bezrukov, D. Gorbunov, and M. Shaposhnikov, On initial conditions for the Hot Big Bang, JCAP 06 (2009) 029, [0812.3622].
  * (34) J. Garcia-Bellido, D. G. Figueroa, and J. Rubio, Preheating in the Standard Model with the Higgs-Inflaton coupled to gravity, Phys. Rev. D 79 (2009) 063531, [0812.4624].
  * (35) D. G. Figueroa, Preheating the Universe from the Standard Model Higgs, AIP Conf. Proc. 1241 (2010), no. 1 578–587, [0911.1465].
  * (36) E. I. Sfakianakis and J. van de Vis, Preheating after Higgs Inflation: Self-Resonance and Gauge boson production, Phys. Rev. D 99 (2019), no. 8 083519, [1810.01304].
  * (37) M. P. DeCross, D. I. Kaiser, A. Prabhu, C. Prescod-Weinstein, and E. I. Sfakianakis, Preheating after Multifield Inflation with Nonminimal Couplings, I: Covariant Formalism and Attractor Behavior, Phys. Rev. D 97 (2018), no. 2 023526, [1510.08553].
  * (38) M. P. DeCross, D. I. Kaiser, A. Prabhu, C. Prescod-Weinstein, and E. I. Sfakianakis, Preheating after multifield inflation with nonminimal couplings, II: Resonance Structure, Phys. Rev. D 97 (2018), no. 2 023527, [1610.08868].
  * (39) M. P. DeCross, D. I. Kaiser, A. Prabhu, C. Prescod-Weinstein, and E. I. Sfakianakis, Preheating after multifield inflation with nonminimal couplings, III: Dynamical spacetime results, Phys. Rev. D 97 (2018), no. 2 023528, [1610.08916].
  * (40) R. Nguyen, J. van de Vis, E. I. Sfakianakis, J. T. Giblin, and D. I. Kaiser, Nonlinear Dynamics of Preheating after Multifield Inflation with Nonminimal Couplings, Phys. Rev. Lett. 123 (2019), no. 17 171301, [1905.12562].
  * (41) Y. Watanabe and E. Komatsu, Reheating of the universe after inflation with f(phi)R gravity, Phys. Rev. D 75 (2007) 061301, [gr-qc/0612120].
  * (42) L. H. Ford, Gravitational Particle Creation and Inflation, Phys. Rev. D 35 (1987) 2955.
  * (43) B. Spokoiny, Deflationary universe scenario, Phys. Lett. B 315 (1993) 40–45, [gr-qc/9306008].
  * (44) T. Damour and A. Vilenkin, String theory and inflation, Phys. Rev. D 53 (1996) 2981–2989, [hep-th/9503149].
  * (45) P. J. E. Peebles and A. Vilenkin, Quintessential inflation, Phys. Rev. D 59 (1999) 063505, [astro-ph/9810509].
  * (46) M. Peloso and F. Rosati, On the construction of quintessential inflation models, JHEP 12 (1999) 026, [hep-ph/9908271].
  * (47) G. Huey and J. E. Lidsey, Inflation, brane worlds and quintessence, Phys. Lett. B 514 (2001) 217–225, [astro-ph/0104006].
  * (48) A. S. Majumdar, From brane assisted inflation to quintessence through a single scalar field, Phys. Rev. D 64 (2001) 083503, [astro-ph/0105518].
  * (49) K. Dimopoulos and J. W. F. Valle, Modeling quintessential inflation, Astropart. Phys. 18 (2002) 287–306, [astro-ph/0111417].
  * (50) C. Wetterich, Variable gravity Universe, Phys. Rev. D 89 (2014), no. 2 024005, [1308.1019].
  * (51) C. Wetterich, Inflation, quintessence, and the origin of mass, Nucl. Phys. B 897 (2015) 111–178, [1408.0156].
  * (52) M. W. Hossain, R. Myrzakulov, M. Sami, and E. N. Saridakis, Variable gravity: A suitable framework for quintessential inflation, Phys. Rev. D 90 (2014), no. 2 023512, [1402.6661].
  * (53) J. Rubio and C. Wetterich, Emergent scale symmetry: Connecting inflation and dark energy, Phys. Rev. D 96 (2017), no. 6 063509, [1705.00552].
  * (54) D. G. Figueroa and E. H. Tanin, Inconsistency of an inflationary sector coupled only to Einstein gravity, JCAP 10 (2019) 050, [1811.04093].
  * (55) C. Caprini and D. G. Figueroa, Cosmological Backgrounds of Gravitational Waves, Class. Quant. Grav. 35 (2018), no. 16 163001, [1801.04268].
  * (56) T. J. Clarke, E. J. Copeland, and A. Moss, Constraints on primordial gravitational waves from the Cosmic Microwave Background, JCAP 10 (2020) 002, [2004.11396].
  * (57) M. Giovannini, Stochastic backgrounds of relic gravitons: a theoretical appraisal, PMC Phys. A 4 (2010) 1, [0901.3026].
  * (58) L. A. Boyle and A. Buonanno, Relating gravitational wave constraints from primordial nucleosynthesis, pulsar timing, laser interferometers, and the CMB: Implications for the early Universe, Phys. Rev. D 78 (2008) 043531, [0708.2279].
  * (59) D. G. Figueroa and E. H. Tanin, Ability of LIGO and LISA to probe the equation of state of the early Universe, JCAP 08 (2019) 011, [1905.11960].
  * (60) T. Opferkuch, P. Schwaller, and B. A. Stefanek, Ricci Reheating, JCAP 07 (2019) 016, [1905.06823].
  * (61) N. Bernal, A. Ghoshal, F. Hajkarim, and G. Lambiase, Primordial Gravitational Wave Signals in Modified Cosmologies, JCAP 11 (2020) 051, [2008.04959].
  * (62) D. G. Figueroa and C. T. Byrnes, The Standard Model Higgs as the origin of the hot Big Bang, Phys. Lett. B 767 (2017) 272–277, [1604.03905].
  * (63) K. Dimopoulos and T. Markkanen, Non-minimal gravitational reheating during kination, JCAP 06 (2018) 021, [1803.07399].
  * (64) D. Bettoni, A. Lopez-Eiguren, and J. Rubio, Hubble-induced phase transitions on the lattice with applications to Ricci reheating, 2107.09671.
  * (65) D. G. Figueroa, A. Florio, F. Torrenti, and W. Valkenburg, The art of simulating the early Universe – Part I, JCAP 04 (2021) 035, [2006.15122].
  * (66) J. Williamson, Low-storage runge-kutta schemes, Journal of Computational Physics 35 (1980), no. 1 48–56.
  * (67) A. Bazavov and T. Chuna, Efficient integration of gradient flow in lattice gauge theory and properties of low-storage commutator-free Lie group methods, 2101.05320.
  * (68) D. G. Figueroa, A. Florio, F. Torrenti, and W. Valkenburg, CosmoLattice, 2102.01031.
  * (69) O. Bertolami, P. Frazao, and J. Paramos, Reheating via a generalized non-minimal coupling of curvature to matter, Phys. Rev. D 83 (2011) 044010, [1010.2698].
  * (70) R. Kallosh and A. Linde, Universality Class in Conformal Inflation, JCAP 07 (2013) 002, [1306.5220].
  * (71) L. E. Parker and D. Toms, Quantum Field Theory in Curved Spacetime: Quantized Field and Gravity.  Cambridge Monographs on Mathematical Physics. Cambridge University Press, 8, 2009.
  * (72) M. Carpenter and C. Kennedy, Third-order 2n-storage runge-kutta schemes with error control, 1994.
  * (73) M. Carpenter and C. Kennedy, Fourth-order 2n-storage runge-kutta schemes, 1994\.