Lattice Simulations of Non-minimally Coupled Scalar Fields in the Jordan Frame

The presence of scalar ﬁelds with non-minimal gravitational interactions of the form ξ | φ | 2 R may have important implications for the physics of the early universe. We propose a procedure to solve the dynamics of non-minimally coupled scalar ﬁelds directly in the Jordan frame, where the non-minimal couplings are maintained explicitly. Our algorithm can be applied to lattice simulations that include minimally coupled ﬁelds and an arbitrary number of non-minimally coupled scalars, with the expansion of the universe sourced by all ﬁelds present. This includes situations when the dynamics become fully inhomogeneous, fully non-linear (due to e.g. backreaction or mode rescattering eﬀects), 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 ﬁeld when the inﬂaton oscillates following the end of inﬂation.

The presence of scalar fields with non-minimal gravitational interactions of the form ξ|φ| 2 R may have important implications for the physics of the early universe. We propose a procedure to solve the dynamics of non-minimally coupled scalar fields directly in the Jordan frame, where the nonminimal couplings are maintained explicitly. Our algorithm can be applied to lattice simulations 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.

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 φ (either singlet or charged), one can add to the action an operator of the form ξ|φ| 2 R, where R is the Ricci scalar and ξ 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 [1,2], where it is a running parameter that cannot be set to zero at all energy scales. 1 In the case of the Standard Model (SM) Higgs field Φ, the operator ξ|Φ| 2 R is actually the only missing operator of dimension-4 that respects all symmetries of the SM and gravity. The coupling ξ 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, ξ 10 15 [3]. 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. [4][5][6][7][8][9][10][11][12][13].
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. [14][15][16][17][18]). 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. [19]). If f (R) = 0 and f (R) = 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 [20], defined by f (R) = R + αR 2 with α > 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 ξ|φ| 2 R lead to inflationary predictions in excellent agreement with current observational constraints [21,22]. 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 a daniel.figueroa@ific.uv.es b adrien.florio@stonybrook.edu c toby.opferkuch@cern.ch d bestef@physik.uzh.ch 1 Exceptionally, the running vanishes for the conformal value ξ = 1/6. Higgs [23]. It is very interesting that data from cosmological observations clearly favours plateau-like potentials that naturally emerge in these scenarios [24].
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 [25][26][27][28]). This can occur in non-minimally coupled inflaton χ scenarios, where preheating into other degrees of freedom, e.g. scalar fields {φ}, can be realized very efficiently when couplings of the form g 2 χ 2 φ 2 or ξφ 2 R are considered. Preheating scenarios considering a simple monomial inflaton potential and a preheat scalar field nonminimally coupled to gravity ξφ 2 R were first considered in [29] and later on in [30], where an inflaton-preheat field coupling g 2 χ 2 φ 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. [29], and we keep that terminology here. Preheating following inflation due to higher order curvature terms f (R) = R + α n R n has been studied considering geometric preheating effects in [31,32]. Preheating after Higgs-Inflation was originally studied in [33][34][35], considered in more detail in [7,36], and lastly for modified setups, like R 2 -Higgs inflation, in [37]. Preheating in multi-field inflationary scenarios with N scalar fields {φ j } and couplings ξ j φ 2 j R has also been extensively studied [38][39][40][41][42]. Finally, preheating with a non-minimal gravitational coupling f (φ)R with f a general function of φ, has been also considered, see e.g. [43]. All of the above instabilities due to the presence of non-minimal gravitational couplings can be generically regarded as "gravitational reheating" mechanisms 2 , 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 [44] (see also [45]). Namely, a massless scalar field φ non-minimally coupled to gravity is excited towards the end of inflation. 3 The inflaton potential is chosen such that there is a sustained kination dominated era with stiff equation of state 1/3 < w ≤ 1 after inflation. As a consequence, the energy stored in φ (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 [47][48][49][50][51][52][53][54][55]. The original gravitational reheating mechanism, however, was shown in Ref. [56] to be inconsistent with BBN/CMB constraints [57,58] due to an excess amount of gravitational wave production 4 . 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 [44,45,47]), but rather experiences a tachyonic instability due to the change in sign of R ∝ (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. [64], 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 [62], see also [65,66]. In [62], 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. 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 is non-singular. 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 Ω 2 = 1 − ξ(φ/m p ) 2 which appears to be non-invertible for φ 2 = m 2 p /ξ and ξ > 0 (here m p 2.4 × 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 Ω 2 (φ(x)), with φ(x) often treated as a quantum field. Some works evaluate the map using the vacuum expectation value φ 2 , but then it is not clear that the Einstein frame description fully captures the physics of the theory originally written in the Jordan frame, especially in the case where the initial conditions are determined purely through quantum fluctuations.
In this paper, we introduce a technique for solving the system directly as written in the original Jordan frame, avoiding the need to perform any conformal transformation. In particular, we are able to solve the dynamics of an arbitrary scalar field φ with a non-minimal coupling to gravity ξφ 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 with η an "α-time" variable related to cosmic time by dt = a(η) α dη. Here α 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) where primes indicate derivatives with respect to η. 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 {ϕ m } minimally coupled to gravity, together with a scalar field φ non-minimally coupled to gravity. Without loss of generality, the action of this system reads where 1 2 m 2 p R is the standard Hibert-Einstein term, 1 2 ξRφ 2 represents a non-minimal gravitational interaction of φ, and V (φ, {ϕ m }) encompasses both the self interactions of φ as well as its non-gravitational interactions with the minimally-coupled matter sector. The term L 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 φ, we obtain the following equation of motion for φ where = g µν ∇ µ ∇ ν and ∇ µ 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 φ. Using the α-time metric given in Eq. (1), the above equation becomes Equivalently, we can think of this gravitational interaction as part of an effective potential that includes all together the non-minimal coupling to gravity, the non-gravitational interactions with the minimallycoupled matter sector, as well as the self-interactions of φ. For convenience, we can then think of a Lagrangian for φ given by The Einstein equations are obtained by varying Eq. (3) with respect to g µν with where T m µν and T φ µν have been defined as the energy-momentum tensors of the minimally-coupled matter fields and the non-minimally coupled scalar field φ, respectively. In particular, one finds (see Appendix B) the energy-momentum tensor of the non-minimally coupled field to be The trace of T φ µν , defined as T φ = g µν T φ µν , 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 where G = g µν G µν = (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 where V ,φ = ∂V /∂φ. Notably, if ξ = 1/6, then for V = 0 or V ∝ φ 4 we find that T φ µν is traceless, i.e. T φ = 0, as a consequence of the conformal invariance of S φ = d 4 x √ −gL φ in these cases. Given the FLRW metric in Eq. (1), the consistency of the Einstein equations requires that T µν takes the form of the energy-momentum tensor of a perfect fluid T µ ν = diag {−ρ(η), p(η), p(η), p(η)}. We note that while fields can develop large spatial inhomogeneities, the homogeneous and isotropic pressure and energy density p(η) and ρ(η) 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 µν 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 α-time where we defined H = a /a, which is related to the cosmic time Hubble rate H as H = H/a α . We define the energy density and pressure as with . . . denoting volume averages. With these definitions, the explicit expressions for the energy density and pressure of the non-minimally coupled field ρ φ and p φ are found to be 5 5 The volume averages of the total divergence terms ∇ 2 φ 2 = ∇ · ∇φ 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.
In principle, one can solve for the scale factor a(η) 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 Inserting the expression for T φ 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 This expression for R can be directly related to the evolution of the scale factor using Eq. (2). This leads to the differential equation 6 that together with the equation of motion for φ 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φ with f * some typical field amplitude and ω * a characteristic (inverse) time scale of the problem to be studied. The choice of f * and ω * 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. [68] for examples). If the matter sector simply comprises of a set of scalar fields {ϕ m }, these are normalized as in Eq. (22). We note that rescaling the coordinates by ω * naturally induces the following rescaling in R It is also natural to introduce rescaled energy densities and pressure Next, we reduce the order of the equation of motion forφ by introducing a conjugate momentum variable as 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 [68]. If the matter sector is comprised of scalar fields, we simply introduce a set of conjugate momenta {π ϕ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 functionalK φ , as follows 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(η) as we arrive to a system of coupled first-order differential equations depending on another kernel functional, To close the system, an expression forR is needed in both kernels K φ , K a . Using Eq. (20), we can writẽ where we have used T m = 3p m − ρ m and introduced the volume-averaged kineticẼφ K and gradientẼφ G energy densitiesẼφ In summary, Eqs. (26) and (28), together with the expression forR in Eq. (29) (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 ξ i φ 2 i in Eq. (29).

III. LATTICE FORMULATION
In order to evolve our system of equations Eqs. (26), (28) and (29) 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 δx = L/N . We work with finite differences and use the following notation for the forward and backward derivatives where f is an arbitrary scalar function defined on the lattice sites n = (n 1 , n 2 , n 3 ), andî 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 We are now in a position to define the evolution equations by introducing the discrete kernels which we have already written in terms of natural field and spacetime variables, c.f. Eq. (22). We have also introduced dimensionless discrete derivatives∇ given by Eq. (31) in terms of the dimensionless lattice spacing δx =L/N = ω * δx, withL = ω * L. At this point, it is important to realize thatR =R[φ,π φ , {φ m }, {π ϕm }] depends on all fields and conjugate momentum variables, and hence the kernel for the non-minimally coupled fieldφ 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. [68] 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. [69,70]. 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 (the resulting convergence is shown in Appendix D). In terms of rescaled variables, it reads where we have dropped the ∇ 2 φ 2 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 CosmoLattice [68,71], which can perform user-friendly and versatile field theory simulations. These new algorithms will be made publicly available in a future update of CosmoLattice.

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 φ 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 [29] illustrated in Fig. 1. The fact that R becomes oscillatory can be seen from the traced Einstein equations, assuming that the homogeneous inflaton field χ initially dominates the energy density of the universe, such that One illustrative example is the case where V inf (χ) = 1 2 m 2 χ 2 , in which case R can be approximated in cosmic time (α = 0) as In this expression, it is manifest that R oscillates between positive and negative values due to the harmonic oscillations of χ. In general, for an inflaton potential with a minimum around the origin and an arbitrary power law behavior V inf (χ) ∝ |χ| 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 2 φ,eff = ξR, will still alternate periodically between positive and negative values. As a consequence of the periodic tachyonic stages (m 2 φ,eff < 0), initial quantum vacuum fluctuations of the spectator field φ can be exponentially amplified if the strength of its non-minimal coupling ξ is large enough. The amplification may persist until the effective tachyonic mass of φ is fully screened by its own self-interactions, or until the energy of φ 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 A a linear analysis of the initial instability of the mode functions of φ, then in Section IV B we present an analysis of the evolution of the system once the dynamics become non-linear.

A. Initial Conditions via Linear Analysis
Our procedure will consist of computing the power spectrum of the φ 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 χ and a light spectator field φ with a non-minimal coupling to gravity, similar to Refs. [29,31,32,72] with In this theory, the inflaton χ and the spectator field φ interact only gravitationally through the non-minimal coupling ξRφ 2 . During slow-roll inflation, we have a quasi de-Sitter phase where R ≈ 12H 2 and H ≈ constant (its time derivative is slow-roll suppressed). This means that for ξ > 0, the spectator field has a heavy effective mass m 2 φ,eff ≈ 12ξH 2 during inflation. We assume that this effective mass dominates over V (φ), such that the potential V (φ) can be neglected during inflation. This, combined with the fact that the non-minimally coupled spectator field φ is energetically subdominant during inflation, justifies the use of a linear analysis. It will be convenient to work in conformal time (α = 1) where the metric is conformally flat and quantization proceeds as in Minkowski space. In that case, we can write the action for φ in terms of the canonically normalized field ϕ = aφ We then canonically quantize ϕ aŝ where [â k ,â † k ] = (2π) 3 δ(k − k ) and the modes are normalized such that ϕ k ϕ * k − ϕ k ϕ * k = i. The mode functions ϕ k (τ ) obey the equation of motion given by Eq. (42) in momentum space We assume that the evolution of each mode starts far inside the Hubble radius, namely −kτ 1, where the curvature is negligible. In that case, the modes should approach the Bunch-Davies vacuum We are interested in the power spectrum ∆ ϕ (k, τ ), which in terms of the two-point function is defined as The power spectrum of the original field φ is then related as For our numerical results (which involve integrating over many e-folds of inflaton), it proves easier to solve Eq. (44) in cosmic time (α = 0), where it readsφ with the Bunch-Davies initial condition now expressed as where we have used τ ≈ −1/(aH), given that H changes very slowly. As previously mentioned, the non-minimally coupled spectator field φ 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 with the evolution of the homogenous inflaton governed by the standard Klein-Gordon equation of motion in cosmic timeχ In our numerical analysis we consider an observationally viable inflationary model inspired by α-attractors [73], with inflaton potential parametrised as [21] V inf (χ) = Λ 4 tanh p cχ m p , with p = 4, 6 , which flattens out for |χ| m p /c (where c is a dimensionless parameter), and takes a power-law form V ∝ χ p for |χ| m p /c. We take c = 0.1 which reproduces the observed value of the scalar perturbations at CMB scales for V inf (χ CMB ) = (1.6 × 10 16 GeV) 4 and saturates the upper bound on the scale of inflation [21] corresponding to Λ = 1.79 × 10 16 GeV.
We solve Eq. (48) 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) = β, with β 1 a penetration factor. Larger values of β better approximate the Bunch-Davies initial condition, but also increase simulation time, so as a compromise we choose β = 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 β inside the horizon, we require approximately ∆N log 10 3 β ≈ 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 χ i such that we obtain 14 e-folds of inflation (this corresponds to χ i = 9.23(10.87)m p for p = 4(6), respectively). Following this procedure, we numerically integrate Eqs. (48), (50) and (51) 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 ξ = 50.
For comparison, our numerical results for φ can be compared to the predicted power spectrum during inflation in pure de-Sitter space, which was computed analytically in Ref. [62] as Here z = k/(aH * ), µ 2 = 12(ξ − 3/16), and H * is taken to be the Hubble rate at the end of inflation. The approximate equality holds for ξ > 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 ∝ k 2 , indicating that they are not excited. The transition between power laws occurs around k/(aH * ) ≈ µ, where we have µ ≈ √ 12ξ for ξ 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. (48), as the tachyon is regularized (on a per mode basis) when k/(aH) ≈ |ξR/H 2 | ≈ O(1) √ ξ √ 6ξ, where we used that tachyonic growth happens for −6 ≤ R/H 2 < 0 and assumed ξ 1. Correspondingly, the peak in the power spectrum, which clearly emerges about ∼1 e-fold after the end of inflation, is observed at k/(aH) ≈ O(1) √ ξ. 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.

B. Non-linear Lattice Analysis
After a clear peak in the power spectrum emerges in the linear analysis, but still before any interactions of φ or its backreaction onto the background dynamics are relevant, we move to solving the system on the lattice. In particular, we treat φ as a classical field whose initial fluctuations are drawn randomly from a Gaussian distribution with power spectrum given by ∆ φ 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 where Eq. (56) determines the background evolution as discussed in Section II and R is defined as with F (φ) given in Eq. (20). We adapt the continuum equations to the lattice following Section III and Appendix C. All simulations with λ = 0 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 ω * = H i , where H i is the Hubble rate at the start of the linear analysis 14 e-folds before the end of inflation. Note that we are using time-steps of H i δt = 0.01 in the evolution with k IR /H i = 2.5 × 10 −3 (4 × 10 −3 ) for ξ = 10 (50 or 100). In the case λ = 1 × 10 −5 , we used lattices of size N = 512 and k IR /H i = 1 × 10 −2 . Let us first consider the case where V (φ) = 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 ≤ R/H 2 ≤ 12. While R is negative, the non-minimally coupled field has a tachyonic effective mass and can experience exponential growth. The structure of this growth can be seen in the power spectrum of the left-hand panel of Fig. 3, where different line colors indicate the number of efolds post inflation. Here we see the growth of the peak agrees well between the linear (dashed) and lattice (solid) results so long that the energy density in the NMC field is sub-dominant. Departure from the linear analysis can be more easily seen in the right-hand panel of Fig. 4 where we see the expectation value to overcome the expansion of the universe (as in the ξ = 50, 100 cases) since there we neglected the contribution of φ to the background evolution. This is not the case in the lattice analysis, where for these large values of ξ the backreaction of φ on the background evolution asymptotically drives ξ φ 2 /m 2 p 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ξ − 1) φ 2 drives initially R to a large positive value, as represented by the green (blue) spikes for ξ = 100 (50) in the left panel of Fig. 4. This results in a large positive effective mass squared ξR, that induces a restoring force for φ, opposing its growth. The field velocity φ 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 Parameter values for the inflationary potential are c = 0.1, Λ 1.79 × 10 16 GeV, 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 ξ 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 ξ but for p = 6, where the energy density of the inflaton redshifts faster than radiation.
spike, R oscillates with a small damped amplitude, so the successive tachyonic mass stages cannot overtake anymore the expansion of the universe, and ξ φ 2 /m 2 p approaches asymptotically a constant value. Returning to the power spectrum (left-hand panel of Fig. 3) at N = 2 differences in both the peak and also the UV tail of the spectrum arise. The origin of the additional structure in the peak of the lattice results is simply resultant from the Ricci scalar remaining positive once the backreaction occurs. Subsequently the NMC just behaves as a free oscillator and is no longer driven. If we define the equation of state in terms of the total pressure and energy densities w = p(η)/ρ(η), then w can be written in terms of R by combining Eqs. (2), (13) and (14) where H = H/a α . 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 ρ is always constrained to be positive definite by Eq. (13), it can be seen from Eq. (58) that R/H 2 > 12 corresponds to a large negative pressure, namely p < −ρ, which violates the classical energy conditions due to the non-minimal interaction of φ with the gravitational field. In particular, the dominant contribution to the pressure during these spikes comes from the (1 − 4ξ) φ 2 term in Eq. (18), which is always negative for ξ > 1/4. Before turning to the evolution of the energy density, a comment on the effect of a quartic interaction term in the scalar potential is in order. In the right-hand panel of Fig. 3 we show the power spectrum for the case where ξ = 100 with a quartic λ = 10 −5 . As expected re-scattering of modes leads to additional power in the UV spectrum while also screening the effect of the tachyon, namely the growth of the peak is diminished already at N = 1.1.
Finally, in Fig. 6 we show the evolution of the energy density of both the inflaton χ (red) and the non-minimally coupled spectator field φ (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 ξ = {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 (φ) = 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 (φ). We consider the cases of V (φ) = m 2 φ 2 /2 with mass m/H i = 0.02 (thin purple line) as well as V (φ) = λφ 4 /4 with a quartic λ = 10 −5 (thin gray line). The case V (φ) = m 2 φ 2 /2 exhibits similar behavior to the V (φ) = 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 ξ 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 φ in curved spacetime, be it the SM Higgs or otherwise, inevitably acquires a non-minimal coupling to gravity of the form ξ|φ| 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, however, the equivalence of these two frames in situations where the initial conditions are set by quantum fluctuations is unclear.
In contrast, in this work we have developed an approach to solve the dynamics of non-minimally coupled scalars in an expanding universe directly as written 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 CosmoLattice [68,71] 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 φ 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 φ. While the effective squared mass is negative, tachyonic growth of the non-minimally coupled field can occur if the value of ξ 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 φ reaches an O(1) value of the total energy density within an 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 ξ 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, 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 CosmoLattice, which will allow anyone to perform user-friendly and versatile cosmological simulations involving non-minimally coupled scalar fields.

ACKNOWLEDGMENTS
We would like to thank Nicolás Loayza for useful numerical assistance and for validating a number of results presented here. In addition, TO would also like to thank Valerie Domcke, Alan Guth, David I. Kaiser, 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 α-time Given a flat "α-time" FLRW background described by the metric the Christoffel symbols can be computed from The non-vanishing Christoffel symbols are found to be where primes indicate derivatives with respect to η. From this, the components of the Ricci tensor are found to be From this, the Ricci scalar R = g µν R µν can be computed directly. We find (A6) Appendix B: Energy-momentum tensor of a non-minimally coupled scalar The energy momentum tensor for the non-minimally coupled φ is defined in Eq. (9) as The remaining variation reads A useful identity one can show by explicit computation from the expression of R in terms of the metric is the following for any scalar function f and ∇ ν is the covariant derivative associated to g µν . Applying this to Eq. (B2), we obtain and putting it altogether we get where we have defined ≡ ∇ σ ∇ σ = g ρσ ∇ ρ ∇ σ . We would like to compute the trace of this energy-momentum tensor for φ. For a moment, let us write down the result working in d + 1 dimensions. Then, we find where G = g µν G µν = (1 − d)R/2 is the trace of the Einstein tensor with respect to the metric. This expression can be further simplified using φ 2 = 2φ φ + 2∂ µ φ∂ µ φ and the equation of motion Eq. (5) for φ which gives φ φ = ξRφ 2 + φV ,φ . This leads to the following expression Note that the coefficient of ∂ µ φ∂ µ φ + ξRφ 2 vanishes for ξ = (d − 1)/(4d), which is indeed the conformal value of ξ in d + 1 dimensions. Now setting d = 3, we get which indeed leads to Eq. (12) for T φ = g µν T φ µν quoted in the main text. Let us comment on some specific cases of Eq. (B7) which reproduce known results. Consider the conformal value for ξ in d + 1 = 4 dimensions, ξ = 1/6. Then, we find which vanishes identically in the scaleless cases of V = 0 or V ∝ φ 4 , as expected. For a quadratic potential V = m 2 φ 2 /2, we get which reproduces the result of Ref. [74].

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. [68]. Consider a vector x(t) of M -variables x(t) = (x 1 (t), . . . , x M (t)) T and a system of first order differential equations of the typė Then, a RK method of order s is characterized by a one-step iteration of the type with This iteration effectively split the time interval δt into s subintervals δt = s i=1 c i δ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 Explicit RK methods have the property b ij = 0 for all j ≥ 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: (C7) 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 storelve (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 [69] (see also Ref. [70] for a recent application in lattice QCD). This method hinges on rewriting of Eqs. (C2) to (C4) as with ∆ y (i) = A i ∆ y (i−1) + δt , K y (i−1) , 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. [69,70] 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: RK2ME : We have also implemented the following third order method from Ref. [75] 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: Finally, fourth order 2N -storage schemes also exist, we refer the interested reader to Refs. [70,76] 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 y (3) is already at second order accuracy in δt. At the extra memory cost of saving the previous solution x n in case the update fails, one can then easily turn it into an adaptive time-step RK scheme. As reviewed in Refs. [70,75], this can be achieved by estimating the distance ∆ in some norm between y (3) and y (4) . If this distance is smaller than some requested tolerance , the update is accepted; if not it is rejected and the step is repeated. In both situation, the time step is updated to δt new = 0.95 · ∆ 1/3 δt old . (C14) This updated always decrease δ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 ∆ is to compute the Euclidean distance between the solutions ∆ = | y (3) − y (4) | = B 4 |∆ y (4) | (C15) with | y| = L i=1 y 2 l for the L-component vector y = (y 1 , . . . , y L ) T . Note that the efficiency of such an adaptive scheme varies from model to model and needs to be studied on a case-by-case basis.
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: ∆φ, ∆πφ, {∆φ m }, {∆πφ m }, ∆a, ∆π a . We can then implement a generic s-stage 2N -storage RK method as follows The final piece involves implementing Eq. (36), 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 ∆ is computed as a sum over all the fields and all lattice points of the type