Starting inflation from inhomogeneous initial conditions with momentum

We investigate the circumstances under which cosmic inflation can arise from very inhomogeneous initial conditions using numerical relativity simulations. Previous studies have not considered cases with non-zero momentum density due to technical challenges with solving the coupled Einstein constraint equations. Here we address these, introducing and comparing several different ways of constructing cosmological initial conditions with inhomogeneous scalar field and time derivative profiles. We evolve such initial conditions with large inhomogeneities in both single- and two-field inflationary models. We study cases where the initial gradient and kinetic energy are much larger than the inflationary energy scale, and black holes can form, as well as cases where the initial scalar potential energy is comparable, as in scenarios where inflation occurs at nearly Planckian densities, finding large-field inflation to be generally robust. We consider examples of initial conditions where a large scalar field velocity towards non-inflationary values would prevent inflation from occurring in the homogeneous case, finding that the addition of large gradients in the scalar field can actually dilute this effect, with the increased expansion and non-vanishing restoring force leading to inflation.


Introduction
Cosmic inflation [1][2][3][4][5] has been proposed as a solution to some of the fine-tuning problems of Standard Big Bang cosmology, namely, the horizon and flatness problems.But in order for inflation to be a successful explanation for the observed homogeneity of the universe on large scales, it should be able to arise from generic, inhomogeneous initial conditions.Inhomogeneities will not necessarily spoil inflation in models where inflation naturally begins at nearly Planckian densities, [5][6][7][8][9], though this has not be studied in detail.The most recent observations of the cosmological microwave background motivate studying of the problem of initial conditions in low energy scale inflation, in which the potential energy density is much below the Planckian scale [8][9][10][11][12][13].Here, we focus on the effects of large inhomogeneities on the onset of inflation, both in scenarios where it occurs at nearly Planckian and sub-Planckian energy scales, using evolutions in fully nonlinear general relativity.
This question has been studied using tools from numerical relativity in a number of papers [14][15][16][17][18][19][20][21][22][23][24][25][26][27][28], complementing work evolving inhomogeneous fields on homogeneous spacetimes [29][30][31][32][33], and using analytic techniques [34][35][36] (see [37] for a short review on the topic).Focusing on more recent work, [23] showed that large field inflation is robust to simple inhomogeneous initial conditions even when the initial gradient energy is many orders of magnitude larger than the vacuum energy density, provided that the universe is initially expanding everywhere and that the scalar field range remains within the slow-roll regime.Reference [23] also included cases with large inhomogeneities that give rise to the formation of black holes, but showed that inflation can succeed even then.This happens because while the overdense scalar field regions collapse into black holes, the underdense regions evolve into voids that become dominated by the inflationary potential energy at later times such that inflation may begin.This line of research was then extended in two ways: [24,27,28,38] expanded the classes of inflationary models under investigation and [24,26,28] the classes of inhomogeneities.In [24], inhomogeneities in both the scalar field profile and the extrinsic curvature were explored.In particular, the initial expansion was assumed to take the following simple form K(⃗ x) = −Cϕ(⃗ x)+K 0 where C > 0 is a free parameter and K 0 is an integration constant and the initial velocity of scalar field was given by some constant such that the momentum constraint is trivially satisfied.In this ansatz, the initial hypersurface contains regions of local expansion and contraction.It was shown that as long as the spacetime is on average initially expanding, then inflation will occur in some patch even if other parts of the spacetime might collapse.Reference [26] investigated inhomogeneities in the transverse traceless part of the extrinsic curvature, A T T ij ̸ = 0, and found that these initial perturbations, roughly corresponding to gravitational radiation, initially reduce the total number of e-folds as the amplitude is increased, but that this reduction saturates and that in general the amplitude of the scalar perturbations remains the main driver of the onset of inflation.The effect of initial inhomogeneities in the scalar kinetic energy, but with constant field profile, such that the scalar momentum density was initially zero, was studied in [28].
In order to carry out such evolutions of various classes of initial conditions to determine whether they eventually lead to an exponentially expanding spacetime, one must begin with initial data that satisfies the constraint part of the Einstein equations.In practice, this requires specifying the values of various metric and matter components on the initial time slice in such a way that they satisfy the Hamiltonian and momentum constraint equations.However, all of the abovementioned studies are similar in that they rely on special choices of initial conditions that ensure the momentum constraint is trivially satisfied.The addition of momentum in the initial conditions is not only natural, but also interesting to pursue since, if the initial velocity of the inflaton were large enough, it could prevent the onset of inflation.Solving the momentum constraint, in addition to the Hamiltonian constraint is highly non-trivial, not only because it requires solving three additional coupled non-linear elliptic equations, but also because it is challenging to separate freely specifiable versus constrained degrees of freedom in a manner where the underlying physical interpretation of the free data is transparent, while also ensuring that a unique solution exists.If we over-restrict the system, no solutions to the constraint equations will exist, also known as the problem of existence.On the other hand, if we do not restrict the system sufficiently, there may be multiple nearby solutions, in which case our numerical solver might not converge to a single one.This socalled problem of local uniqueness is very relevant in cosmology [39], and has implications for how matter source terms in the momentum and Hamiltonian constraint are specified, as discussed below.In addition, a standard trick to choose which variables to set, and which to solve for, while ensuring the solution is unique, is to perform a conformal decomposition of the energy and momentum density.However, this does not extend straightforwardly to scalar fields since the energy and momentum density are functions of the scalar field and its time and spatial derivatives.
In this work, we solve both the Hamiltonian and momentum constraint equations for non-trivial inhomogeneities in the momentum density of the scalar field, considering scenarios where the time derivative of scalar field has a large homogeneous value, as well as scenarios where it has large spatial variations.We solve the equations using the conformal thin-sandwich (CTS) formalism [40], making using of the code described in [41], but without using the usual conformal rescaling of the matter terms.We find that our new scheme gives greater control over the initial conditions in the matter sector.We then study the evolution of several classes of initial conditions that have non-trivial inhomogeneities in the scalar field momentum density and have not been previously studied in the literature.We consider scenarios where the length scale of the inhomogeneities is comparable to the initial effective Hubble radius and (i) the initial scalar gradient and kinetic energy are comparable and much larger than the inflationary energy scale; (ii) the initial scalar gradient, kinetic energy, and potential energy are all comparable (as might arise in large field inflationary models where inflation begins at nearly Planckian densities); and (iii) in two-field inflationary models where the scalar fields are non-interacting.Our results for (i) are, broadly speaking, in agreement with previous analyses that assumed a vanishing initial velocity profile for scalar field.However, we also show that when the initial kinetic energy of scalar field is such that a homogeneous universe will fail to inflate, the addition of gradient energy can slow the scalar field before it reaches the end of the inflationary plateau.We attribute this to an increased initial expansion rate and non-negligible pullback force due to the presence of gradients in the inflaton.Our results suggest that large gradients can mitigate the disruptive effect of a non-zero initial scalar velocity profile.Our results for (ii) show that in cases where the kinetic, gradient, and potential energy are all comparable, the universe will rapidly transition to exponential expansion, typically without forming collapsing regions.Finally, we extend our methods to the study of the effects of adding inhomogeneities to cosmological scenarios where the universe undergoes two stages of inflation.Taken together, our results suggest inflation can arise from highly inhomogeneous conditions.
The remainder of this paper is as follows.We discuss the inflationary models we use in section 2. Our approach to solving the constraint equations on the initial time slice is outlined in section 2.1 and 2.2.We comment on the existence and uniqueness of our solutions in 2.3.Our numerical methods and diagnostics for evolving the Einstein equations and matter are described in section 2.4.Our results are presented in section 3, and we conclude in section 4. In appendix A, we discuss our numerical implementation in more detail, while in appendix B, we provide more details on the construction of our initial data.
We work in four spacetime dimensions, with metric signature (− + ++); we use lowercase Greek letters (µ, ν, ...) to denote spacetime indices and Latin letters (i, j, k, ..., although t is reserved for the time coordinate index) to denote spatial indices.The Riemann tensor is

Theory
In this paper, we consider inflationary theories with one or two canonical scalar fields, ϕ and θ, both minimally coupled to Einstein gravity, such that the action is given by where g is the determinant of the spacetime metric, R the Ricci scalar and V (ϕ, θ) the potential energy function.In most of the cases, we will consider single-field inflation by setting θ = 0.In all cases, we assume that the scalar fields do not interact, such that The Einstein equations that follow from the action (2.1), are then where the stress-energy tensor is and similarly for θ.The equation of motion for each individual field is then where ′ is used to denote the partial derivative with respect to corresponding scalar field.In general, each field will undergo inflation when its potential energy density is greater than the sum of its kinetic and gradient energy density.Assuming the universe is spatially homogeneous, ϕ will undergo inflation if V ϕ > 0 and the slow-roll parameters satisfy: in which case the field is slowly rolling.The condition ϵ V ≪ 1 implies that the potential energy density driving inflation is roughly constant, ρ V ϕ ≡ 2V ϕ ≈ const., while η V ≪ 1 ensures that inflation persists for a number of e-folds.

Initial conditions: metric
We wish to specify initial data on a spacelike hypersurface Σ t parameterized by t that is consistent with the Einstein equations.We formulate the problem using the Arnowitt-Deser-Misner (ADM) formalism, and decompose the metric into a spatial metric γ ij , lapse N , and shift vector β i as and write the extrinsic curvature as where the Lie derivative L is taken with respect to the timelike unit normal to slices of constant coordinate time n µ = 1/N, −β i /N .The initial data for the metric and matter sector cannot be freely and independently specified, but must satisfy the momentum and Hamiltonian constraint equations where R, and D i are, respectively, the extrinsic curvature scalar, Ricci scalar, and covariant derivative associated with γ ij , ρ ≡ n µ n ν T µν is the energy density, and p i ≡ −γ iµ n ν T µν is the momentum density as measured by an Eulerian observer with four-velocity n µ .In the 3+1 decomposition, initial data for the Einstein and matter equations are then a set of 20 functions representing N , β i , γ ij , K ij , ρ, and p i on the initial time slice that together satisfy the constraint equations.
We use a modified version of the CTS method [40], implemented in [41], as a prescription to separate freely specifiable from constrained degrees of freedom in such a way that the underlying physical interpretation of the free data is transparent.We outline the key features below, and refer the reader to [41] for more details.In the CTS formalism, the approach is to perform a conformal decomposition of the spatial metric and extrinsic curvature in order to introduce quantities that can be specified in a more well-motivated way.Introducing the conformal factor Ψ, we can define a number of conformal quantities as where the time-derivative of the conformal metric γij ≡ Ψ 4 γij −1 3 γ ij γ kl γkl is traceless by construction.The conformal lapse Ñ is related to the lapse by Ñ ≡ Ψ −6 N , while R and D are the Ricci scalar and covariant derivative with respect to γij , and L is the corresponding conformal Killing operator, defined by With these definitions, the constraint equations become an elliptic condition for Ψ, and three elliptic equations for β i : In the CTS formalism, one then usually proceeds by specifying initial data for the free data γij , γij , K, Ñ , ρ, p i (2.16) (or some conformally rescaled version of ρ and p i ), and solving the four elliptic equations for the shift β i and conformal factor Ψ. We emphasize that in the CTS method, one is in principle free to choose any values for the free data for which a solution can be found.We now outline and motivate our choice of free data for the case of inhomogeneous conditions in an initially expanding universe filled with two (unless otherwise stated) scalar fields and periodic boundary conditions 1 .
For simplicity, we choose the spatial metric on the initial time slice to be conformally flat, and set the trace of the extrinsic curvature K to be constant, while the transverse-traceless part γij is set to zero (loosely this is equivalent to setting the gravitational wave background to zero).In principle, K is a free parameter, representing a uniform expansion rate across the initial hypersurface.However, we also have to satisfy an integrability condition for the Hamiltonian constraint For a periodic domain, if we approximate the conformal factor to be roughly unity, this requires K2 /3 to be close to the initial energy density averaged over the initial hypersurface.For two scalar fields, this reduces to where ϕ and θ are the values of scalar fields on the initial time slice and ∂ t ϕ and ∂ t θ are specified up to some conformal factor.The minus sign gives us an initially uniformly expanding universe with Hubble parameter H 0 = −K/3 (a positive K would imply a contracting universe).In the special case of a homogeneous scalar field profile, this choice of initial data gives a Friedman-Lemaitre-Robertson-Walker (FLRW) solution.
We now describe the choice of free data in the matter sector.We will restrict the discussion to a single field, say ϕ, but keep in mind that the same definitions and assumptions apply to the other field θ.

Initial conditions: matter
In the case of a scalar field, the quantities to be specified on the three dimensional spatial hypersurface are the scalar field ϕ and its time derivative ∂ t ϕ.We thus rewrite the energy density at some time and any point in the hypersurface as the sum of the kinetic, gradient, and potential energy densities, and the momentum density as where η is the conjugate momentum 2 defined by Clearly, the technical challenge here is to choose which quantities to specify in such a way that not only does a unique solution exist, but also such that the values for the quantities can be physically motivated.The fundamental quantities we freely specify are the initial scalar field profile and a non-trivial velocity profile defined below.We comment on the existence and uniqueness of our solutions in section 2.3.The inhomogeneous initial conditions for the scalar field profile are chosen to be such that on the initial time slice, ϕ(t = 0, ⃗ x) = ϕ 0 + δϕ x sin(kx) + δϕ y sin(ky) + δϕ z sin(kz) (2.22) where ⃗ x = (x, y, z) is the spatial coordinate of hypersurface labelled by the time coordinate t, k = 2π/L is the wavenumber and δ ⃗ ϕ = (δϕ x , δϕ y , δϕ z ) is a measure of the amplitude of the initial inhomogeneities which in general can be different in each direction, but we will only consider δϕ = δϕ x = δϕ y = δϕ z .The maximum total amplitude of the fluctuations For simplicity, we do not consider different amplitudes in the different spatial directions, or additional scalar field variations at smaller wavelengths.In [23], including these was found not to qualitatively affect the results.We fix the (coordinate) length of the simulation domain in each direction to be equal to the wavelength L, and consider various ratios of this lengthscale to the initial Hubble scale (including gradient and all other energy contributions).The initial conditions of the scalar field depend on the amplitude of the inhomogeneities δϕ and the potential V (ϕ 0 ) which sets the background inflationary Hubble scale.
We explore the effects of inhomogeneities on both single-field and two-field inflationary models.When studying single-field inflation, we consider two types of inflationary potentials, both of which are large field models, i.e., the inflaton needs to traverse a super-Planckian range in field space for there to be a sufficient period of inflation.The first potential we consider is a simple quadratic potential, (2.23) Although strongly disfavored by the most recent observations [10,11], we do not expect our conclusions to be dependent on the shape of the potential, since the key feature is the flatness of potential and super-Planckian distance in field space to the minimum of potential (see [27] for a study of the effects of the potential shape on inhomogeneous inflation).To confirm this, we do, however, also look at the T-shaped potential, describing a class of so-called α-attractors, motivated by supergravity and string theory [42,43] and consistent with the most recent BICEP/KECK data [13].The two-field model we will consider is given by where M ≫ m.Here, inflation can start at M 2 ϕ 2 = O(1) if all gradient and kinetic terms are much smaller than V ϕ (ϕ) [9].When the period of inflation driven by ϕ ends, a second stage of inflation driven by θ begins.Our choice of two-field model is motivated by the recently proposed α-attractor generalization [44,45] of the hybrid inflation scenario [46].Although this scenario describes an inflaton field θ, interacting with a Higgs-type field ϕ, one can show that the fields become decoupled in the large field limit, in which case the inflationary model reduces to the two-field model studied here.
Going back to the initial conditions on the velocity profile the scalar field, we find three distinct ways of specifying well motivated initial data, which we now describe.

Uniform initial time derivative of scalar field.
The first method is to specify the time derivative of the scalar field, which for simplicity we assume to be constant on the numerical domain where δϕ ∼ O(1).The sign of the parameter ω determines whether the scalar field moves up or down the inflationary potential.In particular, when ω < 0, the field is driven towards the minimum of the inflationary potential, which could cause an early end to inflation.We will thus give particular attention to this scenario.
We emphasize that, although the time derivative of the scalar field is constant, after the Hamiltonian and momentum constraints are solved for Ψ and β i , the reconstructed physical energy and momentum density are inhomogeneous.While the reconstruction of the physical energy densities is unambiguous, the lack of ability to specify a particular configuration for ρ and its individual components is inconvenient.This motivates another way of specifying initial data.
Here, the trick is to specify a rescaled conjugate momentum, defined as This ensures that if the conformal gradient energy, γij ∂ i ϕ∂ j ϕ and kinetic energy η2 have a chosen ratio, then so will the physically reconstructed quantities following the solution of the constraints.In other words, the ratio of the initial kinetic and gradient energies is not modified by solving for the shift and conformal factor in CTS formalism.
The ansatz for η is chosen to be spatially varying, although with fixed negative sign on the entire domain η = −kδϕ cos 2 (kx) + cos 2 (ky) + cos 2 (kz). (2.28) such that the scalar field initially moves down the inflationary potential.We note that this initial data breaks the symmetry of the energy density about the origin.
Finally, instead of solving the momentum and Hamiltonian constraints using the CTS formalism, a third approach is to, following [47,48], choose initial data such that the momentum constraint is analytically satisfied.We then use the Hamiltonian constraint to solve for the conformal factor.We outline the main assumptions of this construction here, and refer the reader to appendix B for more details.
The initial data for the spatial metric and extrinsic curvature scalar is specified the same way as above, but the lapse and shift are chosen to be N = 1 and β i = 0.The initial data for the scalar field profile is obtained by specifying ϕ(t = 0, ⃗ x) according to (2.30) below, and the conjugate momentum by specifying where it is important to note that η has a different scaling with the conformal factor from η introduced above (2.29) and that the last step assumes N = 1 and β i = 0. Given those assumptions, the momentum constraint is then solved by the following ansatz and a particular solution for Âij given in eqn.(B.4).This allows the reconstructed ∂ t ϕ(t = 0, ⃗ x) to take both positive and negative values in the domain.For this study, we will further assume δ η = δ ηx = δ ηy = δ ηz and δ φ = δ φx = δ φy = δ φz .

Local uniqueness and existence
We now briefly discuss how the different ways of specifying scalar field initial data relate to considerations of how the Hamiltonian and momentum constraint couple, and issues of local uniqueness. 3First, we recall the more familiar case with fluid matter initial data, in order to compare and contrast to the scalar field case.In the former, where the stress tensor is an algebraic function of the fluid variables, one typically specifies a conformal energy density ρ ≡ Ψ n ρ and a conformal momentum density pi ≡ Ψ 10 p i .For the energy density, a value n > 5 is chosen based on uniqueness considerations, which we discuss below, while the exponent of the conformal momentum density is chosen to remove the Ψ dependence in the momentum equation (2.15) [40].In the special case where K is chosen to be constant, the momentum constraint then becomes independent of Ψ, and the constraint equations decouple in the sense that one can first solve the momentum constraint for β i , and then using the solution for β i , solve the Hamiltonian constraint for the conformal factor.In contrast, in the scalar field case matters are complicated by the fact that the energy and momentum depend on the gradient of the scalar field.As can be seen from eq. (2.19), specifying a conformal energy density would involve solving an additional nonlinear partial differential equation to determine the scalar field, and we are typically interested in specifying the initial scalar field profile as the quantity of physical interest anyway.However, noting that pi = −2Ψ 6 ηγ ij ∂ j ϕ is fixed by choosing ϕ (and hence ∂ j ϕ) and say η ≡ ηΨ 6 at t = 0, then this choice of scalar field initial data achieves the same outcome as above, and (since we always choose K to be constant here) the momentum constraint decouples from the Hamiltonian constraint.This is the case considered above when constructing initial data such that the momentum constraints are analytically satisfied (2.29), and is also the case where there are uniqueness results in the mathematical literature for a scalar field on a compact manifold [53].As discussed in Ref. [39], issues with the non-uniqueness of the momentum constraint can arise in this case when the kernel of the operator ∂ j L is non-trivial, which occurs with typical choices of initial data (e.g. a conformally flat metric).These are side-stepped with the analytical solution used here.Fixing ∂ t ϕ, as in eq.(2.29), also removes the dependence of the momentum constraint on Ψ, since Ψ 6 η = ∂ t ϕ + β i ∂ i ϕ / Ñ .However, it introduces an extra term involving β i , so that the previous analyses do not directly apply.For the conformal rescaling (2.27) that fixes the ratio of kinetic to gradient energy, the momentum constraint does not decouple and, again, previous analyses do not apply.
Finally, we consider the local uniqueness of the Hamiltonian constraint.For simplicity, we will ignore the momentum constraint and restrict to the choices for the metric free data discussed in Sec.2.1 (e.g.conformally flat, and so on).Say we have some solution Ψ 0 to the CTS form of the Hamiltonian constraint (2.14), and we are interested in nearby solutions Ψ = Ψ 0 + u where u is small compared to Ψ 0 .Linearizing, we have where, in terms of the functional derivative of the kinetic energy term, If q(⃗ x) > 0, then, using the maximum principle, one can show that u = 0 everywhere, and Ψ 0 is (locally) unique [54,55].For the case of fluid matter mentioned above, the choice of ρ is motivated by ensuring that ρ makes a positive contribution to q.In the scalar case, clearly, this will not be true in general.When fixing η or ∂ t ϕ, the functional derivative term in eq.(2.33) becomes 7Ψ −8 0 η2 /4, and hence makes a positive contribution to q.When specifying η, the contribution is −η 2 /4, and hence negative (as noted above, in this case the momentum constraint does not decouple as well).In either case, for the parameters considered here, we find that q(⃗ x) > 0. In part, this can be attributed to our choice of K through eq.(2.18).A sufficient condition for q > 0 is that when fixing η or ∂ t ϕ, where here ⟨. ..⟩ denotes the coordinate volume average.When fixing η instead, the last term becomes +η 2 /5.In all the cases we consider in this study, we find Ψ 0 ∼ 1, and the gradient and kinetic energy terms, which appear with smaller prefactor on the right-hand side of (2.34) compared to the left-hand side, are approximately equal in magnitude, and either larger, or approximately equal, to the scalar potential energy term.This argument can also be straight-forwardly extend to multiple non-interacting scalar fields, in which case the kinetic, gradient and potential energy densities in (2.34) would simply be the sum of energy densities of individual fields.
Though the above discussion provides some motivation and guidance, in the end there is no mathematical result guaranteeing the existence and uniqueness of solutions to the constraint equations for all the cases we consider.However, the important thing for the purposes of this study is that we find we are able to numerically construct convergent solutions.In fact, there are several examples in the literature where the constraint equations have been solved without issue when it was known that more than one solution exists [56][57][58][59][60].

Numerical implementation
We construct initial data satisfying the constraint equations by numerically solving the CTS equations, discretized with second-order accurate finite differences, using a multigrid algorithm.Details can be found in [41].We then evolve the Einstein equations in the generalized harmonic formulation using the code described in [61,62], and assuming periodic boundary conditions.See appendix A for more details on our numerical methods.In order to characterize our results we make use of the following diagnostic quantities.We compute the total and individual energy densities defined according to (2.19).We define a fiducial local Hubble expansion rate which allows us to compute the corresponding local number of e-folds of expansion N as where the integration is along the integral curve of n µ , and τ is the proper time given by the lapse through dτ = N dt (see [23,48]).Note that in a homogeneous (i.e.FLRW) spacetime this quantity is related to the initial and final scale factor as N = Hdτ = log(a/a 0 ).The average of a quantity over the spatial volume is denoted by where γ = det γ ij .
In some cases (in particular, for smaller values of k/H 0 ), there is prompt gravitational collapse in some regions and black holes form.To characterize the boundaries of black holes in our dynamical setting, we use the concept of an apparent horizon, defined as the outermost marginally outer trapped surface.We excise a region interior to the black hole from the computational domain, and do not evolve the equations there (see appendix A for more details on the implementation).From the area of the apparent horizon A we define the irreducible mass, M irr = A/16π.We find that the angular momentum of the black holes is negligible, and hence the irreducible mass closely approximates the Christodoulou mass.To make contact with observations, we compute the power spectrum of scalar fluctuations, which restoring Planckian units and assuming slow-roll inflation is given by evaluated at the time when the cosmic microwave background fluctuations crossed the horizon.
We also characterize the parameters we choose for the background cosmological solutions (to which we add large inhomogeneities) by quoting the scalar spectral index and the tensor-toscalar ratio, which in the slow-roll approximation are given by (2.39)

Results
We examine the conditions under which an initially expanding universe, with either a single scalar field, or two non-interacting scalar fields (see section 2), transitions to exponential expansion.We construct initial data using the approach described in sections 2.1-2.2.
We study two regimes of interest, distinguished by the energy scale at which inflation takes place.We first consider single-field inflationary models where the universe is initially dominated by the gradient and kinetic energy of scalar field, η 2 ∼ γ ij ∂ i ϕ∂ j ϕ ∼ 10 3 V (ϕ 0 ), in which case the initial expansion rate is large compared to the inflationary Hubble rate (we refer to this as "low energy" or sub-Planckian inflation, see section 3.1).We find that one or more (depending on the symmetry of the total energy density on the initial time slice) regions of the domain undergo gravitational collapse, leading to the formation of black holes.Ignoring the interior of the black hole(s), the spacetime is still expanding on average, and we find that, independently of the type of initial data, the gradient and kinetic energy dilute away until the numerical domain becomes dominated by inflationary energy.
We also explore inflationary models where the amount of kinetic, gradient, and potential energy are all comparable on the initial time slice, and their sum is order unity in units of the wavelength of the initial homogeneities (which is also comparable to the initial Hubble radius) ρL 2 ≲ O(1) (see section 3.2).Here, we refer to this scenario as high energy, or nearly Planckian energy inflation, though, of course, our analysis is entirely classical and ignores quantum effects.We first consider single-field inflationary models, where we find that the solution quickly transitions to exponential expansion.Independently of the way initial data is constructed, the late-time evolution always consists of an exponentially expanding universe with a rapidly decreasing gradient energy and kinetic energy that asymptotes to its corresponding slow-roll value.We then extend our study to a two-field inflationary model where the energy content is such that the first stage of inflation may start at nearly Planckian energies, while the second phase will occur at sub-Planckian energies.For technical reasons discussed below, we are not able to evolve this scenario past the first stage of inflation.However, already at this stage we find large patches where gradients have become negligible and the spacetime is well described by a homogeneous evolution set by the local inflaton values.Extrapolating thusly indicates that these regions should undergo a prolonged period of inflation driven by the second field.
In the following, we quantify our results using the methods described in section 2.4, by computing, for instance, the ratio of kinetic and gradient energy to the potential energy.In all the cases we study, we consider perturbations such that the ratio of the wavenumber of the initial scalar field content to the expansion rate on the initial time slice k/H 0 is of the order of one, meaning our simulations probe the strong-field regime.As here we focus on models of large-field inflation, where ϕ 0 is large in Planckian units, and consider fluctuations of the the inflaton that are of order one, we note that we do not consider scenarios where the perturbations around the average scalar field value that exceed the distance in field space to the end of inflationary plateau.Previous studies [23,24,27] have found that inflation can fail to happen in such scenarios.

Low energy inflation
We first consider solutions where the potential energy density is initially subdominant.To begin with, we describe an illustrative case where the time derivative of the scalar field is uniform on the initial time slice, i.e. given by eq.(2.26).We choose the potential to be given by an α-attractor potential eq.(2.24) with ϕ 0 = 3.1, such that, in the absence of inhomogeneities (δϕ = 0), and if the time derivative of the scalar field were zero, the combination of parameters (ϕ 0 , m, α) would result in 60 e-folds of inflation, a scalar power spectrum of ∆ 2 R ∼ O(1), a scalar spectral index of n s = 0.97, and a tensor-to-scalar ratio of r = 0.001 for modes that cross the horizon 60 e-folds before the end of inflation.However, we pick an initial value of ∂ t ϕ such that, in a homogeneous Universe, the scalar field would evolve off the plateau of the inflationary potential without inflation occurring.The gradient and kinetic energy are chosen to be comparable on the initial time slice, but 800 times larger than the potential energy density.Figure 1 shows the individual volume-averaged energy densities (left panel) as a function of the effective scale factor a = exp(⟨N ⟩).We also plot the kinetic φ2 FLRW and potential energy ρ V ϕ,FLRW densities when solving the homogeneous FLRW equations in the absence of gradient energy and specifying the initial value of φFLRW to be the time derivative of the scalar field on the initial time slice given by eq.(2.26).
We find that in some regions, the maximum energy density quickly increases and black holes form (at the times indicated by the grey shaded areas in figure 1).As was pointed out in previous studies [23,24], the formation of black holes can be motivated using the hoop conjecture [63].Indeed, the hoop conjecture predicts that, if the mass of an overdensity exceeds the mass of a black hole of the same size, then the overdensity will collapse to a black hole.Following this argument, one can expect black holes to form when 4  3 πk −3 ρ ∼ k −1 /2 which is equivalent to k/H 0 ∼ 1, or δϕ ∼ 1 in Planck units.Here, since the total energy density on the initial time slice is symmetric with respect to positive and negative values of the scalar field, two identical black holes form in the periodic domain.These black holes only create locally collapsing regions, and, after removing their interior from the domain of integration, the spacetime is still expanding on average.The gradient and kinetic energy are diluted until the spacetime is dominated by the potential energy, after which inflation starts.
As the effective scale factor of the spacetime increases, the proper distance between the black holes also increases, as their density is also diluted by inflation.Figure 1.Left: The volume-averaged kinetic η 2 , gradient ρ grad , and potential ρ V ϕ energy densities plotted against the averaged measure of the effective scale factor for the case where the time derivative of the scalar field is uniform on the initial time slice and the potential has a T-shape (2.24) with (ϕ 0 , m, α, ω, δϕ) = (3.0916,0.14, 1/6, −7.1, 0.9).For comparison, we also show the evolution of the kinetic energy density φ2 FLRW , and potential energy density ρ V,FLRW when solving the corresponding homogeneous FLRW equations in the absence of gradient energy and by specifying the initial value of φFLRW to be the time derivative of the scalar field on the initial time slice.The energy density for a radiation dominated universe is shown by the dotted yellow line.Middle: The volume average of inflaton field value in comparison to the homogeneous solution.Right: The value of scalar field at the locations ⃗ x 1 = L(0, −1/4, 1/4) and ⃗ x 2 = L(1/8, −1/4, 1/4) when solving full equations of motion, and when solving equation for damped harmonic oscillator (3.3).
Surprisingly, we find that the addition of gradient energy allows inflation to occur in cases where it would not in the equivalent homogeneous case.We will now demonstrate that the reasons are two-fold.On the one hand, there is some resistance from the gradient pressure, which acts as a restoring force, and hence tends to pull the scalar field up the potential and away from the minimum.On the other hand, the addition of gradient energy results in an increase in the initial expansion rate, as a result of which the oscillations of the inflaton are damped.We now confirm this picture, and propose a simple toy-model which allows us to reproduce our results and gain some intuition.Let us consider the Klein-Gordon equation on the background of an FLRW spacetime, where we have included the gradient terms in addition to the usual friction term, and the dot here denotes the time derivative with respect to the proper time τ .We further assume that the gradient energy can be approximated by a homogeneous energy density initially equal to the volume average, and that scales with an inverse power p of the scale factor where V 0 ≡ V ϕ (ϕ = ϕ 0 ) is the initial potential in the absence of inhomogeneities, and r gives the ratio of the initial value of gradient energy density to the inflationary energy density r ≡ ⟨ρ grad,0 ⟩ 2V 0 .In addition, we assume that the gradient energy redshifts in the same way as radiation, i.e. p = 4.In order to treat the spatial derivative term, we make the simplifying approximation that ∂ 2 ϕ ∼ (ik/a) 2 ∆, where a is the spatially homogeneous scale factor, and we have introduced the deviation from the background solution ∆ ≡ ϕ(t, x, y, z) − χ 0 (t), where χ 0 (t) is the solution to the equation of motion (3.1) assuming the third term vanishes.Writing this in terms of the homogeneous number of e-folds, dN = Hdτ , the equation of motion (3.1) then becomes where, since H 2 = ρ/3 with ρ given by (2.19), using (3.2), and further assuming η 2 ∼ ρ grad , the coefficients are and where H 2 V = 2V 0 /3 is the initial inflationary Hubble rate.This expression takes the form of a damped harmonic oscillator, where H 2 is analogous to the mass, A represents the Hubble friction, B the restoring pullback force, and the potential gradient the driving force pushing the inflaton down the potential.It is now clear that, for r > 0, the additional gradient energy not only increases the initial friction as well as the mass term, but also provides an additional pullback force for points away from the average value of the inflaton.Note that, as the universe expands, the contribution from the gradient energy to the friction and restoring terms is exponentially suppressed, such that if the scalar field remains within the inflationary part of the potential during the initial dynamical phase, then it can transition into slow-roll inflation.
In the right panel of figure 1, we compare the results of numerically integrating the damped harmonic equation (3.3) to the behavior of the full simulations at two arbitrarily located points.We find that the results are qualitatively the same and conclude that this simple toy-model can be used to build intuition about the dynamics of the full solution.Similarly, the left panel shows that the gradient and kinetic energy of the full evolution decreases as a −4 , as would be expected for a radiation dominated universe.We also note that the arguments given above suggest that other significant contributions to the expansion rate from inhomogeneities, for example from a gravitational wave energy density, as studied in [24], should also mitigate the disruptive effect of a small nonzero initial scalar field velocity.
We find that qualitatively similar results hold when starting from initial conditions constructed by specifying a rescaled version of the conjugate momentum such that η = Ψ 2 η (2.28), or such that the momentum constraint constraint is analytically satisfied η = Ψ 6 η (2.30).In both cases, the rescaled conjugate momentum is chosen such that the gradient and kinetic energy are comparable on the initial time slice, but ∼ 1000 times larger than the potential energy density.The former is demonstrated in the right panel of figure 2, where we choose the potential given by an α-attractor potential (2.24), and an initial scalar field profile given by (2.28), with parameters such that, in the absence of inhomogeneities (δϕ = 0), and if the conjugate momentum of the scalar field were zero, the combination of parameters (ϕ 0 , m, α) would result in 60 e-folds of inflation, a scalar power spectrum of ∆ 2 R ∼ 0.05, a scalar spectral index of n s = 0.97, and tensor-to-scalar ratio of r = 0.006 for modes that cross the horizon 60 e-folds before the end of inflation.We find that the gradient and kinetic energy decrease until the spacetime is dominated by the inflationary energy after which inflation may start.Unlike for initial data with a constant velocity profile, here the energy density on the initial time slice is no longer symmetric with respect to positive and negative values of the scalar field, and therefore the black holes are no longer symmetric about the origin.Similarly, the left panel of figure 2 shows the individual energy densities for a solution constructed according to (2.30) and with quadratic potential such that in the absence of inhomogeneities and conjugate momentum, the combination of parameters (ϕ 0 , m) would result in 60 e-folds of inflation, ∆ 2 R ∼ 1 × 10 −3 , n s = 0.97, and r = 0.13 for modes that cross the horizon 60 e-folds before the end of inflation.Here, again, two regions collapse to black holes, but the universe keeps expanding, and eventually transitions to exponential expansion.These results emphasize that our conclusions are generic, and not sensitive to the particulars of how the scalar field momentum is chosen.

Single-field inflation
We next consider single-field inflationary models where inflation may naturally begin near the Planck scale.Thus, solutions where the sum of the kinetic, gradient and potential energy of the inflaton is order unity in Planck units.We further restrict our study to initial conditions where the kinetic, gradient, and potential energy are initially comparable, and the potential is quadratic in the scalar field, V ϕ (ϕ) = m 2 ϕ 2 .In figure 3, we show the volume-averaged total energy density (left panel) and the volume-averaged expansion rate (right panel) computed from (2.35), as a function of the effective scale factor, for the three types of initial data described in section 2.2.We specify the initial data such that, in the absence of gradient and kinetic energy, the homogeneous initial value of ϕ 0 = 11 and other parameters would result in 60 e-folds of inflation, ∆ 2 R ∼ O(1), n s = 0.97, and r = 0.13 for modes that cross the horizon 60 e-folds before the end of inflation.Independent of the type of initial data used, as the scale factor increases, the average expansion rate and total energy density decreases (although at different rates depending on the initial data), until the kinetic and gradient energy is subdominant, leading to a universe dominated by the inflationary potential, hence undergoing exponential expansion.This is illustrated in figure 4, where we show the individual average energy densities.Although the gradient energy quickly becomes negligible, the kinetic energy asymptotes to a constant value, consistent with the value the kinetic energy would asymptote to when solving the homogeneous FLRW equations in the absence of gradient energy and specifying the initial value of φ2 FLRW to be the volume average value of the kinetic energy of inhomogeneous solution on the initial time slice.
Though we evolve for long enough that the spacetime is dominated by the inflation potential energy and exponential expansion, we do not evolve to the end of inflation, as growing gradients from different regions in the domain inflating at different rates eventually lead to large truncation error.
We confirm that inhomogeneities do not prevent inflation in models of chaotic inflation occurring near the Planck scale [6][7][8].Although the simplest models of inflation (including the quadratic potential we study here), are now strongly disfavored by the most recent observations from BICEP/Keck and Planck [10,11], they could still appear in the context of more complicated inflationary models containing more than one scalar field.We discuss an example of such model in the following section.

Two-field inflation
Up to this point, we have only considered single field models.However, our methods can be applied equally well to more complicated inflationary models.To give an example, we consider a simple theory describing two non-interacting scalar fields, θ and ϕ, such that the potential is given by (2.25), which combines a quadratic dependence for ϕ with an α-attractor term for θ.We further assume that M ≫ m.Here, the potential can give rise to two stages of inflation.The first stage of inflation is driven by the quadratic potential, and may start at nearly Planckian energies and be relatively short.Once inflation driven by ϕ ends, the second stage of inflation driven by the α-attractor potential, and lasting more than 60 e-folds, may begin.Here we are not able to evolve this model past the first stage of inflation, likely due to limited numerical resolution and the gauge choice.However, already at this stage we find that gradient energy has become subdominant, and by comparing the evolution in different regions to the homogeneous evolution with the same scalar field values, we can approximately extrapolate their later behavior.
We consider values for ϕ 0 and θ 0 such that, in the absence of inhomogeneities (δϕ = 0 and δθ = 0) and kinetic energy on the initial time slice, the first stage of inflation would last one e-fold, and the second stage of inflation would last 60 e-folds and give rise to a scalar power spectrum with ∆ 2 R ∼ 0.1, n s = 0.97, and r = 0.006.The remaining free data is chosen such that the total energy density is initially on the order of 1/L 2 , the potential energy is The volume-averaged total energy density (left) and expansion rate (right) plotted against a volume-averaged measure of the effective scale factor for cases with a quadratic potential given by (2.23).The quantities ⟨ρ V,ϕ ⟩ and ⟨H V ϕ ⟩ respectively represent the volumeaveraged measures of the potential energy density and inflationary expansion rate as a function of time.Both quantities can be seen to asymptote to the values expected for potential energy dominated inflation.The solid line represents a solution for initial data constructed using (2.26) and (ϕ 0 , m, ω, δϕ) = (11.0,0.23, −5.83, 0.79).The dashed line represents initial data given by (2.27) with (ϕ 0 , m, δϕ) = (11.0,0.13, 0.85).Finally, the dash-dotted curve represents initial data given by (2.29) and ( φ0 , m, σ, δ φ, δ η) = ( √ 2 11.0, 0.12, 0.1, 0.68, 2.21).
The volume-averaged kinetic η 2 , gradient ρ grad and potential ρ V ϕ energy densities plotted against the averaged measure of the effective scale factor for the cases shown in figure 3.For comparison, the evolution of the kinetic energy density when solving the FLRW equations in the absence of gradient energy and specifying the initial value of φ2 FLRW to be the volume average value of the kinetic energy of inhomogeneous solution on the initial time slice, is also shown (dashed yellow line).
1/7 of the total energy, and the remaining fraction is provided by the kinetic and gradient energy of the fields θ and ϕ, with each term, for each field contributing equally to the total energy content.We choose an initially uniform time derivative of the scalar field for ease of implementation, but based on the previous section we do not expect our results to depend on the initial velocity profile of the scalar fields.The conjugate momentum on the initial time slice is such that, if one were to solve the FLRW equations in the absence of gradient energy, but specifying the initial values for ( φFLRW , θFLRW ) to be the volume average of conjugate momenta of the corresponding scalar fields on the initial time slice, then the first stage of inflation would only last half an e-fold and the second stage 35.The left panel of figure 5 shows the individual contributions to the total energy density for each scalar field in our simulation compared to when solving the FLRW equations in the absence of gradient energy and by specifying the initial value for the time derivative of scalar field to be the volume average value of the conjugate momentum of inhomogeneous solution on the initial time slice.
We first study the first stage of inflation driven by the quadratic potential, shown in the left panel of figure 5, in detail.In line with our earlier observations, we find that the average kinetic energy of the scalar field ⟨η 2 ϕ ⟩ (dashed red line) does not decrease as quickly as in the homogeneous case (dashed orange line).We also note that the average potential energy of ϕ (dash-dotted green line) increases at first, as does the average scalar field ⟨ϕ⟩ initially, shown in figure 6, likely pulled up by the gradients in the field.As the universe expands, the gradient energy decreases and eventually, when the gradient energy is small enough, the potential energy starts decreasing.
. The gradient ρ grad , kinetic η 2 , and potential energy ρ V ϕ contributions to the volume average energy density as a function of the effective scale factor for each individual scalar field (ϕ left panel and θ middle panel) and both scalar fields combined (right panel).Left: Individual contributions to average energy density for scalar field with quadratic potential driving first stage of inflation.For comparison, the evolution of the kinetic energy density η 2 FLRW and potential energy density ρ V,FLRW when solving the FLRW equations in the absence of gradient energy, and specifying the initial value for the time derivative of scalar field to be the volume average value of the conjugate momentum of inhomogeneous solution on the initial time slice, is also shown.Middle: Same as the left panel, but for a scalar field with the α-attractor potential driving second stage of inflation.Right: The sum of the volume-averaged gradient and kinetic energy densities for each scalar field separately ⟨ρ grad + η 2 ⟩ and the volume-averaged potential energy densities for each scalar field.
Unfortunately, we were not able to able to evolve this scenario past the first stage of inflation.We attribute this to a lack of resolution to resolve large gradients developing between different regions of the domain inflating at different rates.We next illustrate these different regions and conjecture that if we had sufficiently large resolution, then the effectively causally disconnected Hubble patches will keep inflating, albeit at different rates, until the second stage of inflation kicks in.
Figure 6 shows the range of values the scalar field driving the first stage of inflation and the conjugate momentum take, as a function of the effective scale factor.This already suggests that different patches will inflate at different rates, which in turn will lead to the formation of gradients in the field values across the numerical domain.
In particular, the region where sharp gradients develop in the field values, which we eventually can no longer resolve due to the finite numerical resolution of our grid, corresponds to where the scalar field ϕ stops inflating first.Figure 7 shows the scalar field ϕ, its conjugate momentum η ϕ , the extrinsic curvature K, and the gradient energy ρ grad at the location where these unresolved features develop (indicated by φ0 ), but shortly before these become severe enough that the evolution must be halted (indicated by vertical dashed grey line in figure 9).We find that the extrinsic curvature is larger than in the rest of the domain, and that the gradient energy is negligible, as expected.We also show the value of the scalar field θ and its conjugate momentum η θ for the field driving the second stage of inflation.These quantities indicate that the scalar field driving the second stage of inflation is still inflating at that location.
As the universe expands, the left and middle panels of figure 5 show that the average gradient energies of ϕ and θ are quickly diluted until they are negligible compared to the average inflationary energy of either field.This suggests that when the gradient energy becomes negligible in a neighbourhood surrounding a point of interest, then one can use the FLRW equations together with the local value of scalar field and its time derivative to compute the evolution of scalar field at that location.We first show that this approximation is valid for the location where unresolved gradients eventually develop.We then use this to extrapolate the behavior of both scalar fields, at different locations, past the point we are able to evolve.
In the left panel of figure 9, we plot the values of the respective scalar fields at the location where unresolved features eventually develop ( φ0 , θ0 ) as a function of the local scale factor (purple dots).The smallest value of the scale factor for which we show the values of the scalar fields corresponds to the time at which the gradient energy becomes negligible in a neighbourhood surrounding the point of interest.We also show the evolution of the scalar field when solving the FLRW equations (dash-dotted lines), ignoring the gradient terms, and fixing the initial values for the scalar field and its time derivative to be the corresponding values for the inhomogeneous scalar fields and conjugate momenta at that spatial point.As expected, we find close agreement between the numerical and FLRW solution, and make use of the FLRW solution to extrapolate the behavior of scalar fields past the point we were able to numerically evolve.For the spatial point where ϕ stops inflating first, we find that θ will still inflate for ∼ 250 e-folds at that particular point.
As another representative example, we also show the scalar field, its conjugate momentum, the extrinsic curvature, and the gradient energy at the location where ϕ obtains its maximum value (indicated by φmax in figure 8).Similarly, we find that the gradient energy is negligible in a neighbourhood of that point and plot the scalar field values and their corresponding FLRW solutions as a function of the local scale factor in the right panel of figure 9.
We extrapolate the FLRW solution, and deduce that ϕ and θ will inflate for an additional 3 and 250 e-folds respectively.
These results suggest that other parts of the domain will also keep inflating, although at different rates determined by the local value of the scalar field and its conjugate momentum at that location.We thus conjecture that, if we were able to continue the evolution, ϕ would fluctuate around zero in regions where it has reached the bottom of the potential (just like in the homogeneous case).We leave a more detailed analysis to confirm this to to future work.
Assuming the above picture, as the universe keeps expanding, the potential energy of the α-attractor potential eventually dominates the quadratic potential, after which the second stage of inflation starts.The middle panel of figure 5 shows that the gradient and kinetic energy of the scalar field driving the second stage of inflation becomes smaller than the potential energy even before the α-attractor potential dominates over the quadratic potential.Furthermore, we find that (shortly before large gradients make the evolution inaccurate) the wavelength of the perturbations ∼ L is about 60 times larger than the volume averaged Hubble radius of scalar field driving second stage of inflation.This suggests that by the time the first stage of inflation ends, the universe is effectively made up of patches that are homogeneous on Hubble scales and, at least in most regions, the second stage of inflation may generate a nearly-scale invariant spectrum consistent with observations.For comparison, we also show the evolution of ϕ FLRW , the scalar field when solving the FLRW equations in the absence of gradient energy, and specifying the initial value for the time derivative of scalar field to be the volume average value of the conjugate momentum of inhomogeneous solution on the initial time slice.Right: Similar quantities, but for the conjugate momentum of the scalar.

Discussion and conclusion
In this paper, we have undertaken the first study of the robustness of large field models of inflation to large initial gradients and kinetic terms in the initial scalar field profile.We proposed three different ways of specifying scalar field configurations with non-vanishing momentum density when solving the Hamiltonian and momentum constraint equations using the CTS formalism.Solving the full system of coupled constraint equations allowed us to examine the effects of several initial conditions not previously studied in the literature.Our results are summarized as follows:  • We studied the evolution of single-field inflationary models where the initial gradient and kinetic energy densities of the inflaton are much larger than potential energy.We demonstrated that the presence of significant gradient energy can actually negate the effect of a large time derivative in the scalar field, allowing inflation to occur in cases where, in the homogeneous case, the kinetic energy would drive the scalar field off the inflationary plateau before exponential expansion occurred.We explained this in terms of the increased initial expansion rate and restoring pullback force in the presence of large gradients, deriving a simple toy-model to demonstrate this.This implies that rather natural initial gradients in the initial conditions will mitigate the effects of a non-zero initial field velocity of the inflaton.
• We also performed the first study of the impact of inhomogeneities on single-field inflationary models where inflation may start at nearly Planckian densities.We found that, in cases where the gradient, kinetic and potential energy densities are all comparable at nearly Planckian densities, the universe will rapidly transition to exponential expansion.In these cases we did not find black hole formation.
• In addition to single-field inflationary models, we performed a preliminary study of the effects of inhomogeneities on cosmological scenarios where the universe undergoes two stages of inflation, the first one at nearly Planckian energies and the second one at sub-Planckian scales consistent with observations.Although we were not able to evolve the spacetime to the point where the second stage of inflation would start, and there-fore make detailed conclusions, we find at this time large regions where the gradients are negligible and the spacetime is well approximated by the homogeneous evolution determined by the local inflaton values.We leave a more detailed study of these models for future work.
Exploring the two-field inflationary models considered here would require developing new numerical methods.As mentioned above, we believe one of the main difficulties lies in being able to resolve large gradients developing between different regions of the numerical domain inflating at different rates.A better choice of gauge, as well as working with higher resolution or adaptive mesh refinement may help address this, as well as allow us to evolve the single-field inflationary models for longer.
In this study, we only considered two families of inflationary potentials out of the multitude that have proposed in the literature.However, we verified that, starting from the same initial conditions, these two different potentials gave similar qualitative results and we do not expect our results to depend strongly on the details of the potential, since the key features are the flatness of the potential and the super-Planckian distance in field space that must be traversed to reach the minimum of the potential.For two-field inflationary models, for simplicity in this first study, we focused on theories where the scalar fields are non-interacting.It would be interesting to consider the evolution of two-field inflationary models where the scalar fields do interact and explore whether there are still robust to generic initial conditions (see e.g.[8] and [38]).Note however, that in the α-attractor generalization [44,45] of the hybrid inflation scenario [46], the fields become decoupled in the large field limit, and the evolution of the initial conditions reduces to the one performed in our paper.
Another direction for future work would be to adapt the methods developed here to study the robustness inflation in the presence of other types of matter.One particular class of models of interest would be matter where the speed of sound of the fluctuations, or equivalently the Jeans length, goes to zero [64], as in, for example, dust.This would be a case especially prone to the gravitational collapse of initial perturbations.
The methods developed here for constructing general relativistic initial data are rather generic, and could be applied to a number of other cosmological scenarios where it is important to include momentum density.It would also be interesting to compare these methods to those proposed in the recent work of [65], where the authors construct initial data with non-trivial matter configurations by using a new scheme based on the Conformal Transverse-Traceless (CTT) formalism to solve for both constraint equations.In this so-called CTTK approach, rather than choosing a constant value for the trace of the extrinsic curvature K, one chooses an initial profile for the conformal factor, and solves an algebraic equation for a spatially varying K.
(www.scinethpc.ca) and Digital Research Alliance of Canada (alliancecan.ca).Calculations were performed on the Symmetry cluster at Perimeter Institute, the Niagara cluster at the University of Toronto, and the Narval cluster at Ecole de technologie supérieure in Montreal.

A Numerical methodology
We solve the equations of motion (2.2) and (2.5) using the generalized harmonic formulation, as described in [61].The numerical scheme we use follows that of [62], which we briefly summarize here.We discretize the partial differential equations in space, using standard fourth-order finite difference stencils, and in time, using fourth-order Runge-Kutta integration.We control high frequency numerical noise using Kreiss-Oliger dissipation [66].We use constraint damping to control the constraint violating modes sourced by truncation error, with damping parameter values similar to those used in black hole formation using the generalized harmonic formulation [61].We fix the gauge freedom through specifying the source functions H α , choosing damped harmonic coordinates [67,68], as in [23,69].During the expansion phase, we dynamically adjust the time step size in proportion to the decreasing global minimum of 1/N where N is the lapse (this would be N = 1/a 3 in a homogeneous FLRW universe with harmonic time) in order to avoid violating the Courant-Friedrichs-Lewy condition [23,70,71].
Following [61], we track the evolution of any apparent horizons by finding the surfaces where the outer null expansion vanishes, and excise an ellipsoid-shaped region interior to the horizon.We typically set the ratio of the maximum ellipsoid axis to the maximum black hole radial value to be 0.78.
The simulation domain is three dimensional with length L in each direction.The simulations are performed with between 384 and 512 points across each linear dimension.We construct initial data describing an inhomogeneous cosmology as described in section 2.2 and 2.1 using the the conformal thin sandwich formalism, as described in [41].The elliptic equations are solved using a second-order accurate multigrid scheme with a typical resolution of 384 points across each linear dimension.
Finally, we present a convergence test of our code and setup.In figure 10, we show the time evolution of the norm of the constraint violation C α ≡ □x α −H α for the strong-field case considered in figure 2 for different numerical resolutions.For this case, the lowest resolution is 256 across each linear dimension for the evolution and 192 for the initial data code.The medium and high resolutions correspond, respectively, to an increased resolution of 3/2 and 2× that of the lowest resolution, both in the initial data and evolution code.We find that the constraints initially converge to zero at roughly second order and then eventually transition to fourth order convergence.This is because the truncation error is dominated by the second order convergence of the initial data code at first, but eventually the truncation error of fourth order evolution code takes over.The high resolution in the convergence study is equivalent to the resolution we use for most of the cases studied here.

B Modified CTS formalism
In some cases studied here, instead of solving both the momentum and Hamiltonian constraint using the the CTS formalism [41], we construct initial data by following a similar approach to [47,48], generalized to three dimensions, and choosing initial data such that the momentum constraint is automatically satisfied.We then solve for the conformal factor where f 0 , f 1 , e 0 , e 1 , d 0 , d 1 , and σ are parameters to choose and Â11 (x, y, z) = −η 0 f 1 cos(kx) − The conditions on the constants f 0 , f 1 , e 0 , e 1 , d 0 , and d 1 come from substituting the ansatz into (B.1).

Figure 6 .
Figure 6.Left: The minimum ϕ min , maximum ϕ , and volume-averaged values of the scalar ⟨ϕ⟩ driving the first stage of inflation as a function of effective scale factor for the case shown in figure5.For comparison, we also show the evolution of ϕ FLRW , the scalar field when solving the FLRW equations in the absence of gradient energy, and specifying the initial value for the time derivative of scalar field to be the volume average value of the conjugate momentum of inhomogeneous solution on the initial time slice.Right: Similar quantities, but for the conjugate momentum of the scalar.

Figure 7 .Figure 8 .
Figure 7. Two-dimensional spatial slices showing different quantities shortly before strong gradients (around the point indicated by the arrow) make the evolution inaccurate (at the time indicated by vertical dashed grey line in figure 9) for the case shown in figure 5.The first column shows the scalar field (top) driving the first stage of inflation and its conjugate momentum (bottom).The second column shows the corresponding quantities for the scalar driving the second stage of inflation.The third column shows the extrinsic curvature scalar (top) and total gradient energy.

Figure 9 .
Figure 9. Scalar field and conjugate momentum values as a function of the local effective scale factor for two spatial points (indicated by φ0 and φmax on the spatial slices in figures 7 and 8, respectively), chosen such that gradient energy is negligible in surrounding neighbourhood (blue dots).The corresponding homogeneous FLRW solution using the same initial values plotted in the figure is also shown.The vertical dashed grey lines indicate the local scale factors corresponding to the two-dimensional spatial slices shown in figures 7 and 8.