A Bumpy Start to a Smooth Ride: Onset of Inflation amid Backreaction from Inhomogeneities

TL;DR

This study demonstrates that inflation in large-field models remains robust even with significant initial inhomogeneities, showing that nonlinear backreaction effects do not diminish the total volume of initial conditions leading to sufficient inflation.

Contribution

It introduces a nonperturbative Hartree approximation to analyze the impact of inhomogeneities on inflation onset, confirming robustness of inflation under complex initial conditions.

Findings

Inflation is robust for large-field models despite initial inhomogeneities.

The total volume of initial conditions producing sufficient inflation remains conserved with backreaction.

Some initial conditions that fail without inhomogeneities succeed when backreaction is included.

Abstract

We analyze the onset of inflation for a simple single-field model when the system begins with significant inhomogeneities on length-scales shorter than the initial Hubble radius. We incorporate certain nonlinear interactions among the coupled degrees of freedom by using the nonperturbative Hartree approximation. Consistent with recent, more computationally intensive numerical-relativity studies, we find inflation to be robust for large-field models, even when the system begins with significant structure on sub-Hubble scales. We consider the space of initial conditions $(\varphi (t_0), \dot{\varphi} (t_0))$, where $\varphi$ is the vacuum expectation value of the quantized field that drives inflation. Although some regions of $(\varphi (t_0), \dot{\varphi} (t_0))$ that would have yielded sufficient inflation in the absence of inhomogeneities fail to do so when backreaction from inhomogeneities is incorporated, an equal volume of such regions succeeds in producing sufficient inflation which did not do so in the absence of inhomogeneities. For large-field models, in other words, the total volume of the space of initial conditions $(\varphi (t_0), \dot{\varphi} (t_0))$ that yields sufficient inflation is {\it conserved} when we incorporate nonlinear backreaction from inhomogeneities, compared to the case in which inhomogeneities are neglected.

Figures

Equations

$$d{s}^2=-(1+2\Phi )d{t}^2+a^2(t)(1-2\Psi )h_{ij}(\mathbf{x})d{x}^{i}d{x}^j,$$

$$h_{ij}(\mathbf{x})d{x}^{i}d{x}^j=\frac{d{r}^2}{(1-K{r}^2)}+r^2\left(d{\theta }^2+{\sin }^2\theta d{\varphi }^2\right).$$

$$S=\int d^{4}x\sqrt{-g}\left[\frac{M_{\mathrm{pl}}^2}{2}R-\frac{1}{2}{g}^{\mu \nu }{\partial }_{\mu }\varphi {\partial }_{\nu }\varphi -V(\varphi )\right].$$

$$\square \varphi -V_{,\varphi }=0$$

$$G_{\mu \nu }\equiv R_{\mu \nu }-\frac{1}{2}{g}_{\mu \nu }R=\frac{1}{M_{\mathrm{pl}}^2}{T}_{\mu \nu },$$

$$T_{\mu \nu }={\partial }_{\mu }\varphi {\partial }_{\nu }\varphi -g_{\mu \nu }\left[\frac{1}{2}{g}^{\alpha \beta }{\partial }_{\alpha }\varphi {\partial }_{\beta }\varphi +V(\varphi )\right].$$

$$\square \varphi =\frac{1}{\sqrt{-g}}{\partial }_{\mu }\left[\sqrt{-g}{g}^{\mu \nu }{\partial }_{\nu }\varphi \right].$$

$$\varphi (x^{\mu })\to \hat{\varphi }(x^{\mu })=\phi (t)+\delta \hat{\varphi }(x^{\mu }),$$

$$⟨\hat{\varphi }(x^{\mu })⟩\equiv ⟨0|\hat{\varphi }(x^{\mu })|0⟩=\phi (t),⟨\delta \hat{\varphi }(x^{\mu })⟩=0.$$

$$\delta \hat{\varphi }(x^{\mu })=\int d\tilde{\mu }(k)\left[\delta {\varphi }_{k\ell m}(t){\hat{a}}_{k\ell m}{Z}_{k\ell m}(\mathbf{x})+H.c.\right],$$

$$\int d\tilde{\mu }(k)={\int }_0^{\infty }dk\sum _{\ell =0}^{\infty }\sum _{m=-\ell }^{\ell },$$

$${\nabla }^2{Z}_{k\ell m}\equiv \frac{1}{\sqrt{h}}{\partial }_i\left[\sqrt{h}{h}^{ij}{\partial }_j{Z}_{k\ell m}\right]=-k^2{Z}_{k\ell m}.$$

$$\begin{array}{rl}[{\hat{a}}_{k\ell m},{\hat{a}}_{k^{\prime }{\ell }^{\prime }{m}^{\prime }}^{†}]& =\delta (k-k^{\prime }){\delta }_{\ell {\ell }^{\prime }}{\delta }_{m{m}^{\prime }},\\ [{\hat{a}}_{k\ell m},{\hat{a}}_{k^{\prime }{\ell }^{\prime }{m}^{\prime }}]& =[{\hat{a}}_{k\ell m}^{†},{\hat{a}}_{k^{\prime }{\ell }^{\prime }{m}^{\prime }}^{†}]=0\end{array}$$

$${\hat{a}}_{k\ell m}|0⟩=0,⟨0|{\hat{a}}_{k\ell m}^{†}=0$$

$$\hat{\Psi }(x^{\mu })=\int d\tilde{\mu }(k)\left[{\Psi }_{k\ell m}(t){\hat{a}}_{k\ell m}{Z}_{k\ell m}(\mathbf{x})+H.c.\right].$$

$$\begin{array}{rl}(\delta \hat{\varphi }{)}^2& \to ⟨(\delta \hat{\varphi }{)}^2⟩,\\ (\delta \hat{\varphi }{)}^3& \to 3⟨(\delta \hat{\varphi }{)}^2⟩\delta \hat{\varphi },\\ (\delta \hat{\varphi }{)}^4& \to 3⟨(\delta \hat{\varphi }{)}^2{⟩}^2,\\ (\delta \hat{\varphi }{)}^5& \to 15⟨(\delta \hat{\varphi }{)}^2{⟩}^2\delta \hat{\varphi },\end{array}$$

$$\ddot{\phi }+3H\dot{\phi }+V^{(1)}(\phi )+\frac{1}{2}{V}^{(3)}(\phi )⟨(\delta \hat{\varphi }{)}^2⟩=0$$

$$\begin{array}{rl}& \delta {\ddot{\varphi }}_{k\ell m}+3H\delta {\dot{\varphi }}_{k\ell m}\\ & +\left[\frac{k^2}{a^2}+V^{(2)}(\phi )+\frac{1}{2}{V}^{(4)}(\phi )⟨(\delta \hat{\varphi }{)}^2⟩\right]\delta {\varphi }_{k\ell m}\\ & =2(\ddot{\phi }+3H\dot{\phi }){\Psi }_{k\ell m}+4\dot{\phi }{\dot{\Psi }}_{k\ell m},\end{array}$$

$$V^{(n)}(\varphi)\equiv\left(\frac{d^{n}V(\phi)}{d\phi^{n}}\right)\bigg{|}_{\phi=\varphi}.$$

$${\partial }_i\left(\dot{\hat{\Psi }}+H\hat{\Psi }\right)=\frac{1}{2M_{\mathrm{pl}}^2}\left(\dot{\phi }+\delta \dot{\hat{\varphi }}\right){\partial }_i\delta \hat{\varphi }.$$

$${\dot{\Psi }}_{k\ell m}+H{\Psi }_{k\ell m}=\frac{1}{2M_{\mathrm{pl}}^2}\dot{\phi }\delta {\varphi }_{k\ell m}.$$

$$\begin{array}{rl}\bar{\rho }(t)& \equiv \frac{1}{2}{\dot{\phi }}^2+V(\phi ),\\ \delta {\hat{\rho }}_{(1)}(x^{\mu })& \equiv \dot{\phi }\delta \dot{\hat{\varphi }}-{\dot{\phi }}^2\hat{\Psi }-⟨(\delta \dot{\hat{\varphi }}{)}^2⟩\hat{\Psi }+\frac{1}{a^2}\hat{\Psi }⟨(\nabla \delta \hat{\varphi }{)}^2⟩\\ & +V^{(1)}(\phi )\delta \hat{\varphi }+\frac{1}{2}{V}^{(3)}(\phi )⟨(\delta \hat{\varphi }{)}^2⟩\delta \hat{\varphi },\\ \delta {\rho }_{(2)}(t)& \equiv \frac{1}{2}⟨(\delta \dot{\hat{\varphi }}{)}^2⟩+\frac{1}{2a^2}⟨(\nabla \delta \hat{\varphi }{)}^2⟩\\ & +\frac{1}{2}{V}^{(2)}(\phi )⟨(\delta \hat{\varphi }{)}^2⟩+\frac{1}{8}{V}^{(4)}(\phi )⟨(\delta \hat{\varphi }{)}^2{⟩}^2,\end{array}$$

$$⟨(\nabla \delta \hat{\varphi }{)}^2⟩\equiv h^{ij}⟨{\partial }_i\delta \hat{\varphi }{\partial }_j\delta \hat{\varphi }⟩.$$

$$H^2=\frac{1}{3M_{\mathrm{pl}}^2}\left[\bar{\rho }+\delta {\rho }_{(2)}\right]-\frac{K}{a^2}.$$

$$\begin{split}&-6H\left(\dot{\Psi}_{k\ell m}+H\Psi_{k\ell m}\right)+\frac{2}{a^{2}}\left(3K-k^{2}\right)\Psi_{k\ell m}\\ &\quad=\frac{1}{M_{\rm pl}^{2}}\bigg{\{}\dot{\varphi}\,\delta\dot{\phi}_{k\ell m}-\dot{\varphi}^{2}\Psi_{k\ell m}-\langle(\delta\hat{\phi})^{2}\rangle\Psi_{k\ell m}\\ &\quad\quad\quad\quad\quad+\frac{1}{a^{2}}\Psi_{k\ell m}\langle(\nabla\delta\hat{\phi})^{2}\rangle+V^{(1)}(\varphi)\,\delta\phi_{k\ell m}\\ &\quad\quad\quad\quad\quad\quad\quad+\frac{1}{2}V^{(3)}(\varphi)\langle(\delta\hat{\phi})^{2}\rangle\,\delta\phi_{k\ell m}\bigg{\}}.\end{split}$$

$$\begin{array}{rl}\bar{p}(t)& \equiv \frac{1}{2}{\dot{\phi }}^2-V(\phi ),\\ \delta p_{(2)}(t)& \equiv \frac{1}{2}⟨(\delta \dot{\hat{\varphi }}{)}^2⟩-\frac{1}{6a^2}⟨(\nabla \delta \hat{\varphi }{)}^2⟩\\ & -\frac{1}{2}{V}^{(2)}(\phi )⟨(\delta \hat{\varphi }{)}^2⟩-\frac{1}{8}{V}^{(4)}(\phi )⟨(\delta \hat{\varphi }{)}^2{⟩}^2.\end{array}$$

$$\left(2\dot{H}+3H^2+\frac{K}{a^2}\right)=-\frac{1}{M_{\mathrm{pl}}^2}\left[\bar{p}+\delta p_{(2)}\right].$$

$$\begin{array}{rl}\dot{H}& =-\frac{1}{2M_{\mathrm{pl}}^2}\left[\bar{\rho }+\bar{p}+\delta {\rho }_{(2)}+\delta p_{(2)}\right]\\ & =-\frac{1}{2M_{\mathrm{pl}}^2}\left[{\dot{\phi }}^2+⟨(\delta \dot{\hat{\varphi }}{)}^2⟩+\frac{1}{3a^2}⟨(\nabla \delta \hat{\varphi }{)}^2⟩\right].\end{array}$$

$$\begin{array}{rl}& \left[\dot{H}+\frac{2}{3M_{\mathrm{pl}}^2{a}^2}⟨(\nabla \delta \hat{\varphi }{)}^2⟩+\frac{1}{a^2}\left(k^2-3K\right)\right]{\Psi }_{k\ell m}\\ & =\frac{1}{2M_{\mathrm{pl}}^2}\left[\ddot{\phi }\delta {\varphi }_{k\ell m}-\dot{\phi }\delta {\dot{\varphi }}_{k\ell m}\right].\end{array}$$

$$\hat{\mathcal{R}}=\hat{\Psi }+\frac{H}{\dot{\phi }}\delta \hat{\varphi }.$$

Full text

††preprint: MIT-CTP/5111

I Introduction

Several observable features of our Universe today, including its spatial flatness, large-scale homogeneity, and the specific pattern of temperature anisotropies in the cosmic microwave background radiation, are readily explained if our Universe underwent a brief phase of early-universe inflation. (For reviews, see Refs. Guth and Kaiser (2005); Mukhanov (2005); Bassett et al. (2006); Lyth and Liddle (2009); Martin et al. (2014); Guth et al. (2014); Linde (2015); Baumann and McAllister (2015); Martin (2019).) The question remains, however, whether the onset of inflation itself is fairly generic, or whether inflation requires fine-tuned initial conditions. (For reviews, see Refs. Goldwirth and Piran (1992); Brandenberger (2017).)

Many analyses have highlighted the effectiveness of an inflationary attractor: at least for single-field models, if one neglects spatial inhomogeneities, then one can show that large classes of models will flow into inflation, across broad ranges of initial conditions and couplings, even for regions of phase space in which the system is initially dominated by the field’s kinetic (rather than potential) energy. The attractor behavior is especially effective in large-field models with sufficiently flat regions of the potential, for which $\epsilon_{n}\equiv(M_{\rm pl}\,\partial_{\phi})^{n}\ln V(\phi)\ll 1$, where $M_{\rm pl}\equiv 1/\sqrt{8\pi G}=2.43\times 10^{18}$ GeV is the reduced Planck mass Salopek and Bond (1990); Kung and Brandenberger (1990); Muller et al. (1990); Liddle et al. (1994); Mukhanov (2005); Tzirakis and Kinney (2007); Remmen and Carroll (2013); Martin et al. (2014); Azhar and Kaiser (2018); Chowdhury et al. (2019).

These analyses neglected perturbations in the spacetime metric, however, thereby leaving open the question of whether the onset of inflation remains robust even in the presence of significant inhomogeneities. Pioneering efforts began as early as the mid-1980s to numerically simulate the onset of inflation amid such inhomogeneities Albrecht et al. (1985, 1987); Kurki-Suonio et al. (1987); Feldman and Brandenberger (1989); Brandenberger and Feldman (1989); Goldwirth and Piran (1989, 1990); Goldwirth (1991); Laguna et al. (1991); Kurki-Suonio et al. (1993); Goldwirth and Piran (1992); Brandenberger (2017). Yet these early studies did not yield an unambiguous answer to the question of whether the onset of inflation required an initial patch of spacetime to be smooth across a characteristic length-scale as large as (or perhaps considerably larger than) the Hubble radius. (See also Refs. Vachaspati and Trodden (2000); Ijjas et al. (2013); Berezhiani and Trodden (2015).)

Recently two groups have conducted studies of the onset of inflation that implement full $(3+1)$-dimensional numerical relativity; each found that large-field inflation (in single-field models) is strongly robust, even amid significant initial inhomogeneities East et al. (2016); Clough et al. (2017, 2018). (They also found that the onset of small-field inflation can be significantly more robust than previous semi-analytic treatments had suggested, although small-field models generically require more finely tuned initial conditions than large-field models do; see also Ref. Marsh et al. (2018).) In particular, each group found that large-field inflation would start even if the spacetime initially had no Hubble-sized smooth patches. Given the significant computational resources required to simulate the full Einstein field equations in $(3+1)$-dimensions, however, the latest studies have been restricted, to date, to a handful of specific forms of the potential $V(\phi)$ and to limited regions of parameter space East et al. (2016); Clough et al. (2017, 2018).

In this work, we develop a complementary numerical approach. Although by design our approach cannot capture the full range of gravitational effects that could (in principle) affect the onset of inflation, it incorporates certain nonlinear interactions while remaining significantly more efficient than the full-bore numerical relativity studies. It therefore facilitates the study of a wide range of models, across broad regions of phase space and parameter space, while tracking the evolution of as many as $n\sim{\cal O}(10^{2})$ coupled modes, spanning a wide range of length-scales. (It also offers a feasible means of generalizing to multifield models; see also Ref. Easther et al. (2014).) For the case we have investigated so far — a single-field model that yields large-field inflation — our numerical results are consistent with the fully relativistic simulations of Refs. East et al. (2016); Clough et al. (2017, 2018).

We work to linear order in metric perturbations, $\Psi(x^{\mu})$, but incorporate nonlinear interactions among the field fluctuations, $\delta\phi(x^{\mu})$, by adopting the (nonperturbative) Hartree approximation. We therefore capture some of the nonlinear backreaction effects that are absent in studies that work only to linear order in fluctuations. Whereas an ordinary quantum loop expansion corresponds to a power series in the coupling constant $\lambda$ and in $\hbar$, in the Hartree approximation we resum a particular class of Feynman diagrams (the so-called “cactus” or “daisy” diagrams) to all orders. Hence the Hartree approximation moves beyond perturbative series in $\lambda$, $\hbar$, or $\delta\phi$. (See, e.g., Refs. Jackiw (1974); Dolan and Jackiw (1974); Cornwall et al. (1974); Chang (1975); Boyanovsky et al. (1994).)

By working only to linear order in $\Psi(x^{\mu})$, our approach cannot capture extreme gravitational phenomena like the formation of trapped surfaces or the collapse of local regions to form black holes. But the fully relativistic simulations in Refs. East et al. (2016); Clough et al. (2017, 2018) confirm that the formation of such black holes does not interrupt the overall flow of the system into inflation: regions outside of the collapse continue to expand, and the density of any such collapsed regions quickly dilutes. Moreover, as discussed in Ref. Clough et al. (2017), no black hole that forms from pre-inflationary overdensities can grow so large as to encompass the entire Hubble sphere, so even the largest black holes that might form fail to disrupt the overall flow into inflation; see also Ref. Kleban and Senatore (2016). Therefore our simplified treatment can complement the more computationally intensive studies.

Another feature of our approach is that we treat the origin of inhomogeneities as ultimately quantum-mechanical. One of the most significant achievements of inflationary cosmology is to provide a first-principles description of the primordial inhomogeneities that seed large-scale structure as arising from quantum-mechanical fluctuations of matter fields during inflation Guth and Kaiser (2005); Mukhanov (2005); Bassett et al. (2006); Lyth and Liddle (2009); Martin et al. (2014); Guth et al. (2014); Linde (2015); Baumann and McAllister (2015); Martin (2019). We aim to analyze the onset of inflation in a comparable way; after all, if there were a priori reasons to expect some pattern of classical inhomogeneities at arbitrarily early times, we wouldn’t need a mechanism like inflation to account for large-scale structure. Given the high energies involved and the possibility of nonlinear interactions among the quantized matter fields at very early times, we consider scenarios in which the initial state of the quantum fluctuations departs significantly from the minimum-energy (Bunch-Davies) state. The larger quantum fluctuations seed significant spatial inhomogeneities and affect the dynamics of the system.

Our goal in this paper is to establish our formalism, introduce details of our numerical approach, and apply our system to a simple model; we defer detailed applications to a wider range of models, including both large-field and small-field models, single-field and multifield cases, to future work. In Section II we derive the coupled equations of motion for gravitational and matter degrees of freedom within the Hartree approximation. In Section III we discuss initial conditions for the field fluctuations and our UV regularization scheme. Section IV presents our main numerical results, indicating that across broad regions of phase space, a simple model like $V(\phi)=\lambda\phi^{4}/4$ will flow into the inflationary attractor even for initial inhomogeneities as large as $|\Psi(x^{\mu})|\lesssim 0.5$. In particular, we find a shift in the regions of phase space that support sufficient inflation (compared to the case in which we neglect all perturbations), but no reduction in the total volume of the space of initial conditions $(\varphi(t_{0}),\dot{\varphi}(t_{0}))$ that yields sufficient inflation, where $\varphi$ is the vacuum expectation value of the quantized field that drives inflation. (We consider $N\geq 65$ efolds of inflation to be “sufficient” for addressing the usual shortcomings of standard hot big bang cosmology Dodelson and Hui (2003); Liddle and Leach (2003); Mukhanov (2005); Bassett et al. (2006); Lyth and Liddle (2009); Martin (2019).) Concluding remarks follow in Section V. In Appendix A we present more details of our discrete spectrum for the fluctuations. Appendix B provides additional information about how we set initial conditions for the field fluctuations, while Appendix C discusses the initialization of the metric perturbations. In Appendix D we discuss various numerical convergence tests.

II Equations of Motion in the Hartree Approximation

We work in $(3+1)$ spacetime dimensions and use units in which $c=\hbar=1$. We consider scalar metric perturbations around a background Friedmann-Lemaître-Robertson-Walker (FLRW) line-element, and work in longitudinal gauge,

[TABLE]

where $\Phi(x^{\mu})$ and $\Psi(x^{\mu})$ are scalar functions. As usual, the background metric on (comoving) spatial sections may be written

[TABLE]

As we will see below, within the Hartree approximation the anisotropic pressure vanishes, so that $\Phi(x^{\mu})=\Psi(x^{\mu})$.

As is well known, there is no straightforward generalization of longitudinal gauge beyond linear order in $\Psi(x^{\mu})$. In particular, at second order in metric perturbations, first-order scalar perturbations source tensor perturbations, and so on. Gauge subtleties therefore affect any perturbative calculation that aims to move beyond ${\cal O}(\Psi)$ while incorporating only scalar degrees of freedom. (For a recent discussion, see Appendix B of Ref. Giblin et al. (2019).) Hence we restrict our analysis to first order in spatially varying quantities.

The restriction to linear order in metric perturbations, treated in longitudinal gauge, can be reasonably well motivated. As demonstrated in Ref. Giblin et al. (2019), corrections from a fully relativistic treatment compared to a linear treatment in longitudinal gauge typically scale as ${\cal O}(\Psi^{2})$, unlike the case for synchronous gauge, in which relativistic corrections to the linearized treatment can be as large as ${\cal O}(1)$. In fact, as found in Ref. Giblin et al. (2019), linearized scalar metric perturbations in longitudinal gauge tend to exaggerate gravitational effects on length-scales longer than the Hubble radius, $\ell>H^{-1}(t_{0})$, compared to a fully relativistic treatment. Given the prior focus on whether inflation can start amid inhomogeneities with typical length-scales $\ell\leq H^{-1}(t_{0})$ Kurki-Suonio et al. (1987); Laguna et al. (1991); Kurki-Suonio et al. (1993); Goldwirth and Piran (1989, 1990); Goldwirth (1991); Goldwirth and Piran (1992); Brandenberger (2017); Albrecht et al. (1985, 1987); Feldman and Brandenberger (1989); Brandenberger and Feldman (1989); Vachaspati and Trodden (2000); Ijjas et al. (2013); Berezhiani and Trodden (2015), we will be most interested in gravitational effects on sub-Hubble length-scales. In longitudinal gauge, meanwhile, large perturbations, $|\Psi(x^{\mu})|\geq 0.5$, can lead to coordinate singularities. Hence the regime of interest is $|\Psi(x^{\mu})|\lesssim{\cal O}(0.5)$.

We consider single-field models for which the action may be written

[TABLE]

Varying the action with respect to $\phi$ and $g_{\mu\nu}$ yields the coupled equations of motion,

[TABLE]

and

[TABLE]

with

[TABLE]

As usual, the covariant d’Alembertian operator is given by

[TABLE]

In addition to expanding the gravitational degrees of freedom to first order in $\Psi(x^{\mu})$, we also consider fluctuations in the field $\phi$. Upon quantizing the field, we have

[TABLE]

with

[TABLE]

We expand the field fluctuations as

[TABLE]

where “$H.c.$” denotes Hermitian conjugate, the measure $d\tilde{\mu}(k)$ is given by

[TABLE]

and $Z_{k\ell m}({\bf x})$ is an eigenfunction of the comoving spatial Laplacian operator,

[TABLE]

The creation and annihilation operators obey the usual commutation relations

[TABLE]

and satisfy

[TABLE]

for all $k$, $\ell$, and $m$. As we will see, upon expanding the field $\phi$ as in Eq. (8) and working to linear order in $\Psi$, Eqs. (4) and (5) will couple $\Psi(x^{\mu})$ to terms linear in $\delta\hat{\phi}(x^{\mu})$. To remain consistent, we therefore quantize $\Psi(x^{\mu})\rightarrow\hat{\Psi}(x^{\mu})$. Because we are working with single-field models, we only require a single set of operators $\hat{a}_{k\ell m}$ and $\hat{a}^{\dagger}_{k\ell m}$. We expand $\hat{\Psi}$ as

[TABLE]

Given our interest in inflationary models, we focus on weakly coupled systems. For such systems, one expects the $S$ matrix to be dominated by forward-scattering processes. Hard-scattering events, which impart large transverse momenta to the scattered particles, should be relatively rare, such that most of the time the particles will emerge from an interaction with (nearly) the same momentum as they had prior to the interaction. The dominant processes, in other words, involve particles propagating along a given trajectory, but moving with a modified, effective mass induced by the (self-)interactions. In such cases, we may approximate the effects of nonlinear self-interactions by calculating a dressed propagator: we select the dominant subset of Feynman diagrams at each order of perturbation theory that contribute to the effective mass, and sum all members of that subclass to all orders, while neglecting the other terms that would appear in a full expansion of the $S$ matrix. (See, e.g., Section 4.7 of Ref. Mattuck (1976).) Within the Hartree approximation, we construct the dressed propagator by performing an infinite resummation of the so-called “cactus” or “daisy” diagrams Jackiw (1974); Dolan and Jackiw (1974); Cornwall et al. (1974); Chang (1975); Boyanovsky et al. (1994).

The Hartree approximation becomes exact for an $O({\cal N})$-symmetric model of ${\cal N}$ interacting scalar fields in the limit ${\cal N}\rightarrow\infty$. In that case, the amplitudes associated with the set of Feynman diagrams picked out for resummation remain parametrically larger than all other terms in the expansion of the $S$ matrix, order by order Jackiw (1974); Dolan and Jackiw (1974); Cornwall et al. (1974). Yet even in the case ${\cal N}=1$, the diagrams that contribute to the Hartree approximation dominate among the contributions to the effective mass at a given order in perturbation theory, and hence Hartree is often used to incorporate nonperturbative effects even for ${\cal N}=1$, as we will do here. (See, e.g., Refs. Chang (1975); Boyanovsky et al. (1994); Khlebnikov and Tkachev (1997); Bassett and Tsujikawa (2001).)

Operationally, this means that within the equations of motion, all terms that are higher order in the field fluctuations, of the form $(\delta\hat{\phi})^{n}$ for $n\geq 2$, are replaced by Cornwall et al. (1974); Chang (1975); Boyanovsky et al. (1994)

[TABLE]

where the particular coefficients on the right-hand side arise from the combinatorics of the various Wick contractions. (One may continue in a similar manner for $(\delta\hat{\phi})^{n}$ with $n>5$, though for the particular model of interest to us here, the higher-order terms will not be relevant.) Because the background spacetime around which we are perturbing is homogeneous and isotropic, the (dressed) two-point function $\langle(\delta\hat{\phi})^{2}\rangle$ is spatially homogeneous.

The Hartree approximation is nonperturbative, so we make no assumption about the relative magnitude of $\langle(\delta\hat{\phi})^{2}\rangle$ compared to $\varphi^{2}$. In particular, when expanding Eqs. (4) and (5), we retain terms of the form $\langle(\delta\hat{\phi})^{2}\rangle\hat{\Psi}$ as well as $\varphi^{2}\hat{\Psi}$. On the other hand, because we are working only to linear order in $\hat{\Psi}$, we do not include any terms of the form $\langle\hat{\Psi}\,\delta\hat{\phi}\rangle$. Within the Hartree approximation, such terms would arise from summing diagrams involving virtual $\hat{\Psi}$ quanta; yet even the bare propagator for such lines, $\Delta_{F}(x-y)$, is ${\cal O}(\hat{\Psi}^{2})$, and hence remains beyond our approximation. To linear order in metric perturbations, in other words, the perturbations $\hat{\Psi}(x^{\mu})$ do not contribute to the Hartree corrections, though, as we will see, the evolution of $\hat{\Psi}(x^{\mu})$ is sensitive to the nonlinear evolution of $\delta\hat{\phi}(x^{\mu})$.

Lastly, we note that the Hartree approximation is spherically symmetric in $k$-space; it does not include any direct mode-mode coupling, which would arise from convolutions of the sort $\int d\tilde{\mu}(k^{\prime})\,d\tilde{\mu}(q)\,\delta\phi_{k-k^{\prime}-q}\,\delta\phi_{k-k^{\prime}}\,\delta\phi_{k}$. (This is consistent with neglecting scattering events that would impart large transverse momenta.) By neglecting such rescattering effects, the Hartree approximation is less efficient at transferring power between modes of different length-scales than a fully nonlinear analysis would be. (See, e.g., the discussion in Ref. Bassett et al. (1999).) On the other hand, any such terms would require moving beyond linear order in the metric perturbations $\hat{\Psi}(x^{\mu})$ — since they would be higher order in spatially varying quantities — and hence the Hartree approximation is especially well-suited for any study that is restricted to linear order in gravitational degrees of freedom.

Our procedure is to expand Eqs. (4) and (5) to linear order in $\hat{\Psi}$ and to arbitrarily higher order in $(\delta\hat{\phi})^{n}$; implement the Hartree approximation to replace higher-order terms $(\delta\hat{\phi})^{n}$ as in Eq. (16); and discard any remaining terms that are beyond linear order in spatially varying quantities. Expanding Eq. (4), we find a set of terms that are spatially homogeneous, and a set of terms that are linear in spatially varying quantities. Requiring each set to vanish yields the coupled equations of motion:

[TABLE]

and

[TABLE]

where overdots denote derivatives with respect to cosmic time $t$, $H(t)\equiv\dot{a}/a$, and we use the notation

[TABLE]

(Because we have in mind application to a model with $V(\phi)=\lambda\phi^{4}/4$ in this paper, no terms with $V^{(n)}(\varphi)$ appear for $n\geq 5$.) The term in Eq. (18) proportional to $V^{(4)}(\varphi)\,\langle(\delta\hat{\phi})^{2}\rangle$, which contributes to the effective mass for the modes $\delta\phi_{k\ell m}$, would not appear if we had performed a perturbative loop expansion. It appears in Eq. (18) because the Hartree approximation yields a self-consistent gap equation for the dressed propagator Chang (1975); Boyanovsky et al. (1994).

From the $0i$ component of Eq. (5), we find

[TABLE]

The Hartree approximation stipulates that any terms quadratic in the field fluctuations $\delta\hat{\phi}$ should be replaced by the corresponding vacuum expectation value. In this case, the relevant term would be $\langle\delta\dot{\hat{\phi}}\,\partial_{i}\delta\hat{\phi}\rangle$, which will vanish: the resulting integrand is an odd function of $k_{i}$, integrated over symmetric limits. Hence this term vanishes within the Hartree approximation, and we find

[TABLE]

In a similar way, the anisotropic pressure that could arise from $T_{ij}=\partial_{i}\delta\hat{\phi}\,\partial_{j}\delta\hat{\phi}$ (for $i\neq j$) vanishes within the Hartree approximation, which is why the metric perturbations $\Phi(x^{\mu})$ and $\Psi(x^{\mu})$ in Eq. (1) remain equal to each other.

We next expand $T^{0}_{\>\>0}=-\rho$ to find the various contributions to the energy density. We find three distinct contributions: $\bar{\rho}(t)$, which depends only on the (homogeneous) vacuum expectation value of the field, $\varphi(t)$; $\delta\hat{\rho}_{(1)}(x^{\mu})$, which includes all terms that are linear in spatially varying quantities; and $\delta\rho_{(2)}(t)$, which includes all spatially homogeneous terms that arise from the fluctuations:

[TABLE]

where

[TABLE]

We see that $\langle\delta\hat{\rho}_{(1)}\rangle=0$. Upon expanding $G^{\mu}_{\>\>\nu}$ to first order in $\Psi$ and equating the terms from the $00$ component of Eq. (5) that are spatially homogeneous, we find the effective Friedmann equation

[TABLE]

Equating the terms in the $00$ component that are linear in spatially varying quantities, we have

[TABLE]

From the $ij$ component of Eq. (5), we are interested in extracting the spatially homogeneous contributions to the pressure, since these are relevant to the evolution of the background spacetime; in particular, we will use these terms to solve for $\dot{H}$. As usual we may parameterize the pressure as $T^{i}_{\>\>i}=3p$, and hence, adopting notation as above, we find

[TABLE]

(We will not need an explicit expression for $\delta\hat{p}_{(1)}$.) Equating the spatially homogeneous terms in the $ij$ component of Eq. (5) yields

[TABLE]

Combining Eqs. (24) and (27), we find

[TABLE]

From Eqs. (24) and (28), we see that the evolution of the background spacetime depends on the homogeneous field $\varphi$ as well as on contributions from two-point functions of the fluctuations. A welcome feature of the Hartree approximation is that the Hartree-corrected energy-momentum tensor obeys the covariant conservation relation $\langle T^{\mu\nu}\rangle_{;\nu}=0$, ensuring that these evolution equations remain mutually consistent with the equations of motion in Eqs. (17) and (18).

Finally, we may combine our expressions from the $00$ and $0i$ components of Eq. (5) — in particular, Eqs. (21) and (25) — to derive a constraint equation relating the modes $\Psi_{k\ell m}$ to $\delta\phi_{k\ell m}$. Upon making algebraic substitutions from Eqs. (17) and (28), we find

[TABLE]

In our numerical simulations, we use Eq. (29) only to set initial conditions for the modes $\Psi_{k\ell m}(t_{0})$, based on the initial conditions for $\varphi(t_{0})$, $\dot{\varphi}(t_{0})$, $H(t_{0})$, $\delta\phi_{k\ell m}(t_{0})$, and $\delta\dot{\phi}_{k\ell m}(t_{0})$; we then evolve the metric perturbations over time using Eq. (21). Although the source term in Eq. (21) is linear in $\delta\phi_{k\ell m}$, the evolution of $\varphi(t)$, $H(t)$, and $\delta\phi_{k\ell m}(t)$ each incorporates nonlinear backreaction effects from the Hartree corrections.

When working to linear order in $\hat{\Psi}$ and $\delta\hat{\phi}$, it is common to study the evolution of the gauge-invariant comoving curvature perturbation, $\hat{\cal R}(x^{\mu})$, which (for single-field models) takes the form Bassett et al. (2006); Lyth and Liddle (2009)

[TABLE]

(As is well-known, $\hat{\cal R}$ is proportional to the gauge-invariant Mukhanov-Sasaki variable, and is equivalent to the curvature perturbation on uniform-density hypersurfaces, $\hat{\zeta}$, in the limit $k\ll aH$ Bassett et al. (2006); Lyth and Liddle (2009).) Although $\hat{\cal R}$ only remains gauge-invariant for linear gauge transformations, it remains a useful quantity for considering the evolution of perturbations even when we incorporate the nonlinear, nonperturbative Hartree corrections, as we will see in Section IV.

III Setting Parameters

In our numerical simulations, we track the evolution of coupled modes within a finite (comoving) spatial volume; this restriction, in turn, means that for any Gaussian curvature $K$ of the background spatial sections, we have a discrete spectrum of allowable wavenumbers. Then the integral over $dk$ in the measure $d\tilde{\mu}(k)$ defined in Eq. (11) is replaced by a discrete sum, indexed by an integer $n\geq 1$:

[TABLE]

and similarly for $\hat{\Psi}(x^{\mu})$. In place of Eq. (13), the creation and annihilation operators now obey $[\hat{a}_{n\ell m},\hat{a}^{\dagger}_{n^{\prime}\ell^{\prime}m^{\prime}}]=\delta_{nn^{\prime}}\delta_{\ell\ell^{\prime}}\delta_{mm^{\prime}}$, with $\hat{a}_{n\ell m}|0\rangle=\langle 0|\hat{a}^{\dagger}_{n\ell m}=0$ for all $(n\ell m)$.

For the remainder of this paper we consider $K=0$ and evolve the modes within a sphere of comoving radius $R$. (We defer the case of nonzero $K$ to future work.) As described in Appendix A, for $K=0$ the normalized spatial eigenfunctions $Z_{n\ell m}(r,\theta,\phi)$ take the form

[TABLE]

Here $N_{n\ell m}$ is a normalization constant, $j_{\ell}(x)$ is the spherical Bessel function, and $Y_{\ell m}(\theta,\phi)$ is the usual spherical harmonic. We choose to use Dirichlet boundary conditions at $r=R$, which fixes $Z_{n\ell m}(R,\theta,\phi)=0$ for all $(\theta,\phi)$, which in turn restricts the allowable wavenumbers $k_{n\ell}$ to a discrete spectrum, namely

[TABLE]

where $x_{n\ell}$ is the $n$th zero of the Bessel function $j_{\ell}(x)$, that is, $j_{\ell}(x_{n\ell})=0$ for $n\geq 1$. (For $\ell=0$, the $k_{n0}$ take the simple form, $k_{n0}=n\pi/R$.)

Within the Hartree approximation, the evolution of $\varphi(t)$ and $H(t)$, as well as the modes $\Psi_{n\ell m}(t)$, depends on the initial conditions for the modes $\delta\phi_{n\ell m}(t_{0})$ and $\delta\dot{\phi}_{n\ell m}(t_{0})$. Because the Hartree approximation replaces higher-order interaction terms in the equation of motion for the fluctuations $\delta\hat{\phi}$ by an effective mass, we may follow many of the usual steps for quantizing a free scalar field in FLRW spacetime, and use this quantization procedure to parameterize initial conditions for $\delta\phi_{n\ell m}(t_{0})$ and $\delta\dot{\phi}_{n\ell m}(t_{0})$.

In the regime of interest, the field fluctuations are nearly massless around $t_{0}$. From Eq. (18), the effective mass is given by

[TABLE]

which is suppressed by the small coupling constant $\lambda$; hence we have $m_{\rm eff}^{2}(t_{0})\ll H^{2}(t_{0})$. Moreover, when we evolve the modes within a sphere of comoving radius $R$, we introduce an infrared cut-off given by $k_{\rm min}=\pi/R$, with $R\sim 1/[a(t_{0})\,H(t_{0})]$. Even for the longest-wavelength modes in our simulation, we therefore have $k_{n\ell}^{2}/a^{2}(t_{0})\gtrsim m_{\rm eff}^{2}(t_{0})$.

We do not assume that the system has attained the mininum energy state at the initial time $t_{0}$, and hence we consider initial conditions for $\delta\phi_{n\ell m}(t_{0})$ and $\delta\dot{\phi}_{n\ell m}(t_{0})$ that depart from the usual Bunch-Davies vacuum state. As described in Appendix B, we parameterize

[TABLE]

where $\alpha_{n\ell m}$, $\gamma_{n\ell m}$, and $\delta_{n\ell m}$ are each real-valued dimensionless constants and $\bar{H}_{0}$ is given by

[TABLE]

From Eq. (22) we note that $\bar{\rho}(t_{0})$ is the energy density associated with the vacuum expectation value at the initial time, $\varphi(t_{0})$. For the quantum fluctuations, the equal-time commutation relation at $t_{0}$ requires

[TABLE]

for all $(n\ell m)$. The Bunch-Davies initial state corresponds to $\gamma_{n\ell m}=1$ and $\delta_{n\ell m}=0$ for all $(n\ell m)$. To depart from the Bunch-Davies initial state, we treat the coefficients $\gamma_{n\ell m}$ and $\delta_{n\ell m}$ as random variables for each mode, drawn from flat distributions within specific ranges. Once $\gamma_{n\ell m}$ is drawn for a given mode, $\alpha_{n\ell m}$ for that mode is fixed from Eq. (37).

In addition to selecting initial conditions for individual modes $\delta\phi_{n\ell m}(t_{0})$ and $\delta\dot{\phi}_{n\ell m}(t_{0})$, we also need to evaluate the various two-point functions that appear in the evolution equations for $\varphi(t)$, $\delta\phi_{n\ell m}(t)$, $H(t)$, and $\dot{H}(t)$; only then can we set initial conditions for the metric perturbation modes $\Psi_{n\ell m}(t_{0})$ and evolve the coupled system forward in time. As discussed in Appendix A, when we evaluate the two-point functions within a sphere of comoving radius $R$, only modes with $\ell=0$ contribute to $\langle(\delta\hat{\phi})^{2}\rangle$ and $\langle(\delta\dot{\hat{\phi}})^{2}\rangle$, whereas only modes with $\ell=1$ contribute to $\langle(\nabla\delta\hat{\phi})^{2}\rangle$, and we find

[TABLE]

Although we have considered the $K=0$ case here, the pattern of which $\ell$ modes contribute to which two-point functions holds for arbitrary $K$, since the $(\theta,\phi)$ portion of the background metric in Eq. (2) does not depend on $K$.

Our expressions for the various two-point functions in Eq. (38) diverge in the limit $n\rightarrow\infty$; this is just the usual ultraviolet divergence that appears in the continuum limit for $k\rightarrow\infty$. Hence we must regularize all sums that appear in the various two-point functions. Since we expand the quantum fluctuations $\delta\hat{\phi}(x^{\mu})$ as sums over modes of comoving wavenumber $k_{n\ell}$, we adopt a simple Gaussian regulator with a comoving UV scale $\kappa$:

[TABLE]

We insert $F_{n\ell}(\kappa,R)$ within the sums when evaluating all two-point functions. For example,

[TABLE]

We use a single UV regulator scale $\kappa$, independent of $\ell$. Once $\kappa$ is fixed, the sums over $n$ may be truncated at some finite number of modes, $n_{\rm max}$.

The regularized two-point functions depend on the UV regulator scale $\kappa$. In general, the two-point function for quantum fluctuations of a nearly massless scalar field in an FLRW background scales as $\langle(\delta\hat{\phi})^{2}\rangle\propto H^{2}$. We therefore parameterize $\kappa=b\bar{H}_{0}$, with $b$ a real, dimensionless constant. To confirm the scaling of the regularized two-point function with $H$, we use the fact that the random coefficients $\alpha_{n\ell m}$ for the mode functions $\delta\phi_{n\ell m}(t_{0})$ vary independently of $n$, so we may take the term $\alpha_{n00}^{2}$ out of the sum in Eq. (40) at $t_{0}$ and replace it by an average value. As detailed in Appendix B, this yields

[TABLE]

(Although only modes with $\ell=m=0$ contribute to the sum in Eq. (40), we draw the random variables $\gamma_{n\ell m}=1/\alpha_{n\ell m}$ from the same distribution for all $(n\ell m)$, and hence the average value $(\alpha^{2})_{\rm avg}$ holds for any $n$, $\ell$, and $m$.) Following similar steps (see Appendix B) we find

[TABLE]

with

[TABLE]

The quantity ${\cal C}$ measures how strongly (on average) the amplitude of the initial state of the quantum fluctuations deviates from the minimum-energy (Bunch-Davies) state. Since $\gamma_{n\ell m}=\alpha_{n\ell m}=1$ and $\delta_{n\ell m}=0$ for all $(n\ell m)$ in the Bunch-Davies state, ${\cal C}_{\rm BD}=2$. Throughout our analysis, we consider quantum fluctuations whose average initial energy density exceeds the Bunch-Davies value by an order of magnitude, with ${\cal C}\simeq 20$, by drawing the random coefficients for each mode from flat distributions within the ranges

[TABLE]

Once $\gamma_{n\ell m}$ is drawn for a given mode, $\alpha_{n\ell m}$ for that mode is fixed from Eq. (37). The ranges in Eq. (44) yield $(\alpha^{2})_{\rm avg}=11.11$, $(\gamma^{2})_{\rm avg}=0.37$, and $(\delta^{2})_{\rm avg}=8.33$.

The actual value of the initial Hubble scale $H_{0}\equiv H(t_{0})$ will be greater than $\bar{H}_{0}$, because $H$ includes contributions from both $\bar{\rho}$ and $\delta\rho_{(2)}$, as in Eq. (24). We define $f\equiv(H_{0}/\bar{H}_{0})>1$; using Eqs. (36) and (42) we find that on average

[TABLE]

We aim to study initial conditions for the system such that $H_{0}\gg H_{\rm infl}$, where $H_{\rm infl}$ is the Hubble scale of the slow-roll inflationary attractor. Yet we also need to keep $H_{0}$ low enough (compared to $M_{\rm pl}$) so that we can begin the simulations with significant power in fluctuations on sub-Hubble scales, $H_{0}<k/a(t_{0})<M_{\rm pl}$, while avoiding trans-Planckian ambiguities. Hence we focus on initial conditions such that $H_{0}\sim 0.1\,M_{\rm pl}$. For perturbations that depart from the Bunch-Davies initial state, with $2\leq{\cal C}\leq 20$, we find $\delta\rho_{(2)}(t_{0})\lesssim\bar{\rho}(t_{0})$ for $\bar{H}_{0}=0.1\,M_{\rm pl}$ if we fix $b=5$. Eq. (45) then yields $f_{\rm avg}\leq 1.4$, corresponding to $H_{0}\lesssim 0.14\,M_{\rm pl}$.

Because we are interested in effects of initial inhomogeneities on length-scales shorter than the initial Hubble radius, we fix the comoving radius $R=1.5\pi\bar{H}_{0}^{-1}>1.5\pi H_{0}^{-1}$. (We set $a(t_{0})=1$.) Then the longest comoving wavelength in the spectrum is $\lambda_{\rm max}=2R>3\pi H_{0}^{-1}$, corresponding to $k_{\rm min}=\pi/R$. This ensures that our longest wavelength modes begin on a superhorizon scale, but that most of our modes are initially subhorizon. Given the form of $F_{n\ell}(\kappa,R)$ in Eq. (39), meanwhile, we find strong suppression of the contribution to the various two-point functions by modes with comoving wavenumber $k_{n\ell}>k_{UV}=3\kappa$, or $k_{UV}=3b\bar{H}_{0}$. Setting $b=5$, we find $k_{UV}\sim{\cal O}(10H_{0})\sim M_{\rm pl}$. To ensure strong numerical convergence we fix $k_{\rm max}=4\kappa$, which corresponds to $n_{\rm max}=30$. Our simulations then involve $4n_{\rm max}=120$ coupled modes: $n_{\rm max}$ each for $\ell=0$ and for $\ell=1$ with $m=-1,0,1$.

We implement our UV regularization via Eq. (39), but do not pursue formal renormalization. For one thing, we are interested in scenarios in which the system begins at high energies $H_{0}\sim 0.1\,M_{\rm pl}$, and we aim to track effects of excited modes with wavenumbers up to $k_{UV}\sim M_{\rm pl}$; hence there are no arbitrarily large hierarchies between the physical energy scales of interest and the natural cut-off scale. More important, by studying initial states for the fluctuations that depart from the usual Bunch-Davies state, as in Eq. (35), any renormalization scheme would depend both on the coupling constants in the Lagrangian and on the particular selection of initial parameters $\alpha_{n\ell m},\gamma_{n\ell m}$, and $\delta_{n\ell m}$ for the modes. Hence any renormalization scheme would change, run by run, with the random draws for these parameters. (See, e.g., Refs. Collins and Holman (2006); Collins et al. (2014).) We therefore defer questions about formal renormalization to future work.

IV Numerical Results

In this section we first consider typical behavior of the coupled system for a particular set of initial conditions, before turning to a more systematic investigation across initial conditions. We study $V(\phi)=\lambda\phi^{4}/4$ with $\lambda=10^{-10}$, and later compare results for $\lambda=10^{-12}$.

When one ignores field fluctuations and metric perturbations, this model yields sufficient inflation, with $N_{\rm infl}\geq 65$ efolds, for $\varphi_{0}\equiv\varphi(t_{0})\sim 15-30\,M_{\rm pl}$, depending on the initial value of $\dot{\varphi}_{0}\equiv\dot{\varphi}(t_{0})$. For $\lambda=10^{-10}$, this corresponds to a value of the Hubble parameter once the system has reached the slow-roll inflationary attractor (with $\dot{\varphi}^{2}\ll V(\varphi)$) of $H_{\rm infl}\sim 10^{-3}\,M_{\rm pl}$; for $\lambda=10^{-12}$, we have $H_{\rm infl}\sim 10^{-4}\,M_{\rm pl}$. We are therefore interested in the behavior of this system in the vicinity of $\varphi_{0}\sim 15-30\,M_{\rm pl}$ when the nonlinear effects of fluctuations are incorporated into the dynamics.

For each value of $\lambda$, we varied $12\,M_{\rm pl}\leq\varphi_{0}\leq 30\,M_{\rm pl}$ in steps of $\Delta\varphi_{0}=0.25\,M_{\rm pl}$, and $-0.25\,M_{\rm pl}^{2}\leq\dot{\varphi}_{0}\leq 0.25\,M_{\rm pl}^{2}$ in steps of $\Delta\dot{\varphi}_{0}=0.01\,M_{\rm pl}^{2}$, for a $73\times 51$ search grid. For each grid point in $(\varphi_{0},\dot{\varphi}_{0})$, we computed $34$ simulations: one with the Hartree corrections turned off (to neglect coupled fluctuations), one with the quantum fluctuations $(\delta\phi_{n\ell m}(t_{0}),\delta\dot{\phi}_{n\ell m}(t_{0}))$ initialized in the Bunch-Davies initial state with ${\cal C}_{\rm BD}=2$, and $32$ in which we initialized the quantum fluctuations with ${\cal C}\simeq 20$ by drawing random coefficients $\gamma_{n\ell m}$ and $\delta_{n\ell m}$ for each mode from the ranges in Eq. (44). This yielded a total of roughly 250,000 individual simulations. The simulations were performed on the Amazon Web Services EC2 service on a 16-core instance, and took a little under two days to complete. Our code was implemented in Python.

Our simulations are initialized in a number of steps. Starting from the values $(\varphi_{0},\dot{\varphi}_{0})$, we construct $\bar{H}_{0}$ as in Eq. (36), from which $R$, $\kappa$, and $k_{\rm max}$ are computed. The spectrum of allowable wavenumbers $k_{n\ell}(R)$ with $\ell=0$ and $\ell=1$ is then constructed. Next, we construct the Bunch-Davies initial conditions for each mode $\delta\phi_{n\ell m}(t_{0})$ and $\delta\dot{\phi}_{n\ell m}(t_{0})$. For perturbed initial data, we draw random values for $\gamma_{k\ell m}$ and $\delta_{k\ell m}$ for each mode, and construct the initial mode data appropriately from Eq. (35). We then compute the appropriate two-point functions, initialize the modes $\Psi_{n\ell m}(t_{0})$ from Eq. (29), and evaluate $\delta\rho_{(2)}(t_{0})$ to construct the actual initial Hubble factor $H_{0}$.

Our simulations evolve the quantities $a(t)$, $\varphi(t)$, and $\dot{\varphi}(t)$. For each wavenumber $k_{n\ell}$ with $\ell=0$ and $\ell=1$, we also evolve two modes: one initialized with $(\delta\phi_{n\ell 0}(t_{0}),\delta\dot{\phi}_{n\ell 0}(t_{0}))=(1,0)$ and a second with $(0,1)$ (along with a corresponding initialization for $\Psi_{n\ell 0}(t_{0})$ for each mode). Given the time evolution of these modes and the initial conditions $\delta\phi_{n\ell m}(t_{0})$, $\delta\dot{\phi}_{n\ell m}(t_{0})$, and $\Psi_{k\ell m}(t_{0})$, every mode can be reconstructed by exploiting the linearity of the equations of motion, since the nonlinear two-point functions are effectively functions of time that are independent of $n$, $\ell$ and $m$. Note that the modes that are evolved remain real. Doing this split makes little difference for the $\ell=0$ modes, but reduces the computational cost threefold for the $\ell=1$ modes. We evolve the system forward in time using a variable time-step RK45 algorithm. We declare that slow-roll begins once $\epsilon<0.1$, and we terminate evolution thereafter at $\epsilon\geq 1$, where

[TABLE]

IV.1 Evolution of Perturbations

In this subsection and the following, we present results from a representative sample of random initializations. We set $\lambda=10^{-10}$ and consider the case $\varphi_{0}=25\,M_{\rm pl}$ and $\dot{\varphi}_{0}=-0.25\,M_{\rm pl}^{2}$, for which the initial kinetic energy in the field $\varphi$ exceeds the potential energy by a factor of 3200. We initialize the fluctuations $\delta\hat{\phi}(t_{0},{\bf x})$ by parameterizing the mode functions as in Eq. (35) and drawing the random initial coefficients $\gamma_{n\ell m}$ and $\delta_{n\ell m}$ for each mode from the distributions in Eq. (44). Then ${\cal C}\simeq 20$ and the energy density in fluctuations $\delta\rho_{(2)}(t_{0})$ begins about ten times greater than in the Bunch-Davies initial state. With these parameters, we find $\delta\rho_{(2)}(t_{0})\simeq\bar{\rho}(t_{0})$ and $H_{0}=0.14\,M_{\rm pl}$, two orders of magnitude greater than $H_{\rm infl}$.

For these initial conditions, the system begins with significant inhomogeneities on length-scales well within the initial Hubble radius. Fig. 1 shows $\Psi(t_{0},{\bf x})$, constructed from modes $\Psi_{n\ell m}(t_{0})$ whose amplitudes are set by Eq. (29). (Further details of how we construct $\Psi(t_{0},{\bf x})$ are given in Appendix C.) The blue cylinder in Fig. 1 has a radius equal to $r_{H}(t_{0})=\pi/H(t_{0})$, such that modes with $k\geq a(t_{0})H(t_{0})$ have wavelengths that fit within the diameter $2r_{H}(t_{0})$. For this choice of initial conditions, the metric perturbations $\Psi(t_{0},{\bf x})$ begin with substantial structure on sub-Hubble length-scales, with spatial inhomogeneities as large as $|\Psi(t_{0},{\bf x})|\simeq 0.4$.

Figures 2 and 3 show $|{\rm Re}(\delta\phi_{n\ell m})|$ and $|{\rm Re}(\Psi_{n\ell m})|$ versus $N\equiv\int Hdt=\ln a$ for early times, for the $\ell=0$ and $\ell=1$ modes with $k_{n\ell}\geq a_{0}H_{0}$. As expected, the modes oscillate with decaying amplitude while inside the Hubble radius, and their amplitudes freeze after Hubble crossing. At later times, after the physical wavelengths of the modes have redshifted to be exponentially larger than the Hubble radius, the amplitudes show a modest secular growth, rising as $|\delta\phi_{n\ell m}|\propto\sqrt{\epsilon}$ and $|\Psi_{n\ell m}|\propto\epsilon$, where $\epsilon$ is given in Eq. (46). This modest late-time growth matches the well-known behavior of perturbations deep in the infrared during the slow-roll regime, as treated in linear perturbation theory. (See, e.g., Section 8.2 of Ref. Mukhanov (2005), as well as Ref. Seery (2010).) Nonetheless, the curvature perturbation, $\hat{\cal R}$, defined in Eq. (30), remains conserved once modes cross outside the Hubble radius. In Fig. 4 we plot the dimensionless power spectrum for the curvature perturbation,

[TABLE]

for $\ell=0$ and $\ell=1$ modes that begin with $k_{n\ell}\geq a_{0}H_{0}$. Consistent with the analytic results in Ref. Senatore and Zaldarriaga (2013), we find that the curvature perturbation remains conserved on super-Hubble length-scales, even when we incorporate nonlinear self-interactions. As shown in Fig. 4, substantial structure on sub-Hubble length-scales at early times damps out before modes cross outside the Hubble radius, producing a smooth patch on horizon scales, and remains exponentially suppressed for the duration of the simulation.

The energy density in fluctuations $\delta\rho_{(2)}$ begins with $\delta\rho_{(2)}(t_{0})\simeq\bar{\rho}(t_{0})$ and then begins to decay, as shown in Fig. 5. Because of the weak coupling $\lambda$, the effective mass for the fluctuations satisfies $m_{\rm eff}(t)\ll H(t)$ at early times, where $m_{\rm eff}$ is given in Eq. (34). While most modes are inside the Hubble radius, with $k/a>H\gg m_{\rm eff}$, their energy density therefore evolves like a gas of (nearly) massless particles, with an equation of state like radiation, $\delta\rho_{(2)}(t)\propto a^{-4}(t)$. At later times, after the modes have crossed outside the Hubble radius and their amplitudes have frozen, $\delta\rho_{(2)}(t)$ becomes constant.

Next we consider the impact of these large initial inhomogeneities on the evolution of the Hubble parameter $H(t)$, shown in Fig. 6. The figure shows $H(t)$ for the same initial values $\varphi_{0}=25\,M_{\rm pl}$ and $\dot{\varphi}_{0}=-0.25\,M_{\rm pl}^{2}$, when we neglect fluctuations (blue); when we initialize the fluctuations in the Bunch-Davies state, with ${\cal C}_{\rm BD}=2$ (yellow); and for a particular simulation in which we initialized the system with large initial fluctuations, ${\cal C}=20$ (pink). The energy density associated with $\varphi$, $\bar{\rho}(t)$, is dominated at early times by the kinetic energy of $\varphi$, and hence it decays as $\bar{\rho}(t)\propto a^{-6}(t)$. When we neglect fluctuations, we therefore find $H(t)\propto a^{-3}(t)$ at early times. On the other hand, for large initial fluctuations with ${\cal C}=20$ and hence $\delta\rho_{(2)}(t_{0})\simeq\bar{\rho}(t_{0})$, we find $H(t)\propto[\bar{\rho}(t)+\delta\rho_{(2)}(t)]^{1/2}\propto a^{-2}(t)$ at early times, while most fluctuations remain sub-Hubble and $\delta\rho_{(2)}(t)$ decays like radiation. For fluctuations that begin in the Bunch-Davies initial state, with ${\cal C}_{\rm BD}=2$ and $\delta\rho_{(2)}(t_{0})\sim 0.1\,\bar{\rho}(t_{0})$, we find an evolution for $H(t)$ intermediate between these two cases. (The authors of Ref. East et al. (2016) likewise found the volume-averaged quantities $\rho_{\rm avg}(t)\propto[a_{\rm avg}(t)]^{-4}$ and $H_{\rm avg}(t)\propto[a_{\rm avg}(t)]^{-2}$ at early times in their numerical simulations of large-field models with significant initial inhomogeneities. See also Ref. Chowdhury et al. (2019).)

The system begins to inflate, with $\ddot{a}>0$, once $\epsilon<1$. For the set of initial conditions we consider here, inflation begins by $N\sim 2$, and the system enters a phase of slow-roll inflation ($\epsilon<0.1$) by $N\sim 3$. During slow-roll, $\bar{\rho}$ is dominated by $V(\varphi)$, while $\delta\rho_{(2)}$ continues to redshift like radiation until most of the modes have crossed outside the Hubble radius, by $N\sim 7$ (as shown in Fig. 5). Hence after slow-roll inflation begins, the system evolves with $\bar{\rho}\simeq V(\varphi)\gg\delta\rho_{(2)}$ and $H(t)$ settles onto a nearly constant value at $H_{\rm infl}\sim 10^{-3}\,M_{\rm pl}$.

As shown in Fig. 7, for the case with large initial fluctuations, the system reaches the slow-roll inflationary attractor while most of the power in fluctuations remains inside the Hubble radius. In the presented simulation, when slow-roll inflation begins (with $\epsilon\leq 0.1$), all of the modes that had begun inside the Hubble radius still remain inside the Hubble radius. Modes with comoving wavenumber up to the UV regulator scale $\kappa=5\bar{H}_{0}=M_{\rm pl}/2$ remain inside the Hubble radius for another 3 efolds after slow-roll begins, and the shortest-wavelength mode in the simulation, with $k_{\rm max}=4\kappa=2M_{\rm pl}$, crosses outside the Hubble radius more than 4 efolds after the system reaches the slow-roll attractor. Hence the early-time dynamics, during which the system enters a phase of slow-roll inflation, occurs with substantial inhomogeneity on sub-Hubble length scales. For this set of initial conditions, in other words, inflation is robust even amid large initial inhomogeneities and with initial conditions for $\varphi(t)$ dominated by kinetic energy.

Furthermore, as shown in Fig. 8, we find that for this set of initial conditions inflation actually persists considerably longer when we include large initial inhomogeneities ($N_{\rm infl}\simeq 69$ efolds of inflation) than when we ignore inhomogeneities ($N_{\rm infl}\simeq 54$ efolds of inflation). (The authors of Ref. Easther et al. (2014) found similar examples in their study of the onset of inflation in multifield models, when fluctuations $\delta\phi(x^{\mu})$ and $\delta\psi(x^{\mu})$ of the two fields were included.) Remarkably, significant initial inhomogeneities extended the duration of inflation in this case, enabling this set of initial conditions $\left(\varphi_{0},\dot{\varphi}_{0}\right)$ to yield sufficient inflation, with $N_{\rm infl}>65$.

IV.2 Trajectories in Phase Space

We can understand the nontrivial effects of large initial inhomogeneities and their nonlinear backreaction on the evolution of $\varphi(t)$ and $H(t)$ by examining the evolution of the system through the phase space $(\varphi(t),\dot{\varphi}(t))$. As discussed in Refs. Remmen and Carroll (2013); Chowdhury et al. (2019), for single-field models and vanishing spatial curvature $K$, the variables $\varphi(t)$ and $\dot{\varphi}(t)$ define an effective phase space for the evolution of the spatially homogeneous system. Obviously $(\varphi(t),\dot{\varphi}(t))$ no longer serves as a proper phase space for the full dynamical system when we incorporate the coupled degrees of freedom $\delta\hat{\phi}(x^{\mu})$ and $\hat{\Psi}(x^{\mu})$, but studying the behavior of the system within $(\varphi(t),\dot{\varphi}(t))$ facilitates comparison with the case in which we neglect fluctuations.

In Fig. 9, we plot the evolution of the system $(\varphi(t),\dot{\varphi}(t))$ when we fix $\varphi_{0}=25\,M_{\rm pl}$ and select $\dot{\varphi}_{0}=\pm 0.25\,M_{\rm pl}^{2}$. Dashed lines show the evolution of the system for $\epsilon>0.1$, before the system reaches the slow-roll inflationary attractor, and solid lines show the evolution once the system has entered slow-roll. We consider cases in which we neglect fluctuations (blue); in which we initialize the fluctuations in the Bunch-Davies state, ${\cal C}_{\rm BD}=2$ (yellow); and in which we initialize the system with large fluctuations, ${\cal C}=20$ (pink). As Fig. 9 makes clear, the value of the field when the system reaches the slow-roll attractor, $\varphi(t_{\rm sr})$, depends on the magnitude of the initial inhomogeneities. In particular, the field $\varphi$ traverses a shorter distance before arriving at the slow-roll attractor when we incorporate fluctuations, compared to when we neglect fluctuations: either less far “up the hill” toward higher values of the potential for an initial field velocity $\dot{\varphi}_{0}>0$, or less far “down the hill” for $\dot{\varphi}_{0}<0$. This effect becomes more pronounced as the size of initial inhomogeneities increases.

We can make sense of this result analytically, using the scaling relations for $H(t)$ identified in the previous subsection. For these initial conditions and the coupling $\lambda=10^{-10}$, the system begins with the kinetic energy in $\varphi$ greatly exceeding the potential energy, and hence $\bar{\rho}(t)\propto a^{-6}(t)$ before the system enters slow-roll. This is equivalent to

[TABLE]

at early times. In that limit, Eq. (17) reduces to $\ddot{\varphi}+3H\dot{\varphi}\simeq 0$, which we may integrate as

[TABLE]

When we neglect fluctuations, $H(t)\propto a^{-3}(t)$ at early times, and we find

[TABLE]

where the subscript “$h$” indicates evolution of the homogeneous system. On the other hand, when we include large initial fluctuations with $\delta\rho_{(2)}(t_{0})\simeq\bar{\rho}(t_{0})$, then $H(t)\propto a^{-2}(t)$ at early times, which yields

[TABLE]

where the subscript “$q$” indicates the evolution of $\varphi(t)$ when we incorporate effects from the coupled quantum fluctuations. Clearly the field $\varphi$ will traverse a greater distance during early times when the fluctuations are neglected, as in Eq. (50), than when their effects are included, as in Eq. (51). The solutions for $\dot{\varphi}(N)$ in Eq. (48) and for $\varphi(N)$ in Eqs. (50) or (51) closely match the trajectories shown in Fig. 9 for the relevant cases, even though the curves in Fig. 9 come from our full numerical simulations.

We may use Eqs. (50)-(51) to estimate the values $\varphi_{h}(N_{\rm sr})$ and $\varphi_{q}(N_{\rm sr})$ at the time $N_{\rm sr}$ when the system reaches the slow-roll attractor. We estimate $N_{\rm sr}$ by setting $\dot{\varphi}^{2}(N_{\rm sr})/2=V(\varphi(N_{\rm sr}))$. For $\lambda=10^{-10}$ and $(\varphi_{0},\dot{\varphi}_{0})=(25\,M_{\rm pl},-0.25\,M_{\rm pl}^{2})$, we have $\bar{H}_{0}=[\dot{\varphi}_{0}^{2}/(6M_{\rm pl}^{2})]^{1/2}$ for the homogeneous case and $H_{0}=\sqrt{2}\,\bar{H}_{0}$ for the case with $\delta\rho_{(2)}(t_{0})\simeq\bar{\rho}(t_{0})$. These yield $\varphi_{h}(N_{\rm sr})=21.4\,M_{\rm pl}$ and $\varphi_{q}(N_{\rm sr})=23.7\,M_{\rm pl}$, again closely matching the numerical results shown in Fig. 9.

After the system reaches the slow-roll attractor, $\bar{\rho}(t)\sim{\rm constant}$ while $\delta\rho_{(2)}(t)$ continues to redshift like radiation until most of the modes have crossed outside the Hubble radius, so the dynamics become dominated by $\bar{\rho}\gg\delta\rho_{(2)}$. In that regime, we may use the usual slow-roll approximation to estimate the duration of inflation,

[TABLE]

where $\varphi_{\rm sr}=\varphi(N_{\rm sr})$, and $\varphi_{\rm end}=\varphi(N_{\rm end})$ is determined by the condition $\epsilon(N_{\rm end})=1$. Again using the usual slow-roll estimate for late times, $\epsilon\simeq(M_{\rm pl}^{2}/2)(V^{(1)}(\varphi)/V(\varphi))^{2}$, we find $\varphi_{\rm end}=\sqrt{8}\,M_{\rm pl}$. Given our estimates of $\varphi_{q}(N_{\rm sr})$ and $\varphi_{h}(N_{\rm sr})$, we then find $N_{\rm infl}=69.2$ efolds of inflation when we incorporate large initial quantum fluctuations, and $N_{\rm infl}=56.2$ efolds when we neglect fluctuations — a close match to the behavior shown for the full numerical results in Fig. 8.

We can thus understand the most significant effect of the coupled fluctuations on the evolution of the system. Large fluctuations raise the initial value of the Hubble parameter compared to the case with no fluctuations, $H_{0}>\bar{H}_{0}$, thereby increasing the initial Hubble drag on the field $\varphi(t)$. Even more significant, backreaction from the fluctuations changes the scaling of $H(t)$ with $a(t)$ at early times, slowing the rate at which $H(t)$ falls, which further increases the effect of Hubble drag on the evolution of $\varphi(t)$. The backreaction dampens $\varphi$’s motion as the system evolves toward the slow-roll inflationary attractor, such that $|\varphi_{q}(N_{\rm sr})-\varphi_{0}|<|\varphi_{h}(N_{\rm sr})-\varphi_{0}|$. Once the system reaches the attractor, the duration of inflation will be governed by the value $\varphi_{q}(N_{\rm sr})$. For an initial velocity “up the hill,” with $\dot{\varphi}_{0}>0$, $\varphi_{q}(N_{\rm sr})<\varphi_{h}(N_{\rm sr})$, and the system will spend less time evolving along the inflationary attractor than in the absence of fluctuations. For an initial velocity “down the hill,” with $\dot{\varphi}_{0}<0$, $\varphi_{q}(N_{\rm sr})>\varphi_{h}(N_{\rm sr})$, and the system will spend more time evolving along the inflationary attractor than in the absence of fluctuations.

IV.3 Phase Space of Initial Conditions

We turn now to discuss the effects of the coupled fluctuations on the evolution of the system across the phase space of initial conditions $(\varphi_{0},\dot{\varphi}_{0})$ for $\lambda=10^{-10}$, as we vary $12\,M_{\rm pl}\leq\varphi_{0}\leq 30\,M_{\rm pl}$ and $-0.25\,M_{\rm pl}^{2}\leq\dot{\varphi}_{0}\leq 0.25\,M_{\rm pl}^{2}$. To investigate the phase space behavior for the perturbed initial conditions, we construct averages from the 32 samples evolved at each point in $(\varphi_{0},\dot{\varphi}_{0})$. Figure 10 shows the average value of $\delta\rho_{(2)}(t_{0})$ at each grid point compared to $\delta\rho_{(2)}(t_{0})$ for Bunch-Davies initial conditions, confirming that for the ranges of coefficients in Eq. (44), we find initial energy densities about ten times greater than for the Bunch-Davies state. In Fig. 11, we plot the average of the initial value $\Psi_{\rm rms}(t_{0})\equiv[\langle\hat{\Psi}^{2}(t_{0})\rangle]^{1/2}$ at each grid point, confirming that for large quantum fluctuations, initialized such that ${\cal C}\simeq 20$, the system begins with $|\Psi(t_{0},{\bf x})|\lesssim 0.5$.

Across $(\varphi_{0},\dot{\varphi}_{0})$, when the system begins with large initial inhomogeneities, the system reaches the slow-roll inflationary attractor ($\epsilon\leq 0.1$) while significant power remains in fluctuations on sub-Hubble scales. Figure 12 shows the average value of the ratio $N_{\rm sr}/N_{\kappa}$ at each grid point for fluctuations that begin with ${\cal C}\simeq 20$, where $N_{\kappa}$ is the time when the mode with comoving wavenumber equal to the UV regulator scale $\kappa$ crosses outside the Hubble radius, $\kappa=aH$. In all simulations, $N_{\rm sr}/N_{\kappa}\leq 0.55\pm 0.06$. The ratio drops to zero at $\dot{\varphi}_{0}=0$, because at those locations in phase space the system begins on the slow-roll attractor, and hence $N_{\rm sr}=0$.

Results for the duration of inflation $N_{\rm infl}$ across these cases are shown in Fig. 13. Consider first the case in which the effects of the coupled fluctuations are neglected. For a given value $\varphi_{0}$, large initial velocities $\dot{\varphi}_{0}>0$ prolong the duration of inflation compared to the case with $\dot{\varphi}_{0}=0$, whereas large initial velocities $\dot{\varphi}_{0}<0$ decrease the duration of inflation. When initial inhomogeneities are included they backreact on $H(t)$, increasing the effect of Hubble drag on $\varphi(t)$, thereby affecting the shape of the contours of constant $N_{\rm infl}$ within $(\varphi_{0},\dot{\varphi}_{0})$. In particular, the effects of large $|\dot{\varphi}_{0}|$ are more quickly damped than when the fluctuations are neglected, so that the field $\varphi$ spends less time evolving along the slow-roll inflationary attractor than the corresponding case without fluctuations for $\dot{\varphi}_{0}>0$, and more time along the attractor for $\dot{\varphi}_{0}<0$.

Figure 14 shows the critical line in $(\varphi_{0},\dot{\varphi}_{0})$, to the right of which the system yields $N_{\rm infl}>65$ efolds of inflation, for the homogeneous system (when we neglect fluctuations), for fluctuations in the Bunch-Davies initial state with ${\cal C}_{\rm BD}=2$, and for larger initial fluctuations with ${\cal C}\simeq 20$. For the cases with large initial fluctuations, we plot the critical line based on averaging across $32$ simulations per grid point, as well as $\pm 2\sigma$ contours.

Consistent with the analysis in the previous subsection, the effects of large initial fluctuations are symmetric for $\pm|\dot{\varphi}_{0}|$, and most significant for large $|\dot{\varphi}_{0}|$. As $|\dot{\varphi}_{0}|$ increases, the initial value $\bar{H}_{0}$ increases; the greater initial energy scale $\bar{H}_{0}$, in turn, yields more initial energy density in fluctuations, $\delta\rho_{(2)}(t_{0})$, which raises $H_{0}>\bar{H}_{0}$ and affects the scaling of $H(t)$ with $a(t)$ at early times. The larger fluctuations also seed larger initial inhomogeneities, $\Psi(t_{0},{\bf x})$.

Although the effects of the nonlinear dynamics of the fluctuations are most pronounced at large $|\dot{\varphi}_{0}|$, the volume of the (projected) phase space of initial conditions $(\varphi_{0},\dot{\varphi}_{0})$ that yields sufficient inflation is conserved. The same volume of the $(\varphi_{0},\dot{\varphi}_{0})$ phase space that yields sufficient inflation which is lost in the region with $\dot{\varphi}_{0}>0$, compared to the homogeneous case, is gained in the region with $\dot{\varphi}_{0}<0$. This provides a useful quantitative measure of the robustness of single-field inflation to large initial inhomogeneities.

IV.4 Varying the Coupling Constant

The main impact of lowering the coupling constant from $\lambda=10^{-10}$ to $\lambda=10^{-12}$ is to increase the ratio $H_{0}/H_{\rm infl}$ from ${\cal O}(10^{2})$ to ${\cal O}(10^{3})$. That prolongs the time during which the enhanced Hubble drag from the coupled fluctuations affects the evolution of $\varphi(t)$, compared to the case in which the fluctuations are neglected. We divided the phase space of initial conditions $(\varphi_{0},\dot{\varphi}_{0})$ into the same grid as for the $\lambda=10^{-10}$ case, and considered cases in which we neglected fluctuations, began with fluctuations in the Bunch-Davies initial state, and began with larger initial fluctuations, with the random coefficients for each mode $\delta\phi_{n\ell m}(t_{0})$ and $\delta\dot{\phi}_{n\ell m}(t_{0})$ drawn from the ranges in Eq. (44). Because the contributions to $\bar{\rho}(t_{0})$ and $\delta\rho_{(2)}(t_{0})$ that are proportional to $\lambda$ remain subdominant for $\lambda=10^{-12}$ as for $\lambda=10^{-10}$, the initial values for $(\delta\rho_{(2)}(t_{0}))_{\rm avg}/(\delta\rho_{(2)}(t_{0}))_{\rm BD}$ and for $\Psi_{\rm rms}(t_{0})$ are little changed from the results shown in Figs. 10 and 11. Likewise, we again find that throughout $(\varphi_{0},\dot{\varphi}_{0})$, $N_{\rm sr}/N_{\kappa}\leq 0.53\pm 0.04$, as shown in Fig. 15.

Figure 16 shows contours of constant $N_{\rm infl}$ in $(\varphi_{0},\dot{\varphi}_{0})$ with $\lambda=10^{-12}$ for the three cases of interest: no fluctuations, Bunch-Davies intial state (${\cal C}_{\rm BD}=2$), and larger initial fluctuations (with ${\cal C}\simeq 20$). In Fig. 17 we plot the critical line in $(\varphi_{0},\dot{\varphi}_{0})$ that yields $N_{\rm infl}=65$ efolds of inflation for each of the three cases. For the cases with large initial fluctuations, we plot the critical line based on averaging across $32$ simulations per grid point, as well as $\pm 2\sigma$ contours. The results are comparable to the case with $\lambda=10^{-10}$. Although the effects of the coupled fluctuations are more dramatic with the smaller coupling — because the system takes more time to reach the inflationary attractor, and hence the enhanced Hubble drag in the presence of coupled fluctuations acts longer on the evolution of $\varphi(t)$ — the effects across the projected phase space $(\varphi_{0},\dot{\varphi}_{0})$ are once again symmetrical for $\pm|\dot{\varphi}_{0}|$, indicating that the total volume of the $(\varphi_{0},\dot{\varphi}_{0})$ phase space that yields sufficient inflation is conserved.

V Conclusions

In this paper we have analyzed the onset of inflation for a simple single-field model, $V(\phi)=\lambda\phi^{4}/4$, when the system begins with significant inhomogeneities on length-scales shorter than the initial Hubble radius. We incorporate certain nonlinear interactions among the coupled degrees of freedom by using the nonperturbative Hartree approximation, which resums an infinite set of Feynman diagrams involving the self-interacting quantum fluctuations $\delta\hat{\phi}(x^{\mu})$ to construct a dressed propagator. By initializing the system in an excited state, with the energy density of fluctuations approximately ten times greater than in the minimum-energy Bunch-Davies state, our simulations begin with significant spatial inhomogeneities, parameterized by the scalar metric perturbation $|\Psi(t_{0},{\bf x})|\lesssim 0.5$.

The energy density in fluctuations, $\delta\rho_{(2)}(t)$, backreacts on the evolution of the Hubble parameter $H(t)$, affecting the scaling of $H(t)$ with $a(t)$. This backreaction, in turn, leads to increased Hubble drag on the evolution of the vacuum expectation value of the field, $\varphi(t)$, affecting how quickly $\varphi(t)$ arrives at the slow-roll inflationary attractor, compared to the case in which one neglects fluctuations.

The impact of inhomogeneities on the evolution of the system is largest for initial conditions that yield the greatest initial value of $H(t_{0})$, since the initial energy density in fluctuations scales as $H^{4}(t_{0})$. Compared to those regions of $(\varphi_{0},\dot{\varphi}_{0})$ that yield $N_{\rm infl}\geq 65$ efolds of inflation when one neglects inhomogeneities, we find some regions that fail to yield sufficient inflation when we incorporate inhomogeneities, and an equal volume of regions that succeed in producing $N_{\rm infl}\geq 65$ but which had failed to do so in the absence of inhomogeneities. (See also Ref. Easther et al. (2014).) In other words, the total volume of the space $(\varphi_{0},\dot{\varphi}_{0})$ that yields sufficient inflation is conserved when we incorporate nonlinear backreaction from inhomogeneities, compared to the case in which we neglect inhomogeneities.

Although we have analyzed the system numerically for a specific form of $V(\phi)$, the arguments about the robustness of inflation for such large-field models do not depend on our choice of $V(\phi)$. All that enters into our semi-analytic argument is that the system can begin with large initial quantum fluctuations, such that $\delta\rho_{(2)}(t_{0})$ is comparable to (or greater than) the initial energy density associated with the vacuum expectation value, $\bar{\rho}(t_{0})$. For weakly coupled models — as required for large-field inflation, in order to produce a spectrum of primordial density perturbations consistent with observations — we generically expect that $\delta\rho_{(2)}(t)$ will evolve at early times with an equation of state like that of a gas of (nearly) massless particles, scaling as $\delta\rho_{(2)}(t)\propto a^{-4}(t)$ while most of the power in fluctuations is on sub-Hubble scales. This behavior for the fluctuations contrasts with the scaling of $\bar{\rho}(t)$ at early times, when the energy density associated with $\varphi(t)$ is dominated by kinetic energy, such that $\bar{\rho}(t)\propto a^{-6}(t)$. The backreaction of $\delta\rho_{(2)}(t)$ on $H(t)$ ensures that $\varphi(t)$ will traverse less distance en route to the slow-roll inflationary attractor than in the absence of inhomogeneities, thereby accounting for the differences we observe in the duration of inflation, $N_{\rm infl}$.

In our numerical analysis we initialize fluctuations at $t_{0}$ across the range of wavenumbers from $k_{\rm min}/a(t_{0})\lesssim 2H_{0}/3$ up to $k_{\rm max}/a(t_{0})=30\,k_{\rm min}/a(t_{0})\sim M_{\rm pl}$. Including any modes with $k>k_{\rm max}$ in the spectrum would lead to trans-Planckian ambiguities. One could nonetheless imagine initializing additional modes at later times — to simulate the notion that modes which had begun with $k/a(t_{0})\gg M_{\rm pl}$ at $t_{0}$ later redshifted to $k/a(t)\leq M_{\rm pl}$ — but we do not expect such additional, trans-Planckian modes to change the qualitative behavior of the system. In our current framework, the system consistently arrives at the slow-roll inflationary attractor while most of the initial power in sub-Hubble fluctuations remains within the Hubble radius. The energy density associated with any modes that might be initialized at later times $t>t_{0}$ would be less than $\delta\rho_{(2)}(t_{0})$, since $H(t)<H(t_{0})$ and $\delta\rho_{(2)}(t)\propto H^{4}(t)$. Moreover, after a few efolds of slow-roll inflation, we expect that any newly initialized fluctuations should begin in the Bunch-Davies state, rather than in the more-energetic initial states that we consider here.

For next steps, we aim to generalize our formalism to include the evolution of systems with nonzero spatial curvature $K$, to consider small-field as well as large-field models, and to extend the formalism to multifield models (akin to Ref. Easther et al. (2014), but incorporating the coupled metric perturbations). Across each of these cases, we believe the approach we have developed here can complement the computationally intensive numerical-relativity approaches of Refs. East et al. (2016); Clough et al. (2017, 2018).

Appendix A Discrete Spectrum for Mode Functions

As described in Section III, we expand quantized field fluctuations $\delta\hat{\phi}\left(x^{\mu}\right)$ and metric perturbations $\hat{\Psi}\left(x^{\mu}\right)$ in eigenfunctions $Z_{n\ell m}\left(\mathbf{x}\right)$ of the comoving spatial Laplacian. Within a comoving spatial volume of finite size, Eq. (12) then takes the form

[TABLE]

with positive integer $n\geq 1$. We restrict attention to a finite sphere of comoving radius $R$ within a spatially flat FLRW background spacetime ($K=0$). We select Dirichlet boundary conditions $Z_{n\ell m}(R,\theta,\phi)=0$, which causes $\mathbf{\nabla}^{2}$ to have a negative definite spectrum as desired, and for the resulting Sturm-Liouville system to yield a complete, orthonormal basis:

[TABLE]

In spherical polar coordinates, Eq. (53) becomes

[TABLE]

where

[TABLE]

The eigenfunctions of the operator ${\bf L}^{2}$ are the familiar spherical harmonics $Y_{\ell m}(\theta,\phi)$, which satisfy ${\bf L}^{2}Y_{\ell m}=\ell(\ell+1)\,Y_{\ell m}$. Solutions to Eq. (55) may then be written in the form

[TABLE]

where $N_{n\ell m}$ is a normalization constant, $j_{\ell}(x)$ is a spherical Bessel function of order $\ell$, and the boundary conditions require

[TABLE]

The requirement of Eq. (58) yields a discrete spectrum of allowable wavenumbers,

[TABLE]

where, as noted below Eq. (33), $x_{n\ell}$ is the $n$th zero of the Bessel function $j_{\ell}(x)$ for $n\geq 1$. Making use of the orthonormality properties of the $Y_{\ell m}(\theta,\phi)$, Eq. (54) then becomes

[TABLE]

Upon using Eqs. (11.49) and (11.50) of Ref. Arfken and Weber (2005), we find

[TABLE]

with (as usual) $j_{\ell}(k_{n\ell}R)=0$. Note that the normalization $N_{n\ell m}$ is independent of $m$. The basis functions $Z_{n\ell m}(r,\theta,\phi)$ provide spectral convergence to any spatial configuration as we take $n_{\rm max}\to\infty$.

Given the spatial eigenfunctions $Z_{n\ell m}({\bf x})$ and the properties of the operators $\hat{a}_{n\ell m}$, $\hat{a}^{\dagger}_{n\ell m}$, we may evaluate various two-point functions. For example, we have

[TABLE]

The term $\langle(\delta\hat{\phi})^{2}\rangle$, which appears throughout the equations of motion in the Hartree approximation, is Eq. (62) in the limit ${\bf x}\rightarrow{\bf y}$. As our slicing of spacetime has been chosen such that spatial slices are homogeneous and isotropic (this holds even inside the finite sphere due to completeness), we may evaluate Eq. (62) at ${\bf x}={\bf y}={\bf 0}$. From

[TABLE]

we note that only the $\ell=0$ mode remains nonzero at the origin, with amplitude $j_{0}(0)=1$. We further note that $Y_{00}(\theta,\phi)=1/\sqrt{4\pi}$ for $\ell=m=0$, and hence we find

[TABLE]

as in Eq. (38), prior to applying the UV regularization. The same steps yield $\langle(\delta\dot{\phi})^{2}\rangle$ as in Eq. (38).

The remaining two-point function of interest is the contribution to $\delta\rho_{(2)}$ from the spatial gradients. We find

[TABLE]

Using Eq. (63), we see that only the term with $\ell=1$ will contribute to the first term in brackets within Eq. (65). From the properties of $Y_{1m}(0,\phi)$ and $N_{n\ell m}$ in Eq. (61), we then find

[TABLE]

For the last two terms in brackets in Eq. (65), we require that the terms be well-behaved in the vicinity of the origin. Given $h^{\theta\theta}=1/r^{2}$ and the properties of the $Y_{\ell m}(\theta,\phi)$, the only contribution to the second term in brackets that will remain regular (and nonzero) near ${\bf x}\rightarrow{\bf 0}$ comes from $\ell=1$, and we find

[TABLE]

Proceeding similarly, we find

[TABLE]

Combining Eqs. (65) - (68), we find

[TABLE]

as in Eq. (38).

Appendix B Initial Conditions for the Field Fluctuations

To establish the initial conditions for the mode functions $\delta\phi_{n\ell m}(t_{0})$ and $\delta\dot{\phi}_{n\ell m}(t_{0})$, we first consider the case of $R\rightarrow\infty$, for which the spectrum of allowable wavenumbers is continuous, with $0\leq k<\infty$. As noted in Section III, within the Hartree approximation higher-order interaction terms among the fluctuations $\delta\hat{\phi}$ are replaced by an effective mass, so we may proceed by following most of the steps for quantizing a free, massive scalar field in FLRW spacetime.

The equal-time commutation relation for a free scalar field stipulates

[TABLE]

where

[TABLE]

is the momentum canonically conjugate to $\delta\hat{\phi}$. Upon expanding $\delta\hat{\Pi}(x^{\mu})$ in a series of complex mode functions $\delta\Pi_{k\ell m}(t)$ and creation and annihilation operators akin to Eq. (10), and using the commutation relations for $\hat{a}_{k\ell m}$ and $\hat{a}^{\dagger}_{k\ell m}$ in Eq. (13), Eq. (70) imposes a constraint on the mode functions:

[TABLE]

We rescale the modes $\delta\phi_{k\ell m}(t)=v_{k\ell m}(t)/a(t)$ and introduce conformal time, $d\tau\equiv dt/a(t)$, so that

[TABLE]

where primes denote $d/d\tau$ and ${\cal H}\equiv a^{\prime}/a$. Then Eq. (72) becomes

[TABLE]

a Wronskian condition that we will use when setting initial conditions at $\tau_{0}$.

In terms of $v_{k\ell m}(\tau)$, Eq. (18) takes the form

[TABLE]

where the source term $S_{k\ell m}$ depends on $\Psi_{k\ell m}$ and $\Psi_{k\ell m}^{\prime}$, and the frequency is given by

[TABLE]

with the effective mass $m_{\rm eff}^{2}$ given in Eq. (34). The effective mass is suppressed by a small coupling constant $\lambda$ and hence we expect it to remain subdominant around $\tau_{0}$, while, on dimensional grounds, we expect $a^{\prime\prime}/a\sim{\cal O}(H^{2})$ around $\tau_{0}$. When we numerically evolve the modes within a region of finite comoving radius $R$ in our simulation, we introduce an infrared cut-off given by $k_{\rm min}=\pi/R$, with $R\sim 1/[a(\tau_{0})H(\tau_{0})]$. Hence even for the longest-wavelength modes in our simulation, we expect $k^{2}\gtrsim\{a^{2}(\tau_{0})\,m_{\rm eff}^{2}(\tau_{0}),a^{\prime\prime}(\tau_{0})/a(\tau_{0})\}$. For setting initial conditions, we therefore approximate $\omega_{k}(\tau_{0})\sim k$.

We do not assume that the system has attained its minimum energy state at the initial time $\tau_{0}$, and hence we consider initial conditions for the modes that could depart from the usual Bunch-Davies vacuum state. (Note that as we are describing a weakly interacting system, the Bunch-Davies state that we construct is only an approximation of the lowest-energy state, which we do not attempt to find.) We parameterize the initial conditions for the modes as

[TABLE]

where $\alpha_{k\ell m}$, $\beta_{k\ell m}$, $\gamma_{k\ell m}$, and $\delta_{k\ell m}$ are each real-valued, dimensionless constants. (In the limit $\omega_{k}\rightarrow k$, the Bunch-Davies state corresponds to $\alpha_{k\ell m}=\gamma_{k\ell m}=1$, $\beta_{k\ell m}=\delta_{k\ell m}=0$.) The Wronskian condition of Eq. (74) then becomes

[TABLE]

Without loss of generality, we may set $\beta_{k\ell m}=0$ for all modes; then Eq. (78) fixes $\alpha_{k\ell m}=1/\gamma_{k\ell m}$. Setting $a(\tau_{0})=1$, Eqs. (73) and (77) yield

[TABLE]

We use $\bar{H}_{0}$ (rather than $H_{0}$) in Eq. (79) because, as a practical matter, we must first select a point in $(\varphi_{0},\dot{\varphi}_{0})$ and then initialize the fluctuations. The parameter $\bar{H}_{0}$, defined in Eq. (36), is determined by the selection $(\varphi_{0},\dot{\varphi}_{0})$; only after the fluctuations are initialized can we evaluate $\delta\rho_{(2)}(t_{0})$ and thereby include their contribution to $H_{0}$. (For the regimes of interest, $H_{0}/\bar{H}_{0}\lesssim 1.4$, so the difference is not significant.) When we evolve the system within a finite sphere of comoving radius $R$, the continuous variable $k$ in Eqs. (78) and (79) is replaced by the discrete spectrum of wavenumbers $k_{n\ell}$ with $n\geq 1$, which yields Eq. (35).

Given these initial conditions for the modes, we may estimate quantities of interest such as $\langle(\delta\hat{\phi}(t_{0}))^{2}\rangle$ and $\delta\rho_{(2)}(t_{0})$. We derive our expressions in the continuum limit using the UV regulator function $F_{n\ell}(\kappa,R)$ of Eq. (39). For the initial value of the two-point function, we begin with the expansion for $\delta\hat{\phi}(x^{\mu})$ in Eq. (10), use Eq. (13), make use of the properties of $Z_{k\ell m}({\bf x})$ at ${\bf x}\rightarrow{\bf 0}$, and substitute $\delta\phi_{k\ell m}(t_{0})$ from Eq. (79) to write

[TABLE]

The third line follows because the random coefficients $\alpha_{k\ell m}$ vary independently of $k$, so we may take the term $\alpha_{k00}^{2}$ outside the integral and replace it with an average value; moreover, since we draw the variables $\gamma_{k\ell m}=1/\alpha_{k\ell m}$ from the same distribution for all $(k\ell m)$, the same average value $(\alpha^{2})_{\rm avg}$ holds for any $k$, $\ell$, and $m$. Upon using $\kappa=b\bar{H}_{0}$, we arrive at Eq. (41).

We proceed similarly to estimate $\delta\rho_{(2)}(t_{0})$. Using Eqs. (10), (13), (22), (39), and (79), we find

[TABLE]

where we have used $a(t_{0})=1$ and assumed that $m_{\rm eff}^{2}(t_{0})\ll H_{0}^{2}$ in the regimes of interest, so that the contributions to $\delta\rho_{(2)}(t_{0})$ are dominated by the kinetic energy and spatial-gradient energy of the fluctuations. Given $\kappa=b\bar{H}_{0}>\bar{H}_{0}$, we find Eq. (42).

The magnitudes of $\langle(\delta\hat{\phi}(t_{0}))^{2}\rangle$ and $\delta\rho_{(2)}(t_{0})$ depend on the random variables $\alpha_{n\ell m}$, $\gamma_{n\ell m}$, and $\delta_{n\ell m}$ for each mode. We draw $\gamma_{n\ell m}$ and $\delta_{n\ell m}$ from flat distributions within specific ranges; once $\gamma_{n\ell m}$ is selected, $\alpha_{n\ell m}$ is fixed by the quantization condition to be $\alpha_{n\ell m}=1/\gamma_{n\ell m}$. If we draw $\gamma_{n\ell m}$ from the range $\{A,B\}$ and $\delta_{n\ell m}$ from the range $\{-C,C\}$, then we expect

[TABLE]

In order to study large initial fluctuations, with ${\cal C}\equiv\left(\alpha^{2}+\gamma^{2}+\delta^{2}\right)_{\rm avg}\simeq 20$, we therefore use the ranges in Eq. (44).

Appendix C Initializing the Metric Perturbations

In this appendix, we discuss subtleties in the initialization of $\Psi_{n\ell m}(t_{0})$, and the techniques we use to avoid these issues numerically.

Given an instantiation of $\delta\phi_{n\ell m}(t_{0})$ and $\delta\dot{\phi}_{n\ell m}(t_{0})$, the metric perturbations $\Psi_{n\ell m}(t_{0})$ are initialized using Eq. (29). For our numerical simulations, we evolve the system within a sphere of comoving radius $R$, within which Eq. (29) takes the form

[TABLE]

As $\dot{H}$ tends to be large and negative at early times, it is common for there to exist a value for $k_{n\ell}$ at which the coefficient of $\Psi_{n\ell m}$ vanishes at $t_{0}$: we dub the $k_{n\ell}$ value at which this occurs the “$\Psi$ pole.” For $k_{n\ell}$ near the pole, $\Psi_{n\ell m}(t_{0})$ would be initialized with an artificially large initial amplitude, which can cause a single mode to dominate the nonlinear contributions.

The pole exists as an artifact of the longitudinal gauge, and arises because the four-dimensional phase space ($\delta\phi_{n\ell m}$, $\delta\dot{\phi}_{n\ell m}$, $\Psi_{n\ell m}$ and $\dot{\Psi}_{n\ell m}$) has two constraints, Eqs. (21) and (29). When $k_{n\ell}$ is on the $\Psi$ pole, the constraint surface in phase space forces $\delta\phi_{n\ell m}$ and $\delta\dot{\phi}_{n\ell m}$ to be strictly proportional to each other. In this situation, we do not lose a degree of freedom (the constraint surface remains two-dimensional); it is simply that $\delta\phi_{n\ell m}$ and $\delta\dot{\phi}_{n\ell m}$ cease to be appropriate coordinates on the constraint surface.

We emphasize that this does not arise due to our use of the Hartree approximation, but rather through imposing canonical commutation relations on $\delta\phi_{n\ell m}$ and $\delta\dot{\phi}_{n\ell m}$ at $t_{0}$ under the assumption that they are free fields. Within linear perturbation theory, this issue is typically overcome by using Mukhanov-Sasaki variables, $Q_{n\ell m}\equiv\delta\phi_{n\ell m}+(\dot{\varphi}/H)\Psi_{n\ell m}$, which are good coordinates on the phase space. However, the $Q_{n\ell m}$ are no longer gauge-invariant beyond linear order in perturbations. If one were to apply the (nonlinear) Hartree corrections to the $Q_{n\ell m}$, the resulting expressions that went beyond ${\cal O}(\Psi)$ would not bear any clear relationship to an expansion of the Einstein tensor $G_{\mu\nu}$ beyond ${\cal O}(\Psi)$.

The most correct way to approach this issue would be to employ Dirac’s constrained Hamiltonian formalism to the constrained phase space, but such an analysis is beyond the scope of this work. For our present purposes, we adopt methods to numerically alleviate the issue of the $\Psi$ pole.

We can compute the boost given to a particular mode $\Psi_{n\ell m}(t_{0})$ as a result of being near the pole as follows. Taking $K=0$, $a(t_{0})=1$, and identifying the pole location as

[TABLE]

we have

[TABLE]

where factors of $\bar{H}_{0}$ (the relevant scale) have been inserted to make the coefficient dimensionless. We set a threshold such that any instantiation with a boost of $\bar{H}_{0}^{2}/(k_{n\ell}^{2}-k_{0}^{2})>10$ for any $\ell=0$ or $\ell=1$ modes was too close to the pole for reliable results. For such instantiations, we randomly changed $R=1.5\pi/\bar{H}_{0}$ up or down by up to 10% (which changes the $k_{n\ell}(R)$ spectrum of the modes), and re-initialized all variables. The new position of the $k_{n\ell}$ values typically meant that no modes fell too close to the pole; if necessary, we repeated the process. We emphasize that these shifts to avoid the $\Psi$ pole are only necessary for setting initial conditions at $t_{0}$. Dynamically, we evolve the modes $\Psi_{n\ell m}(t)$ according to Eq. (21), which is well behaved.

When constructing initial surfaces for $\Psi(t_{0},{\bf x})$ (as in Fig. 1), this technique for avoiding the pole proved insufficient, as such surfaces required modes for $\ell=0,1,\ldots,\ell_{\rm max}$, where $\ell_{\rm max}$ was typically taken to be around 12. (Every mode $\Psi_{n\ell m}(t)$ evolves according to Eq. (21), though only modes $\delta\phi_{n\ell m}$ with $\ell=0,1$ contribute to the Hartree corrections.) The spectrum of all modes in such cases forms a rather dense forest, so shifting the spectrum (by adjusting $R$) simply shifts the pole from one mode to another. To ameliorate the pole in this situation, we multiplied $\bar{H}_{0}^{2}/(k_{n\ell}^{2}-k_{0}^{2})$ by a regulating function of $k_{n\ell}^{2}$ with a double zero at $k_{0}^{2}$ that takes a value of unity away from the pole. The transition width was chosen to be sufficiently narrow so that the complete factor never grew sufficiently large for one mode to dominate the plots. This is not a physical regularization, but serves to suppress modes near the pole for the purposes of visualization, and is used only to construct $\Psi(t_{0},{\bf x})$.

In order to construct an initial spatial representation of $\Psi(t_{0},r,\theta,\phi)$, we begin by sampling the initial conditions for $\delta\phi_{n\ell m}(t_{0})$ and $\delta\dot{\phi}_{n\ell m}(t_{0})$, and use these to construct $\Psi_{n\ell m}(t_{0})$ (with regulator, if necessary). We then turn to the mode expansion for $\hat{\Psi}(x^{\mu})$,

[TABLE]

This is a quantum-mechanical expansion (as evidenced by the $\hat{a}_{n\ell m}$ and $\hat{a}^{\dagger}_{n\ell m}$ operators). To construct a classical realization, we demote the operators to classical complex random variables $\tilde{a}_{n\ell m}$ and $\tilde{a}^{*}_{n\ell m}$ with the same statistics as the quantum operators.

For a given mode, the quantum operator is $\hat{\Psi}_{n\ell m}=\Psi_{n\ell m}(t)\,Z_{n\ell m}\,\hat{a}_{n\ell m}+H.c.$, so we have $\langle\hat{\Psi}_{n\ell m}\rangle=0$ and $\langle\hat{\Psi}^{2}_{n\ell m}\rangle=|Z_{n\ell m}|^{2}|\Psi_{n\ell m}(t)|^{2}$. (We only need to go to the two-point function as we have Gaussian statistics.) Demanding the same statistics for a classical function $\tilde{\Psi}_{n\ell m}=\Psi_{n\ell m}(t)\,Z_{n\ell m}\,\tilde{a}_{n\ell m}+c.c.$, we require

[TABLE]

This is achieved by constructing $\tilde{a}_{n\ell m}=x+iy$ for each mode, where $x$ and $y$ are independent Gaussian random variables with zero mean and standard deviation $\sigma_{x}=\sigma_{y}=1/2$. Constructing the spatial slice is then straightforwardly accomplished by randomly sampling $\tilde{a}_{n\ell m}$, after which Eq. (86) may be evaluated. Identical techniques can be used to construct initial surfaces for $\delta\phi(t,{\bf x})$.

Appendix D Numerical Convergence Tests

In this appendix, we discuss the accuracy of our simulations.

First, consider the accuracy of our numerical integration. The RK45 integrator uses a variable time step; early on in the integration, it takes small steps so as to capture the oscillations of the modes, and after the last mode freezes out, it takes increasingly large steps. We cap the maximum time step to have an estimated $\Delta N$ of 0.1, in terms of efolds $N=\ln a$. We employed relative and absolute tolerances of $10^{-10}$. Investigations showed this to give a good trade-off between accuracy, numerical precision and computational time. Our measurements of the duration of inflation, $N_{\rm infl}$, for a given run are completely dominated by uncertainty as to when to start measuring $N_{\rm infl}$ rather than precision from the numerical evolution.

Second, our simulations cannot simulate the entire continuous universe, but only a discrete portion thereof. Our analysis essentially consists of three tuneable parameters describing this discretization: $\kappa$, $k_{\rm max}$ and $R$. The continuum limit takes each of these parameters to approach infinity, but such lies beyond our computational power. As such, we now discuss the convergence of our simulations as each of these parameters is increased. For each of these convergence tests, we look for convergence in the two-point functions and overall efolds for Bunch-Davies simulations, for which comparisons with analytic results are possible.

The outer boundary $R$ controls both the spectrum of the modes and the weight for each mode’s contribution to the two-point functions. As the two-point functions are computed by summing over $k$ modes, taking $R\to\infty$ corresponds to taking the limit of the Riemann sum, yielding the continuous limit. Holding the upper limit $k_{\rm max}$ fixed, we expect spectral convergence with increasing $R$ based on the properties of the basis functions we employ. Unfortunately, the presence of the $\Psi$ pole (see Appendix C) means that as $R$ increases, modes get pushed closer to the pole, eventually leading to divergent values for the two-point functions. Although our trick of shifting $R$ to avoid the pole provides reasonable initializations (see Appendix C), we see consistency in our results as $R$ increases, but not convergence. To properly observe convergence will require resolving the issue of the $\Psi$ pole. As such, we need to find a balance between avoiding the $\Psi$ pole and having sufficient modes both inside and outside the horizon at $t_{0}$.

Nonetheless, we can test for the dependence of various numerical quantities on our selection of $R$, keeping $\lambda=10^{-10}$, $\kappa=5\bar{H}_{0}$, and $k_{\rm max}=4\kappa$ fixed. When we initialize the fluctuations in the Bunch-Davies initial state and parameterize $R=c/\bar{H}_{0}$, we find that quantities such as $\rho_{\rm total}=\bar{\rho}+\delta\rho_{(2)}$ show only modest variation as we vary $c$ between the fiducial value we used throughout our simulations, $c=1.5\pi=4.71$, up through $c=25.85$, as shown in Fig. 18. For initial conditions $(\varphi_{0},\dot{\varphi}_{0})=(25\,M_{\rm pl},-0.25\,M_{\rm pl}^{2})$, we found $N_{\rm total}=63.01\pm 0.78$ as we varied $4.71\leq c\leq 25.85$, indicating minimal dependence of the evolution of the dynamical system on our choice of $R$.

The regulator $\kappa=b\bar{H}_{0}$ is used to constrain the contribution to the two-point functions from small wavelengths. From theoretical analysis, we expect the two-point functions to diverge with increasing $b$, as shown in Section III. The expected divergence is observed numerically. In order to demonstrate results that are immune to changing $\kappa$ requires the implementation of a renormalization scheme, which is beyond the scope of this paper.

Finally, the wavenumber cutoff $k_{\rm max}$ represents the shortest length-scale that we allow to contribute to the two-point functions. So long as $k_{\rm max}$ is chosen to be sufficiently larger than $\kappa$, we observe the expected exponential convergence with increasing $k_{\rm max}$. We chose $k_{\rm max}=4\kappa$, which corresponds to suppressing the initial contributions of the shortest wavelength modes to the various two-point functions by a factor of $\sim 3\times 10^{-4}$.

Acknowledgements

We are grateful to Katy Clough, John T. Giblin, Jr., Alan H. Guth, Tracy Slatyer, Vincent Vennin, and Matias Zaldarriaga for helpful discussions. Portions of this work were conducted in MIT’s Center for Theoretical Physics and supported in part by the U.S. Department of Energy under Contract No. DE-SC0012567. In addition, PF was partially supported by an NSF Graduate Research Fellowship, and KH by MIT’s Undergraduate Research Opportunities Program (UROP).