Beginning inflation in an inhomogeneous universe

Using numerical solutions of the full Einstein field equations coupled to a scalar inflaton field in 3+1 dimensions, we study the conditions under which a universe that is initially expanding, highly inhomogeneous and dominated by gradient energy can transition to an inflationary period. If the initial scalar field variations are contained within a sufficiently flat region of the inflaton potential, and the universe is spatially flat or open on average, inflation will occur following the dilution of the gradient and kinetic energy due to expansion. This is the case even when the scale of the inhomogeneities is comparable to the initial Hubble length, and overdense regions collapse and form black holes, because underdense regions continue expanding, allowing inflation to eventually begin. This establishes that inflation can arise from highly inhomogeneous initial conditions and solve the horizon and flatness problems, at least as long as the variations in the scalar field do not include values that exceed the inflationary plateau.


Introduction
Inflation was originally proposed to account for the high degree of flatness and homogeneity observed today in our universe at horizon scales [1][2][3][4]. However, an oft-repeated criticism is that, while it is true that if inflation starts, it rapidly inflates away any inhomogeneities, in order to begin inflation requires homogeneity over several or many Hubble volumes, because otherwise the expanding universe will be increasingly dominated by inhomogeneities rather than inflationary potential energy. Hence, how can inflation explain homogeneity on Hubble scales today?
Here we address to what degree inflation requires initial homogeneity. Arguably this problem does not appear if inflation can start at nearly Planckian density [4,5]. There are also several methods for solving the problem of initial conditions for low energy scale inflation, see [6][7][8][9][10][11] and references therein. There has also been work establishing "no hair" type results governing the decay of perturbations for spacetimes with a cosmological constant, including important work showing how to construct a large class of solutions by imposing the requirement that the solution approach de Sitter at late times [12,13]. A review of this can be found in [14]. Here we concentrate on single-field inflation models where the energy scale of inflation is orders of magnitude smaller than the Planck energy, which have become more popular in the light of observational data [15]. Since small perturbations certainly do not prevent inflation from beginning, the cases of interest are those with large initial inhomogeneities, outside the regime of linear theory and requiring numerical solutions.
This question has been studied in the past using general relativistic simulations in 1D [16,17] and spherical symmetry [18] (see also [19] for early work in 3D), as well as in the multi-field inflation case but ignoring inhomogeneities in the gravity sector [20]. In this work, we solve the full, nonlinear Einstein equations numerically including cases with large inhomogeneities that give rise to collapsing regions and trapped surfaces. We find that many cosmologies that are initially expanding, highly inhomogeneous and strongly dominated by spatial gradient energy give rise to inflation. This happens because while overdense regions collapse into black holes, the underdense regions evolve into voids that become dominated by inflationary vacuum energy at (much) later times.
Our result fundamentally alters the interpretation of analyses such as [21], which places a lower bound on the size of an inflationary patch when it eventually emerges, requiring it to be larger than the background Hubble radius. We demonstrate that such inflationary patches can arise simply through gravitational collapse and redshifting from conditions that, -1 -

JCAP09(2016)010
one or more Hubble times prior to the onset of inflation, were very inhomogeneous. This is a dynamical mechanism that does not require any acausal interaction between regions initially separated by more than one horizon distance.
This picture holds when the range of initial scalar field values lies entirely on the "inflationary plateau" where the potential is flat, satisfying n ≡ (M P ∂ φ ) n ln V 1 (where M P is the Planck mass), with the number of inflationary e-folds in homogeneous universe scaling as N ∼ ∆φ/M P 1 (assumingφ is small initially, but a closely related statement can be made whenφ is large, see e.g. [11]). In the particular class of so-called cosmological attractors (see [10,11] and references therein), which have potentials that rapidly approach a positive constant for large values of the inflaton field, this plateau includes all but a finite range of values -and hence, with high probability the average value of the field falls in a range where the potential is almost exactly a cosmological constant. Nevertheless, we also consider the complimentary case where the spatially averaged value of the scalar field φ lies on the plateau, but the variations in the field fall outside this plateau. A naïve expectation would be that these would also undergo inflation. Instead, we find that when the fluctuations around the average exceed the distance in field space to the end of the plateau, δφ > ∆φ, inflation sometimes does not take place regardless of the value of φ . This is because the field feels the effect of the potential beyond the plateau in at least some parts of the universe, and -depending on the potential -this can have a tendency to pull φ strongly towards the minimum, which is necessarily off the plateau. This condition is of course completely invisible in near-homogeneous cosmologies where δφ is small or vanishing.
Our conclusions are based on a set of fully general relativistic simulations in 3+1 dimensions. The field content is a single scalar φ with a potential and average value that would allow for 60 e-folds or more of inflation given homogeneous initial conditions. However, we begin with inhomogeneous conditions where the energy in spatial gradients dominates over the potential energy by a factor of 10 3 . We simulate an initially expanding universe with toroidal topology where the nonzero wavenumbers k of the initial scalar field content satisfy k/H 0 ≥ 1 (i.e. have wavelength ≤ 2πH −1 0 ), where H 0 is the expansion rate on the initial time slice. Without simulating them, the behavior of longer wavelength modes can be understood as follows. When k/H 1 (a "long" mode), the mode simply renormalizes the local density; that is, so long as they remain outside the horizon, the effect of long modes is captured by homogeneous cosmologies, which are well understood [22]. In fact, the effect of long modes on a volume that has been expanding at some arbitrary initial time will either be to make it closed, in which case it collapses after a finite time, or open, in which case the volume keeps expanding. In a universe that is "flat (or open) on average," it is easy to see that there must exist at least one region where the effect of long modes is to keep it open at all times. This is illustrated in figure 1. Once a long mode re-enters the horizon its evolution is much more complex, but is then described by our simulations.
Instead of considering the periodicity of our simulations as a convenient boundary condition, a complementary scenario covered by this setup is to imagine that the universe is a torus of length of order 2πH −1 0 , with perturbations with wavenumber of order H 0 , which could all be on the order of the Planck scale [6][7][8][9][10][11]. With this interpretation, our simulations describe the whole universe rather than a limited part.

Methodology
In order to simulate an inhomogeneous cosmology, we solve the Einstein field equations coupled to the inflaton φ, which has equation of motion φ = V . We evolve the field -2 - Figure 1. Let us consider a region (labelled "(4)") of size comparable of the would-be inflationary Hubble patch. If the region is open or flat on average, it means that modes longer than the region will not make the whole region collapse. For each subdivision of region (4), there has to be one region, here labelled (3), where the role of longer modes is to make the region open or flat, while its complement (the part of (4) not within (3)) can be closed. The same applies to further subdivisions, here (2) and (1). In this way, given a region (4) that is flat or open on average, there is a sequence of smaller region where the effect of long modes is to keep it open at all times. Therefore, long modes will not prevent geodesics originating in region (1) from eventually inflating. As our simulations show, modes shorter than the horizon at any given time can produce localized black holes, but do not make the region as a whole collapse. equations in the generalized harmonic formulation using the code described in [23,24]. We use units where 8πG = c = 1 throughout.
For the scalar field initial conditions, we choose a superposition of standing waves witḣ φ(t = 0) = 0 and where k ranges over all the allowed wavevectors in a three-dimensional periodic domain (e.g. for N = 1, the six wavevectors given by plus and minus each of the three coordinate directions), and the θ k values are randomly chosen phases. We always take the length of the domain in each direction to be equal to the wavelength of the longest mode in that simulation L = 2π/k min , and consider various ratios of this to the initial Hubble scale. We construct initial data for the metric by solving the constraint equations using the code described in [25].
We assume (at the initial time only) that the spatial metric is conformally flat γ ij = Ψ 4 f ij , and that the trace of the extrinsic curvature K is constant while the traceless part is zero, so that the momentum constraint is trivially satisfied. The value of the conformal factor Ψ then comes from solving the Hamiltonian constraint.

JCAP09(2016)010
We choose K based on the integral of the Hamiltonian constraint over the periodic domain (ignoring the conformal factor) which we find to give solutions with Ψ ∼ 1. The choice of a negative K will give us an initially expanding universe with H 0 = −K/3. (A positive K, on the other hand, would give an initially contracting universe, which would likely result in a crunching-type solution.) In the special case of a homogeneous scalar field configuration (δφ = 0) this choice of initial data will reduce to a time slice of an FRW solution.
We consider several inflationary potentials, the simplest being a cosmological constant which is flat for |φ| 1/λ and approximates a step function for large λ. The last is the "notch" potential describing a family of cosmological attractors [26,27]. These potentials cover a large enough class to allow our conclusions to be quite generic: in the flat regions they are well described by a cosmological constant V (φ) = Λ or 0, while the different ways in which the potentials tend to zero cover a wide range of possibilities. We will make use of several quantities defined in terms of the timelike unit normal to slices of constant coordinate time n a for characterizing the simulations below. The vector n a can be considered the four velocity of a fiducial observer whose worldline is orthogonal to the constant time hypersurfaces. We will use the energy density ρ = n a n b T ab where T ab is the scalar field stress-energy tensor. We will also calculate the local volume expansion θ := ∇ a n a = −K. This will allow us to define a fiducial local Hubble expansion rate H := θ/3 and corresponding number of e-folds of expansion N found by integrating n a ∇ a N = θ/3 (see e.g. [28]). We denote the volume average of a quantity by X := X √ γd 3 x/ √ γd 3 x where γ = det γ ij . See the appendix for details on the numerical method and resolution study results demonstrating convergence.

Results
We perform a number of simulations to study the conditions under which a gradient energy dominated universe, ρ = 10 3 Λ, will transition to exponential expansion. In addition to this ratio of gradient to potential energy, we will quantify the inhomogeneities in the different cases we consider by referring to the ratio of the wavenumbers of the initial scalar field content to the expansion rate on the initial time slice k/H 0 . We first present results using a cosmological constant, and following that present cases with non-trivial potentials where some part of the scalar field range falls outside the flat region of the potential.
We choose V (φ) = Λ and initial data with scalar field perturbations of a single wavenumber magnitude corresponding to N = 1 in eq. We show several cases with inhomogeneities at a single wavenumber magnitude (N = 1), as well as one case with several different wavenumber magnitudes (N = 5) where the smallest wavenumber magnitude is k/H 0 = 4. For comparison, the slope of these quantities expected in a radiation dominated FRW universe is also shown. Once apparent horizons are found, we ignore their interiors for the purpose of calculating these quantities, which accounts for the discontinuous features evident in the k/H 0 = 1 and 2 cases.
(where a = e N ), as would be expected in an FRW universe. We find similar results for initial configurations with more than one wavenumber (N = 5 is also shown in figure 2). This continues until the kinetic/gradient energy is subdominant, leading to a cosmological constant-dominated universe undergoing exponential expansion. As shown in the top of figure 3, though initially there is spatial variation in the energy density and rates of expansion/contraction, these eventually go away.
For smaller values of k/H 0 we are in the strong-field regime and there is prompt gravitational collapse in some regions. In the cases with k/H 0 = 1 and 2 shown in the bottom of figure 3, the maximum energy density rapidly increases and black holes form at t ≈ H −1 0 (we emphasize that with the gauge choices made here, the time coordinate is different from that of the usual FRW metric). This is expected since the hoop conjecture [29] predicts that a trapped surface will form roughly when the mass of an overdensity is comparable to the mass of a black hole of the same size 4 3 πGk −3 ρ ∼ k −1 /2, which is equivalent to k ∼ H. In terms of -5 -

JCAP09(2016)010
the amplitude of the scalar field fluctuations, this is also equivalent to δφ ∼ M P . Because of the symmetry in the initial data and potential with respect to positive and negative values of the scalar field, two identical apparent horizons form in the periodic domain. At the final time shown in figure 3, the irreducible mass of each apparent horizons is ≈ 2k −1 and ≈ 0.6k −1 for k/H 0 = 1 and 2, respectively. We do not track the interior of the apparent horizons and avoid any issues with singularities that may develop there by excising this region from the computational domain. However, ignoring the black hole interiors, even in these cases there is still expansion on average and the volume-averaged density decreases until the domain becomes vacuum energy dominated. The proper distance between the black holes also increases roughly in proportion to the average scale factor, indicating they are becoming more and more isolated in an exponentially expanding universe (thus we end in a situation similar to [30]).
Similar results to those obtained with a cosmological constant potential should hold in general for cases where the entire range of scalar field falls on flat region of the potential. We have explicitly verified that for eq. (2.4) with λ = 1/ √ 6 [26], and φ 0 = 6.2 (corresponding to N ∼ 60 in homogeneous inflation), the results obtained are indistinguishable from those presented above with k/H 0 = 1. We next consider cases where the range of scalar field values does not satisfy this condition, and hence the average value of the scalar field φ 0 will be important. For simplicity we will concentrate on initial scalar field configurations with N = 1 and k/H 0 = 4 and vary φ 0 .
In figure 4 we show several cases with the potential given by eqs. (2.3) with λ = 200, which gives a very good approximation to a step function. Of these, the one with the smallest average scalar field value (φ 0 = 0.0735) would give ∼ 60 e-folds of inflation before rolling down to a negative value in a homogeneous FRW model. However in the inhomogeneous case we consider here, the value of the scalar field everywhere becomes negative before the gradient energy of the scalar field is negligible, and there is no phase of exponential expansion. This also happens for the larger average value of φ 0 /δφ = 0.5. However, for φ 0 /δφ = 1 the scalar field approaches a positive value everywhere as the potential energy dominates and there is an inflationary period.
As also shown in figure 4, similar results are obtained using the potential given by eq. (2.4) with λ = 100. Here again, for the cases where φ 0 /δφ is small, the scalar field moves away from the inflationary plateau to a region of field space with lower potential energy, in this case the narrow region at the minimum of the potential at φ = 0.
This can be seen analytically as follows. In a homogeneous universe, time derivatives of φ are proportional to the "force" V (φ). When δφ = 0, spatially averaged time derivatives are proportional to the averaged force V . In the case of the step, when δφ > φ 0 , V ∼ V / φ and so φ is pushed towards negative values parametrically faster than if δφ < φ 0 . In general, on a slow-roll plateau, M P V /V ∼ √ < 1. If the fluctuations are large enough to reach the minimum of V , the spatially averaged force is . This is the origin of the statement made in the introduction regarding δφ < ∆φ. From an Effective Field Theory (EFT) point of view, this result can be understood by noticing that the effective potential for the homogeneous mode, which evolves with a time scale of order Hubble, is obtained by integrating out the short (and fast) modes which are classically excited. Due to the strong nonlinearities, the effective potential can acquire a width of order the amplitude of the classical modes.
One might expect that a potential with a sharp symmetric minimum could avoid this tendency. Consider the notch potential, (2.4) with λ large. Then V is suppressed by 1/λ, . Each panel has three subplots which, from top to bottom, show as a function of time: the energy density minus the vacuum energy contribution, normalized by the vacuum energy density; a measure of the effective scale factor (normalized to be unity at t = 0); the effective expansion rate, normalized by the Hubble constant of a vacuum energy dominated universe. We show the maximum, minimum (omitted from the top panel), and volume-averaged value of the given quantities. In the second subplot, in addition to exp(N ), we also show the cube-root of the spatial volume. We ignore apparent horizon interiors when calculating these quantities, which accounts for the discontinuous features in the bottom two panels.  and so one might expect a sharp notch to have only a small effect, which is indeed the case in the absence of gravity. However, with gravity included the combination of Hubble friction and gravitational nonlinearities rapidly pulls φ into the minimum. The rate at which φ decreases is almost independent of λ for large λ. From an EFT point of view, in an expansion in k/(aH) 1, we have that the leading effect vanishes, but higher orders do not.

Conclusions
For a large class of examples, we find that exponential expansion occurs even when the initial gradients are much larger than the potential energy. The possible exceptions are cases when the range of the initial inhomogeneous scalar field values exceeds the inflationary plateau, in which case the outcome depends on the details of the potential. This is illustrated by the potential with a nearly flat plateau for φ > 0, and a relatively sharp "step" down to zero at φ < 0 given by (2.3) with λ large: when δφ φ 0 the large V at the step near φ = 0 has a strong effect on the time evolution of φ , rapidly propelling it to negative values and preventing inflation from beginning. We provide an analytical understanding of this behavior. We note, however, that this scenario does not apply to generic examples of large-field inflation where φ 0 is large in units of M P , but the fluctuations about this value are order one or smaller.
When the fluctuations in the inflaton field are contained within a flat region of the potential, the potential will essentially act like a cosmological constant, in which case in any region of the universe there are two possible outcomes: recollapse after a finite time, or expansion until the vacuum energy dominates and inflation begins. With our initial conditions -large fluctuations of wavelength ≤ 2π/H 0 in expanding universes that are flat (or open) on average -we find that both types of regions occur, with the former leading to black holes, and the latter eventually dominating the volume. Thus the amplitude of the initial-state fluctuations in φ is irrelevant except in determining the time scales. In particular, we find that in order for inflation to start somewhere, there is no need to assume a Hubble-sized homogeneous initial patch.
An important issue is how natural it might be to have an inflaton potential and initial field range that satisfy the criteria laid out above. That question is beyond the scope of this work and we do not address it here, though we refer the reader to [11,[31][32][33] and references therein which discuss how cosmological attractor configurations with the desired properties, such as approximate shift symmetric potentials, can be constructed in supergravity, Higgs inflation, and other contexts. Additionally, exploring a broader class of initial conditions, including multiple fields, or the effects of homogeneous kinetic energy, will be interesting follow-up work.

A Details of numerical methods
In order to numerically simulate an inhomogeneous cosmology, we solve the Einstein field equations coupled to the inflaton φ in a periodic domain. The field equations are evolved in the generalized harmonic formulation, using the code described in [23,24]. We use the same variation of the damped harmonic gauge [34,35] as in [36]. The inflaton equation of motion φ = V (where V is the potential) is evolved using the same fourth-order finite difference stencils and Runge-Kutta time stepping as the metric. As the simulations evolve and the metric components grow due to expansion, we dynamically adjust the timestep size in proportion to the decreasing global minimum of γ 1/6 /α (where γ is the determinant of the spatial metric and α is the lapse) in order to avoid violating the Courant-Friedrichs-Lewy condition. The simulations are performed with between 256 and 512 points across each linear dimension to establish numerical convergence and estimate error. The convergence of the constraints is shown for the most strong-field case considered above in figure 5. In figure 6, we show the truncation error in the volume-averaged energy density and expansion rate for the k/H 0 = 4 case, which are sub-percent level even at the lowest resolution.