Preheating in Full General Relativity

TL;DR

This paper explores the role of full general relativity in the non-linear dynamics of preheating after inflation, using numerical methods to assess gravitational effects and their impact on early universe processes.

Contribution

It introduces three numerical methods for studying preheating in full general relativity, extending beyond perturbation theory and classical assumptions.

Findings

Perturbation theory fails for horizon-scale modes.

Preheating remains viable without black hole formation.

Full relativity does not necessarily disrupt preheating dynamics.

Abstract

We investigate the importance of local gravity during preheating, the non-linear dynamics that may be responsible for starting the process of reheating the universe after inflation. We introduce three numerical methods that study a simple preheating scenario while relaxing gravitational assumptions, culminating in studying the process in full numerical relativity. We confirm that perturbation theory is no longer valid when one considers modes whose wavelengths are comparable to the size of the horizon at the end of inflation; however, this breakdown does not necessarily lead to a breakdown of the preheating process in non-linear gravity. For the specific model we test we find no evidence for the creation of primordial black holes from the instabilities in this model. Finally, we remark on the opportunity for future numerical study of non-linear gravitational dynamics in the early universe.

Figures

Equations

$${\mathcal{L}}_{\mathrm{m}}=\left(\frac{1}{2}{\partial }^{\mu }{\phi }_I{\partial }_{\mu }{\phi }_I\right)-V({\phi }_I)$$

$$V({\phi }_I)=\frac{1}{2}{m}^2{\phi }_1^2+\frac{1}{2}{g}^2{\phi }_1^2{\phi }_2^2.$$

$$⟨{\left|{\phi }_I(k)\right|}^2⟩=\frac{1}{2\omega }$$

$$g_{\mu \nu }^{\mathrm{FLRW}}=\mathrm{diag}(-1,a^2(t),a^2(t),a^2(t))$$

$$H^2\equiv {\left(\frac{\dot{a}}{a}\right)}^2=\frac{8\pi }{3m_{\mathrm{pl}}^2}⟨\rho ⟩$$

$$\rho \equiv -T_0^0=V({\phi }_I)+\sum _I\left(\frac{1}{2}{\dot{{\phi }_I}}^2+\frac{{\left(\nabla {\phi }_I\right)}^2}{2a^2}\right).$$

$$\ddot{{\phi }_I}+3H\dot{{\phi }_I}-\frac{{\nabla }^2{\phi }_I}{a^2}+\frac{\partial V}{\partial {\phi }_I}=0.$$

$$\ddot{{\phi }_1}$$

$$\ddot{{\phi }_2}$$

$$\mathrm{Var}({\phi }_{\mathrm{I}})=\sqrt{⟨{\phi }_I^2⟩-⟨{\phi }_I{⟩}^2},$$

$$g^{\rm CPT}_{\mu\nu}=\left(\begin{array}[]{cc}-(1+2\Phi)&0\\ 0&a^{2}(t)(1-2\Psi)\end{array}\right)$$

$$\rho =⟨\rho ⟩+\delta \rho ,$$

$${\nabla }^2\Phi -3H{a}^2\left(\dot{\Phi }+H\Phi \right)=4\pi G{a}^2\delta \rho$$

$${\nabla }^2\left(\dot{\Phi }+H\Phi \right)=\sum _I\frac{4\pi }{m_{\mathrm{pl}}^2}{\delta }^{ij}{\partial }_i\left({\partial }_j{\phi }_I{\partial }_0{\phi }_I\right).$$

$${\ddot{\phi }}_I=$$

$$+(4\dot{\Phi }-3H){\dot{\phi }}_I-(2\Phi +1)\frac{\partial V}{\partial {\phi }_I}$$

$$g^{\rm BSSN}_{\mu\nu}=\left(\begin{array}[]{cc}-\alpha^{2}+\beta_{l}\beta^{l}&\beta_{j}\\ \beta_{i}&e^{4\phi}\bar{\gamma}_{ij}\end{array}\right).$$

$${\partial }_t\varphi =$$

$${\partial }_t{\bar{\gamma }}_{ij}=$$

$$+{\bar{\gamma }}_{kj}{\partial }_i{\beta }^k-\frac{2}{3}{\bar{\gamma }}_{ij}{\partial }_K{\beta }^k$$

$${\partial }_{t}K=$$

$$+4\pi \alpha \left(\rho +S\right)+{\beta }^i{\partial }_{i}K$$

$${\partial }_t{\tilde{A}}_{ij}=$$

$$+\alpha \left(K{\tilde{A}}_{ij}-2{\tilde{A}}_{il}{\tilde{A}}_j^l\right)+{\beta }^k{\partial }_k{\tilde{A}}_{ij}$$

$$+{\tilde{A}}_{ik}{\partial }_j{\beta }^k+{\tilde{A}}_{kj}{\partial }_i{\beta }^k-\frac{2}{3}{\tilde{A}}_{ij}{\partial }_k{\beta }^k$$

$${\partial }_t{\bar{\Gamma }}^i=$$

$$-8\pi {\bar{\gamma }}^{ij}{S}_j+6{\tilde{A}}^{ij}{\partial }_j\varphi )+{\beta }^j{\partial }_j{\bar{\Gamma }}^i$$

$$-{\bar{\Gamma }}^j{\partial }_j{\beta }^i+\frac{2}{3}{\bar{\Gamma }}^i{\partial }_j{\beta }^j+\frac{1}{3}{\bar{\gamma }}^{li}{\partial }_l{\partial }_j{\beta }^j$$

$$+{\bar{\gamma }}^{lj}{\partial }_j{\partial }_l{\beta }^i$$

$$\mathcal{H}\equiv 0=$$

$$-\frac{e^{5\varphi }}{12}{K}^2+2\pi e^{5\varphi }\rho ,$$

$${\mathcal{M}}^i\equiv 0=$$

$$S^i=-{\gamma }^{ij}{n}^a{T}_{aj},$$

$$\left({\partial }_t-{\beta }^j{\partial }_j\right)\alpha =-2\alpha K.$$

$${\partial }_t\alpha =-2\alpha \left(K-⟨K⟩\right).$$

$${\partial }_t{\beta }^i$$

$${\partial }_t{B}^i$$

$${\Pi }_I\equiv$$

Full text

I Introduction

If inflation occurred in the early Universe, it must have come to an end—at least locally. While measurements of the cosmic microwave background and large scale structure constrain inflationary models during inflation, the final moments of the inflationary period and the subsequent preheating of the Universe may provide clues for how inflation fits into our models of high energy physics. The lack of a unique mathematical model of inflation is a compelling reason to search for testable predictions from the inflationary and preheating periods. Most models of inflation employ a field (or fields) whose homogeneous value(s) determines the final dynamics of this epoch; in many of these, the field ends up oscillating about a minimum of the potential. Frequently in these scenarios, nonlinear processes take over to accelerate the decay of the inflaton field in a family of scenarios known as preheating, see, e.g., Ref. Amin:2014eta for a review.

While the processes of the field sectors of preheating have been examined closely throughout the past decades, the effects of gravity during preheating have not yet been extensively studied. At the same time, there exist some exciting possibilities from the preheating period. Large, nonthermal, and nonlinear inhomogeneities are characteristic of this period and detailed numerical study is needed to understand the role of this physics.

In this paper we aim to do a first investigation of nonlinear gravity during the violent process of preheating in the early Universe.

To date, cosmological perturbation theory (CPT) has been an effective and predictive way to study departures from a purely homogeneous and isotropic, Friedmann-LeMaître-Robertson-Walker (FLRW), universe. All previous work studying gravity during and after preheating exists in this scenario. The first was done using DEFROST, Frolov:2008hy , where it was shown that local gravitational effects seem to be small, an observation that has been confirmed in other numerical simulations. More recent work continues to study the generation of linear gravitational perturbations during preheating, Huang:2011gf ; Lozanov:2019ylm , where gravitational effects show promise in helping to understand the rich phenomenology of this period.

However, CPT is inherently limited, by construction, to a linearized treatment of gravity and is not able to resolve the strong-field regime. Equally important is that linearized gravity, by definition, does not allow for gravitational modes to couple, hiding potentially important aspects of the preheating process. As such, the search for new physics beyond the perturbative regime strongly motivates the application of full general relativity (GR) to the study of cosmology. Applying nonlinear gravity to cosmological scenarios is a recent advancement in high resolution and accurate numerical simulations. The formalisms employed to do full numerical relativity Pretorius:2005gq ; Campanelli:2005dd ; Baker:2005vv ; Garfinkle:2003bb ; Baumgarte:1998te ; Shibata:1995we have been recently applied to late universe scenarios Giblin:2015vwq ; Bentivegna:2015flc ; Giblin:2016mjp ; Bentivegna:2016stg ; Macpherson:2016ict ; Giblin:2017juu ; East:2017qmk ; Wang:2018qfr ; Macpherson:2018akp ; Macpherson:2018btl ; Giblin:2018ndw as well as preinflationary scenarios East:2015ggf ; East:2016anr ; Clough:2017efm ; Bloomfield:2019rbs and oscillons Muia:2019coe . Other work along these lines include adding gravitational effects for scalar metric perturbations, beyond linear order, BasteroGil:2007mm ; BasteroGil:2010nm during the preheating period. Other works include fully nonlinear gravity, either in one dimension Easther:1999ws or in a restricted gauge in Ref. BasteroGil:2011zza .

In CPT, the linearized Einstein’s equations couple spatial derivatives of the Newtonian gravitational potential $\Phi$ to local inhomogeneities in the energy density, $\delta\rho$. When inhomogeneities grow, the gravitational response may be large.

The idea that preheating can amplify gravitational degrees of freedom dates back to early work Nambu:1996gf ; Taruya:1997iv where the authors noticed that gravitational modes can be excited via parametric instabilities. This work was significantly extended by Refs. Bassett:1998wg ; Bassett:1999mt ; Bassett:1999ta ; Tsujikawa:2002nf , where the authors set out to understand where linear perturbation theory breaks down and point out that nonlinear gravity is necessary to study the gravitational effects of preheating Bassett:1999mt .

More recently, it has been noted that overdensities on fixed subhorizon scales might collapse due to nonlinear gravitational effects Jedamzik:2010dq ; Jedamzik:2010hq —a generic feature of preheating which could lead to the formation of primordial black holes, e.g., Ref. Martin:2019nuw . Additionally, it was shown in Ref. Easther:2010mr that an instability exists for modes near the Hubble scale at the end of inflation. These modes provide the dominant contribution to the density contrast $\delta(t,\vec{x})=\delta\rho(t,\vec{x})/\left<\rho\right>$, causing it to become large. This may signal the breakdown of linear perturbation theory and the need for a nonlinear treatment.

Our primary goal here is to increase the robustness of gravitational approximations during the preheating stage, culminating in studying the problem in full numerical relativity. We begin with a standard implementation of the Grid and Bubble Evolver GABE in an FLRW universe before including localized linear gravity in Newtonian gauge. We then introduce GABERel, an adaptation of GABE Child:2013ria that numerically evolves the full set of Einstein’s field equations on an expanding background. We utilize GABERel to evolve nonlinear inflaton modes in the postinflationary Universe during parametric resonance preheating. This provides a precise numerical treatment of the nonlinear effects that lead to the breakdown of coherent oscillations in the postinflationary Universe and demonstrates the ability of GABERel to move beyond limitations of linear perturbation theory. We confirm prior work that predicts the breakdown of perturbation theory at the end of infation; however, we go further to show that full numerical relativity safely controls these inhomogeneities which do not end up being catastrophic: by either causing the coherent modes of the inflaton to break down or by seeding primordial black holes.

We also present the formalism one can use to calculate the (first-order) gauge-invariant Bardeen potentials in the framework of full GR. We use this as a test of GABERel’s robustness and to compare these potentials to the Newtonian simulation in an equivalent preheating context. We also use GABERel to evolve a single black hole of mass $M$ to assess its ability to resolve strong-field dynamics to demonstrate the accuracy of this new technique.

This paper is organized as follows. In Sec. II we introduce our model of preheating and its numerical implementation in three simulations with increasingly comprehensive treatments of gravitational dynamics. In Sec. III we outline our initialization procedure and give computational details. We present our results in Sec. IV and a discussion in Sec. V. Appendix A contains information regarding how to calculate the Bardeen potentials in a fully general-relativistic simulation. Appendix B outlines methods by which we validate the accuracy of GABERel’s gravitational dynamics.

II Model

We are primarily interested in the effects of gravity on preheating; hence, we study a canonical model of preheating with two scalar fields: the inflaton, $\varphi_{1}$, and a coupled, massless, matter field, $\varphi_{2}$. For convenience we use upper-case latin letters, $I$, $J$, etc, as a field index to simplify notation when we need to sum over fields; we use $\varphi_{I}$ when we need to generally refer to one or all of the scalar fields and repeated field indices imply a sum. The matter Lagrangian for our system is

[TABLE]

with a potential

[TABLE]

Throughout this work, we take the coupling to be $g^{2}=2.5\times 10^{-7}$ and $m=10^{-6}\,\,m_{\mathrm{pl}}$ as a toy model that has been studied extensively in the literature Traschen:1990sw . This is an interesting model both because it has been repeatedly used as a benchmark preheating scenario and is widely recognized, but also because it exhibits broad resonance—see, e.g., Fig. 3 in Ref. Amin:2014eta where linear analysis predicts efficient preheating. Since the dynamics of this model are so well known, it is an ideal case to study as a first try even though this inflationary model is disfavored Aghanim:2018eyx and this specific model of preheating needs extensions to completely deplete the energy in the inflaton, e.g., Ref. Dufaux:2006ee . In this model, when inflation ends the vast majority of the energy density of the Universe is trapped in the inflaton condensate. As the homogeneous mode of the inflaton field oscillates at the minimum of the potential it parametrically amplifies modes of the matter field until the $\varphi_{2}$ particles scatter into the $\varphi_{1}$ field and nonlinear physics leads to the breakdown of the condensate. We want to be careful here to delineate the different phases of the process. The first of these phases, I, is where the $\varphi_{2}$ field is being amplified, but the inflaton condensate is mostly unaffected by any backscattering. The next phase, II, is characterized by backscattering onto $\varphi_{1}$ and an amplification of the variances of both the $\varphi_{1}$ and $\varphi_{2}$ fields. The final phase, III, sees dramatic nonlinear amplifications of both fields; this phase is important as significant power is transferred between different modes. When studying this process on a finite grid, it is important to note that we almost always reach a stage at which significant power is transferred to the smallest-resolvable scales. We will note this moment in the following sections as it becomes more dangerous as the dynamics become more complicated.

Since we are interested in gravitational effects, and want to be able to include any effects that arise when $\ddot{a}=0$, we start our simulations one $e$-folding before inflation ends. In this model, that corresponds to a (homogeneous) field value of $\phi_{0}\approx 0.415\,m_{\rm pl}$ and velocity $\dot{\phi}_{0}\approx-0.154\,\,m_{\mathrm{pl}}$. Inflation ends when the field passes $\phi_{0}\approx 0.201\,\,m_{\mathrm{pl}}$. Throughout the work here we denote the subscript “0” to refer to quantities evaluated at the beginning of the simulation, i.e. we take $a_{0}=e^{-1}$, and the subscript “$*$” will refer to quantities evaluated at the end of inflation, $a_{*}=1$.

In all simulations, we will initialize our scalar field modes to be in the Bunch-Davies vacuum,

[TABLE]

where $\omega=\sqrt{k^{2}+m_{\rm eff}^{2}}$. In an expanding universe, this approximation is only good for modes (well) inside the horizon—however there are known adaptations that work for horizon-sized modes Frolov:2008hy . Since we start our simulations before the end of inflation, it allows the entire (or most of the) box to be subhorizon on the initial slice. In doing so we can trust our initial conditions while still allowing the fields to evolve to be superhorizon by the time inflation ends—knowing that there could be instabilities on the scale of $k/a\sim H_{*}$.

In the next three subsections we present three increasingly comprehensive treatments of the gravitational physics in this model. First we introduce the canonical treatment in a rigidly expanding spacetime, the FLRW limit. We then study the problem in perturbation theory in Newtonian gauge in CPT, followed by an introduction of the fully nonlinear methods using the BSSN formalism Baumgarte:1998te ; Shibata:1995we .

II.1 The FLRW limit

In general, studies of preheating are done in the FLRW regime, where the line element is

[TABLE]

and the dynamics of the scale factor are determined by the 00-component of Einstein’s equations, Friedmann’s equation,

[TABLE]

where the $\left<\ldots\right>$ denote an average over a constant-time hypersurface, and

[TABLE]

As usual each of the scalar fields evolve according to the Klein Gordon equation,

[TABLE]

To be comprehensive, Eq. (7) implies similar—but crucially different—equations of motion for the inflaton and matter field in this preheating model,

[TABLE]

Figure 1 illustrates the importance of these differences in this canonical process for a simulation with $L_{*}=5\,m^{-1}$. While this box is slightly superhorizon at the end of inflation it becomes subhorizon very quickly after; shaded regions in this figure show the three phases of preheating as described above. To track the production of particles across the box, we use the variance of the fields,

[TABLE]

as a way to assess inhomogeneity as a function of time. In the FLRW regime, Friedmann’s equation contains all of the gravitational physics and there is no local gravitational collapse; the fields simply evolve in an expanding background, preheating the Universe.

II.2 Cosmological perturbation theory

To include the effects of local gravity, it is natural to next look to cosmological perturbation theory. In this limit we introduce a scalar perturbation to Eq. (4), in our case in Newtonian gauge Bardeen:1980kt ; Ma:1995ey ; Weinberg:2008zzc ,

[TABLE]

The metric perturbations $\Phi$ and $\Psi$ are the Newtonian gravitational potentials. They are equal when the scalar anisotropic stress vanishes, an assumption we will make. While this is not strictly true for a scalar field source, it is a good approximation when the stress-energy tensor is almost diagonal. We also note that we are not choosing to use the (first-order) gauge invariant Mukhanov-Sasaki variables Mukhanov:1990me , as is commonly done in CPT; here we wish to keep the metric and field perturbations separate with an eye toward going beyond the linear analysis. The energy density can be decomposed at linear order into a background value and a perturbation,

[TABLE]

so that we can study Einstein’s equations order by order. The 00 component yields Friedmann’s equation, Eq. (4), at zeroth order, and

[TABLE]

at first order. It is common to drop the second term of Eq. (13) in the nonrelativistic limit; however, since our sources are scalar fields, we can calculate the second term by simultaneously solving the divergence of the $0i$ terms of Einstein’s equations, $\delta^{ij}\partial_{i}G_{0i}=8\pi/m^{2}_{\rm pl}\,\delta^{ij}\partial_{i}T_{0i}$,

[TABLE]

The effects of local gravity enters into the equations of motion for the scalar fields, $\partial_{\mu}\left(\sqrt{-g}g^{\mu\nu}\partial_{\nu}\varphi_{I}\right)=\sqrt{-g}\,dV/d\varphi_{I}$, resulting in

[TABLE]

In practice we use a spectral method to solve Eqs. (13) and (14) simultaneously. We take the Fourier transform of the right-hand side of these two equations and then invert the Laplacians in momentum space to get the Fourier transform of $\Phi$ as well as the combination $(\dot{\Phi}+H\Phi)$; we then inverse Fourier transform these to get $\Phi$ and $\dot{\Phi}$ in configuration space.

Local gravity during preheating has been studied in Refs. Frolov:2008hy ; BasteroGil:2007mm ; BasteroGil:2010nm , where the size of the Newtonian potential was found to be small. On the other hand, an examination of Eq. (13) shows that the magnitude of the Newtonian potential depends on the volume of the simulation. As the wavelength of a density perturbation approaches to the Hubble scale $k\rightarrow H$, the resulting size of the gravitational response $\Phi$ grows, even if $\delta\rho/\langle\rho\rangle$ is constant. Therefore allowing density perturbations to exist at near-horizon scales (by including more modes in the simulation) amplifies the size of the Newtonian potential.

This is the case in Ref. Easther:2010mr , where there is a mechanism through which the Universe becomes inhomogeneous at the Hubble scale at the end of preheating. In this scenario inflaton modes with $k\sim a_{*}H_{*}$ contribute significantly to the extrinsic curvature as well as the density contrast, eventually spoiling linearity. We demonstrate this effect in Sec. IV, Fig. 6 by including these modes in the CPT simulation; after enough time in the phase of coherent oscillations, roughly $\sim 2$-$3$ $e$-foldings, these modes grow nonlinearly and the validity of the perturbative treatment breaks down. This breakdown, discussed at length in Sec. IV, necessitates a fully general-relativistic treatment of the postinflationary Universe’s gravitational dynamics, raising the possibility of primordial gravitational wave production during preheating as well as the nonlinear growth of density perturbations in a strong gravitational regime Martin:2019nuw .

II.3 Fully nonlinear gravity

The BSSN Formalism—an adaptation of the ADM formalism of GR Arnowitt:1959ah ; Baumgarte:1998te ; Shibata:1995we —utilizes a 3+1 decomposition of the field equations, with the metric parametrized as

[TABLE]

The variables $\alpha$ and $\beta^{i}$ are called the lapse and the shift and parametrize gauge degrees of freedom. The spatial metric $\gamma_{ij}$ is conformally rescaled, $\gamma_{ij}=e^{4\phi}\bar{\gamma}_{ij}$, such that the conformal metric has unit determinant $\mathrm{det}\left|\bar{\gamma}_{ij}\right|=1$. This notation will be general so that quantities written with overbars are related to the conformal spatial metric, $\bar{\gamma}_{ij}$. Then, the field equations can be expressed as a system of first order differential equations in these metric variables,

[TABLE]

In addition to these evolution equations, the BSSN formalism improves numerical stability by defining an auxiliary variable $\bar{\Gamma}^{i}\equiv\bar{\gamma}^{jk}\bar{\Gamma}^{i}_{jk}=-\partial_{j}\bar{\gamma}^{ij}$, which is evolved independently,

[TABLE]

Using the auxiliary variables effectively promotes the “divergence” of the conformal metric to an independent variable so that its second spatial derivatives do not need to be directly calculated via lattice stencils. This promotion makes the problem hyperbolic which, hopefully, causes numerical noise to remain bounded in the simulation. We see this same technique in a moment when we discuss the scalar field evolution.

The $0\mu$ components of the field equations furnish four constraint equations which must be satisfied at all times throughout the system’s evolution to ensure that the numerical simulation remains a valid solution. These read

[TABLE]

in which

[TABLE]

where the vector $n^{a}=(\alpha^{-1},-\alpha^{-1}\beta^{i})$ is normal to the spatial hypersurface.

Satisfying these Hamiltonian, Eq. (II.3), and momentum Eq. (23), constraints—or showing that they are bounded—is a central challenge of performing fully relativistic simulations, which we discuss in Appendix B.

Because the lapse and shift are gauge variables, we are free to define how they evolve; this process fixes a slicing of spacetime which makes it possible to choose the coordinates best fit for different problems. One common choice, geodesic slicing or synchronous gauge, fixes $\alpha=1$ and $\beta^{i}=0$, so an observer at fixed spatial coordinates travels along a geodesic. Another extensively studied choice for the evolution of $\alpha$ is a family of coordinate systems called $1+\log$ slicing

[TABLE]

Here we almost exclusively employ a variant of Eq. (25), where the extrinsic curvature is replaced by its deviation from its average on our spatial hypersurfaces $\langle K\rangle$,

[TABLE]

When using Eq. (26), we will use the hyperbolic Gamma-driver condition for the shift,

[TABLE]

The combination of Eqs. (25), (27), and (28) have proven to be highly robust choices suitable for black hole simulations. The parameter $\eta$ is a constant often related to the total mass in the simulation Baumgarte:2010ndz . In this work we find good numerical results with $\eta\sim 50\,m$.

Following previous work Baumgarte:2010ndz ; Clough:2015sqa , the scalar field equations are split into two first-order differential equations where the momentum of the fields are replaced with $\Pi_{I}$,

[TABLE]

As mentioned earlier when discussing $\bar{\Gamma}^{i}$ in Eq. (II.3), we also promote the spatial derivative of the scalar field to an independent variable, $\psi_{iI}\equiv\partial_{i}\varphi_{I}$, which evolves according to

[TABLE]

In the following section, we discuss how we initialize this system of fully general-relativistic equations of motion and how we implement the evolution on a finite grid.

III Numerical Simulations

III.1 Initial conditions

The initial conditions must satisfy the constraint equations, Eqs. (II.3) and (23), and be physically motivated. We initialize our parametric resonance scenario in correspondence with the standard semiclassical, linearized treatment of this model’s postinflationary dynamics. This choice will facilitate a comparison between GABERel and well-known results from CPT, allowing us to analyze the ability of a fully general-relativistic simulation to evolve the nonlinear dynamics and breakdown of coherent oscillations predicted in Ref. Easther:2010mr .

As stated in Sec. II, the mean power in each mode is set by the Bunch Davies vacuum, Eq. (3). To realize such initial conditions, we randomly generate Gaussian-distributed power spectra for each component of each field’s Fourier mode. The amplitude of each Fourier mode is then drawn from a Rayleigh distribution

[TABLE]

where $\sigma^{2}=L^{3}/4\pi^{4}\omega_{k}$. A random realization of Eq. (32) is then inverse Fourier transformed, giving us an initial field configuration in configuration space.

In order to connect our results with CPT we begin by linearizing the BSSN variables, setting $\bar{\gamma}_{ij}=\delta_{ij}+h_{ij}$, $\tilde{A}_{ij}=0+a_{ij}$, $\phi=0+\delta\phi$, $K=K_{0}+\delta K$, $\alpha=1+\delta\alpha$, and $\beta^{i}=0+\delta\beta^{i}$. We set $h_{ij}=0$, so that gravitational waves do not enter at first order. This choice fixes $\bar{R}=0$. We also choose $\beta^{i}=0$ and $a_{ij}=0$.

To draw a correspondence with Newtonian Gauge perturbation theory, we also have

[TABLE]

and

[TABLE]

Linearizing the Hamiltonian constraint, Eq. (II.3), makes it clear that $K_{0}=-3H$. Working with the Fourier transform of our variables (denoted by overtildes) this constraint is then equivalent to satisfying

[TABLE]

which should look identical to the CPT expression, Eq. (13), if we make the identification that

[TABLE]

The choice $a_{ij}=0$ reduces the momentum constraint to

[TABLE]

which is equivalent to the CPT expression, Eq. (14), with the same identification, Eq. (36).

These identifications mean that any solution to the set of perturbation theory equations, Eqs. (13) and (14), also satisfy the linearized constraints, Eqs. (II.3) and (23), so long as we use the identifications, Eqs. (33), (34) and (36). Of course, this is only strictly true at the linear level, but is sufficient for the simulations here, as can be seen in Appendix B.

III.2 Other numerical details

GABERel solves differential equations in full numerical relativity on a finite, expanding lattice under periodic boundary conditions. In this work, our lattice resolution is $N^{3}$ at $N=64$ and we use a time step $\Delta t=\Delta x/20$ (where $\Delta x=L_{*}/N$ is the initial lattice spacing) with one exception, see Fig. 8. Choosing the time step to be a small fraction of the initial lattice spacing provides good resolution and puts us far from the regime in which causality becomes an issue.

Simulations on a finite grid eventually reach a stage at which significant power is transferred to the smallest resolvable scales (i.e., the Nyquist frequency, $f_{\rm{Nyq}}=2\sqrt{3}\pi N/L$). This is dangerous because the difficulty in resolving such modes leads to a gradual increase in numerical inaccuracy, eventually spoiling the validity of the simulation. One way we alleviate this problem is by initializing this system without very much power on these small scales, thereby delaying the onset of this power growth. We do this by smoothly suppressing power in modes with frequencies larger than the cutoff frequency $f_{c}$, defined by the dimensionless parameter $\xi_{c}\equiv f_{c}/f_{\rm{Nyq}}$, which we always choose to be less than 1. This is done by filtering the power in all modes by a window function on the initial slice,

[TABLE]

Here, $\kappa$ is a parameter between zero and one that determines the sharpness of the window function; in this work we set $\kappa=0.75$ and $\xi_{c}=1/8$. In this way we lower power in the smallest resolvable scales, which do not participate in the initial stages of preheating and make it difficult to satisfy the Hamiltonian and momentum constraints.

IV Results

We begin by examining the effect local gravitational physics has on preheating. In this case (as in Ref. Frolov:2008hy ) we expect to see growth of the Newtonian potential during the preheating stages. We also expect that the statistics of the Newtonian potential will depend on which modes we are able to resolve. As we discussed in Sec. II.2, when the density contrast, $\delta\rho/\left<\rho\right>(k)$ is amplified at the same level, modes near the Hubble scale create larger gravitational deviations from homogeneity than do smaller wavelength disturbances. Figure 2 confirms this by showing the maximum and minimum values of the Newtonian potential as a function of time for three different box sizes. Larger box sizes include longer-wavelength modes which create larger maximum and minimum values of the gravitational potential.

The three box sizes we study in Fig. 2 represent three distinct regimes: the smallest of these sizes, $L_{*}=2\,m^{-1}$, is always subhorizon, that is the longest resolvable wavelength mode is bigger than $H^{-1}$ throughout the simulation. The second case, $L_{*}=5\,m^{-1}$, represents a marginal case where the longest wavelength mode becomes superhorizon just at the end of inflation and then quickly retreats inside the Hubble radius. The third, $L_{*}=11\,m^{-1}$, has a long-wavelength mode that stays smaller than $H^{-1}$ for a few oscillations of the scalar field. We note that the largest of these boxes, $L_{*}=11\,m^{-1}$ has an initial superhorizon mode which mildly breaks our assumptions of Eq. (3). Figure 3 shows a comparison of the physical size ratio of the smallest resolvable wave vector, $k_{\rm min}=2\pi/aL_{*}$ compared to the Hubble scale $H$.

For the three box sizes featured in Figs. 2 and 3, we expect that nonlinear gravity should have a minimal effect on the simulations. The size of the Newtonian potential does not significantly grow, and we expect to only see small changes in the fields resulting from different realizations of the initial vacuum states. This is a good regime, then, in which to validate GABERel and confirm that no nonlinear gravitational effects enter on these scales.

To achieve this validation, we look to see if we can calculate the (first-order) gauge-invariant metric perturbations. Generically, scalar perturbations to the metric can be written as Weinberg:2008zzc

[TABLE]

In terms of these variables, the Bardeen potentials are (see, e.g., Ref. Baumann:2009ds )

[TABLE]

The translation from the generic form of the metric perturbation to these new variables is straightforward, but computationally expensive and we leave the details on how to calculate these from the BSSN variables to Appendix A (see Ref. Giblin:2017ezj for a similar treatment). In CPT, where we have assumed no scalar anisotropic stress, $\Phi=\Phi_{B}=\Psi_{B}$, so we would expect the statistics of these gauge-invariant potentials to be identical between our simulations. To test this, we show a comparison of the power spectra of $\Phi$ (from a CPT simulation) to the dimensionless power spectra, $k^{3}\mathcal{P}_{\Phi}$, as defined by

[TABLE]

of $\Phi_{B}$ and $\Psi_{B}$ from a BSSN simulation in Fig. 4. When using a $1+\log$ slicing, we see very good agreement in these quantities in the regions where we anticipate the results to coincide, despite the radically different methods used to calculate these quantities. The very lowest bins are closest to the horizon where we anticipate CPT and BSSN to disagree, while at higher frequencies nonlinear dynamics mode mix in the BSSN simulations and we see greater power in those simulations. On small scales, where the power is negligible, differences are due in part to resolving small differences in comparatively large quantities. We also note that our two situations also differ in that our CPT simulations explicitly ignore scalar anisotropic stress whereas the BSSN ones do not.

To complete our comparison in Fig. 4, we also run our code in geodesic slicing. This slicing condition is known to violate the linearized Einstein equation Giblin:2018ndw much more quickly than $1+\log$ slicing, which is confirmed by disagreement in high-frequency modes in Fig. 4. This discrepancy is a sign that the geodesic slicing simulations move out of perturbation theory on some scales faster than those in $1+\log$ slicing. Finally, we note that we have to take extreme care when analyzing the output of the BSSN simulations (see Appendix B). In this code, when power is moved to the Nyquist frequency, we see that the Hamiltonian and momentum constraints, Eqs. (II.3) and (23), are no longer satisfied. We take care not to derive physical meaning from these simulations after this time. For all of the BSSN simulations we present here, the constraints start to grow around $t=110\,m^{-1}$ and become unsatisfactorily satisfied before $t=130\,m^{-1}$.

Additionally, one can look at the statistics of the two fields to see agreement between our simulations. Figure 5 shows a comparison of the variances of the two scalar fields—and confirms that there are no notable differences between the simulations we run across box sizes and methods. At some box sizes, geodesic slicing alone exhibits early oscillations in the inflaton variance during the phase of coherent oscillations. This demonstrates the importance of slicing conditions when comparing fully relativistic results to those from CPT and FLRW treatments.

We can now move on to see if there are indications of stronger gravitational instabilities if we include more horizon-sized modes for a longer time in the simulation. For this we consider a box size of $L_{*}=20$, so that we have more modes near the Hubble scale for a longer portion of the simulation. We note that our initial conditions are also slightly inconsistent in this scenario, i.e., a small number of these modes that are not well within the horizon and are not exactly Bunch-Davies, Eq. (3). Figure 6 shows the ratio of the physical size of the lowest-frequency mode to the Hubble scale for a box of size $L_{*}=20$ along with a comparison to the same quantity for $L_{*}=11$.

We see in this case that the lowest frequency mode stays outside the horizon for more of the simulation, allowing many modes to cross into the horizon during the preheating process. These modes do become inhomogenous as the simulation proceeds and, as they are entering the horizon, cause the Newtonian potential to grow to order unity and the simulation to crash. As predicted in Ref. Easther:2010mr , these additional modes become nonlinear during the phase of coherent oscillations; the Universe may undergo a number of $e$-foldings before the contribution to the extrinsic curvature from the shift becomes comparable to that from the homogeneous evolution. This amplification is exactly the effect we see in our CPT simulations.111The authors of Ref. Easther:2010mr work in spatially flat gauge where the shift contributes to the Bardeen potentials. We then rely on comparing the (first-order) gauge invariant potentials either $\Phi$ (from our CPT simulations) or $\Phi_{B}$ and $\Psi_{B}$—to analyze this predicted breakdown of linearity. BSSN simulations, on the other hand, do not rely on the assumptions of CPT, and, hence, the code does not crash at this point. Figure 7 shows a comparison between the CPT simulation and the BSSN simulations for $L_{*}=20\,m^{-1}$ which demonstrates that fully nonlinear gravity is able to resolve the moderate density contrasts that exist as these modes are entering the horizon. It is important to note, though, that when we move to full GR, we do not see a breakdown of the coherent oscillations of the inflation field, despite the fact that the Newtonian potential (or the Bardeen potentials) grows to be close to order unity.

The bottom panel of Fig. 7 shows that statistics of the lapse, $\alpha$ for this simulation. This quantity is important when probing the existence of strong gravity and the possible creation of primordial black holes. In the simulations we have run here we see minimal departure in the lapse, and hence, do not see the creation (or seeds) of primordial black holes. This statement is, however, a model-dependent statement and not a generic prediction for all models of preheating.

We can also look towards the density contrast, $\delta=\delta\rho/\left<\rho\right>$, to assess the onset of nonlinear gravitational physics. This provides an intuitive way of assessing the validity of linear perturbation theory and a consistency check with Fig 6. CPT is built on the assumption that $\delta$ is small. In Fig. 8 we compare two-dimensional slices of $\delta$ throughout a CPT and a BSSN simulation for $L_{*}=20\,m^{-1}$. For both simulations, the density contrast initially grows slowly and remains small, and CPT remains valid. However, around $t\approx 120\,m^{-1}$—right after $\Phi_{k_{\rm{min}}}$ enters the horizon—$\delta$ grows sharply to order 1 in the CPT and BSSN simulations, with the growth being only marginally more pronounced in the CPT simulation. The CPT simulation fails at this point. This confirms the breakdown of linearity predicted in Refs. Bassett:1998wg ; Bassett:1999mt ; Bassett:1999ta ; Tsujikawa:2002nf (and reinforced in Ref. Easther:2010mr ) and shown in Fig. 6: once the horizon-sized modes enter the horizon during preheating, the density contrast quickly exceeds unity and the linear treatment is no longer valid, necessitating a fully general-relativistic treatment. Despite this nonlinear instability, we do not observe strong gravitational collapse or the formation of PBH during this stage of preheating. The growth of the density contrast moves out of the regime of validity of (first-order) CPT after this point, but does not necessarily signal unbounded growth of overdensities. However, neither the CPT nor BSSN simulation is valid shortly after this point. The former because $\Phi$ becomes large and the simulation crashes, the latter because power is moved to the Nyquist frequency and the constraints are no longer satisfied. Still, there is no indication that the gravitational instabilities are large enough, in this model, to continue to collapse. Nevertheless, thanks to the robustness of full numerical relativity, GABERel is able to successfully resolve these nonlinear physical processes that are inaccessible to CPT.

V Discussion

We have explored the effect of adding local gravitational effects to a toy model of preheating and present a rigorous analysis of gravitational effects after inflation. For this model of preheating, we have shown that there are no unexpected results that arise when local gravity—in either perturbation theory or in full GR—is included in the dynamics of the Universe. When local gravity is included, we see the excitation of scalar perturbations; in these cases, the Newtonian potential can get as large as $\Phi\sim 10^{-1}$, however, the subsequent field dynamics do not allow for these perturbations to become nonlinear and there is no indication of compact structures forming. While the Newtonian potential does grow, it does not seem that this is enough to cause nonlinearity in the gravitational part of the system—even when the Newtonian potential becomes order unity. The variations of the Newtonian potential are on large scales and do not seem to cause gravitational collapse.

We have also introduced a new computational tool, GABERel: software that can evolve scalar fields in full numerical relativity. We have validated that this software reproduces expectations from perturbation theory in its regime of validity and can be extended to regions in which perturbation theory fails. Using this tool for a canonical model of preheating shows that nonlinear gravitational effects stabilize the instabilities present in perturbation theory. We see that the phenomena of preheating in this setup shares a background evolution with the unperturbed FLRW analysis for the entirety of the simulation, even when departures grow beyond first-order perturbation theory. It also shows that some choices of coordinates (slicing) do a better job satisfying the linearized Einstein’s equations Weinberg:2008zzc ; Giblin:2018ndw than do others.

Finally, it is important to note that the lack of the formation of collapsed structures is not a generic result. In addition to previous CPT results hinting towards linear gravitational physics, the model of preheating studied here, Eq. (2), is not known to cause large metric perturbations. We envision that stronger instabilities (such as superhorizon or tachyonic instabilities Bassett:1999cg ; Bassett:2000ha ; Finelli:2000ya ; Felder:2000hj ; Felder:2001kt ; GarciaBellido:2001wn ; Kofman:2001rb ; GarciaBellido:2001cb ; Copeland:2002ku ; Suyama:2004mz ; Suyama:2006sr ; Barnaby:2006cq ; Felder:2006cc ; DiazGil:2007qx ; DiazGil:2008tf ; Dufaux:2008dn ; Battefeld:2009xw ; Barnaby:2009wr ; Abolhasani:2009nb ; Braden:2010wd ; Antusch:2015vna ; Adshead:2016iae ; DeCross:2016cbs ; Tranberg:2017lrx ; Adshead:2017xll ) that are closer to the horizon would be ideal candidates for more dramatic gravitational physics, and we look forward to studying more of these models in the future.

VI Acknowledgements

We are extremely grateful for a number of enlightening conversations as this work was being completed. We are indebted to Thomas Baumgarte and James Mertens for discussion of the BSSN formalism and helping with code validation techniques. We are also grateful for many conversations with Peter Adshead, Leia Barrowes, Robert R. Caldwell, David I. Kaiser, and Zachary J. Weiner regarding cosmological perturbation theory and for comments on our drafts. We also acknowledge the undergraduates responsible for probing cosmological perturbation theory in different contexts that eventually yielded the computational expertise to make this work possible: Arthur Conover, Maggie Murphree, Rachel Nguyen, Caroline Popiel, and Christian Solorio. J.T.G. is supported by the National Science Foundation Grant No. PHY-1719652. Numerical simulations here were performed on hardware purchased by the National Science Foundation, Kenyon College, and the Kenyon College department of physics, and we would like to acknowledge the infrastructure support provided by Kenyon College.

Appendix A Bardeen variables

In Sec. IV we relied on the ability to write the Bardeen potentials strictly in terms of BSSN variables. To calculate these we draw a correspondence between the BSSN metric, Eq. (16), and a general form of a metric with scalar perturbations Weinberg:2008zzc

[TABLE]

In terms of these variables, the Bardeen potentials are Baumann:2009ds

[TABLE]

We first set the purely spatial parts of the metric to be equal in both gauges, $g_{ij}^{\rm{BSSN}}=g_{ij}^{\rm{CPT}}$

[TABLE]

There are two relationships that we can get from this. First, by taking the trace of both sides, where repeated lower indices implies a sum, we get

[TABLE]

We can solve this for $\nabla^{2}E$,

[TABLE]

or in Fourier space variables (denoted by overtildes)

[TABLE]

where we have dropped the $3/2$ since it only contributes to the homogeneous mode which must be zero for a perturbation. We can also take the mixed second derivative of Eq. (46), $\partial_{i}\partial_{j}$, where again we sum over repeated downstairs indices,

[TABLE]

which can be rearranged

[TABLE]

and simplified by using Eq. (48),

[TABLE]

[TABLE]

This yields an expression for $\nabla^{2}\Psi$,

[TABLE]

or in Fourier space,

[TABLE]

This allows us to calculate $E$ strictly in terms of BSSN variables and $a$ (which will later be replaced by $\left<\rm{e}^{4\phi}\right>$),

[TABLE]

The last thing we need is $B$, which we get from $g_{0i}$,

[TABLE]

Introduce a sum by differentiating both sides,

[TABLE]

which in Fourier space becomes

[TABLE]

The scalar metric perturbations in terms of $a$ and BSSN variables can be summarized

[TABLE]

where repeated indices are summed over. We define the term common in both Bardeen potentials, Eqs. (44) and (45), as $\mathcal{B}\equiv a^{2}\left(\dot{E}-B/a\right)$. In BSSN variables we then have

[TABLE]

where we the Hubble rate with the average extrinsic curvature $H=-\left<K\right>/3$ as we did when linearizing the Hamiltonian constraint, Eq. (II.3). The Bardeen potentials are then computed strictly in terms of BSSN variables,

[TABLE]

The time derivative in Eq. (65) is very cumbersome, and we do not need to evolve the Bardeen potentials as they are only used to compare our BSSN results with CPT computations. Therefore it is sufficient to store $\mathcal{B}$ from the previous step and calculate $\dot{\mathcal{B}}$ with a Euler method. Incidentally, this means we cannot calculate $\dot{\mathcal{B}}$ until the second time step and also means we introduce a fixed numerical error of order $dx^{2}$ to $\Phi_{B}$.

Appendix B Code Verification

We have run a number of tests on GABERel to verify Szilagyi:1999nu that our gravitational evolution is accurate. Among these are the robust stability test and a Schwarzchild black hole test. In the former we find no exponential noise growth and the constraint violations remain bounded for long-term evolution. For the latter, we simulated a Schwarzschild geometry in trumpet coordinates Dennison:2014sma . In this setup, the conformal factor $\phi$, lapse $\alpha$ and shift $\beta^{i}$ take the form,

[TABLE]

which corresponds to a choice of $R_{0}=M$; N.B. for this test, we use units where $G=1$ and, hence, space and time are in units of $M$. We chose this particular slicing because all of the degrees of freedom of the extrinsic curvature, $K$ and $\tilde{A}_{ij}$, are nonzero throughout the simulation, even though the solution is static. The black hole is initially centered $dx/4$ from a central lattice point. This intentional asymmetry constitutes a more general and robust test of our code’s gravitational dynamics than would otherwise be the case if we impose some artificial symmetry on the system. We evolve the gauge variables, $\alpha$ and $\beta^{i}$, according to $1+\log$ and $\eta=0$ Gamma-driver conditions, respectively, with advective shift terms,

[TABLE]

While the rest of our simulation uses central finite differencing in an RK4 scheme, anything evolved with advective shift terms uses upwind derivative stencils. Gamma driver serves to approximate the more difficult to evolve Gamma-freezing condition, $\partial_{t}\bar{\Gamma}^{i}=0$. Figure 9 shows the stability of our code to this setup. Given our simulation’s periodic boundary conditions and the asymmetric initial position of the singularity, we do not necessarily expect that our solution should approach a perfectly steady state after long-term evolution. Nevertheless we observe no noticeable changes between the initial slice and $t=10\,M$. In this simulation, the box is taken to be $L=8\,M$ with $N=64^{3}$ points.

It’s also important for us to verify that we can trust the code for the simulations presented in Sec. IV. The constraints II.3 and 23 should remain small and bounded throughout the system’s evolution. However, the constraints are dimensionful, and so we must normalize them. Here we will show the Hamiltonian constraint and normalize it by the root sum of the squares of the terms in the constraint,

[TABLE]

We show the evolution of the Hamiltonian constraint in Figs. 10 and 11. It is not surprising that the constraint grows in the third phase of preheating as this is precisely when significant power is transferred to the smallest resolvable scales and numerical error grows. After this point, the numerical evolution is stable but no longer reliable.