Robustness of slow contraction to cosmic initial conditions
Anna Ijjas, William G. Cook, Frans Pretorius, Paul J. Steinhardt, Elliot Y. Davies
PPrepared for submission to JCAP
Robustness of slow contraction tocosmic initial conditions
Anna Ijjas, a,b, William G. Cook, c Frans Pretorius, d Paul J.Steinhardt, d and Elliot Y. Davies d a Max Planck Institute for Gravitational Physics (Albert Einstein Institute), D-30167 Han-nover, Germany b Gottfried Wilhelm Leibniz Universität, D-30167 Hannover, Germany c Theoretisch-Physikalisches Institut, Friedrich-Schiller-Universität, D-07743 Jena, Germany d Department of Physics, Princeton University, Princeton, NJ 08544, USAE-mail: [email protected]
Abstract.
We present numerical relativity simulations of cosmological scenarios in which theuniverse is smoothed and flattened by undergoing a phase of slow contraction and test theirsensitivity to a wide range of initial conditions. Our numerical scheme enables the variation ofall freely specifiable physical quantities that characterize the initial spatial hypersurface, suchas the initial shear and spatial curvature contributions as well as the initial field and velocitydistributions of the scalar that drives the cosmological evolution. In particular, we includeinitial conditions that are far outside the perturbative regime of the well-known attractorscaling solution. We complement our numerical results by analytically performing a completedynamical systems analysis and show that the two approaches yield consistent results. Corresponding Author a r X i v : . [ g r- q c ] J u l ontents ε > The cosmic initial conditions problem presents an unusual theoretical puzzle. Typically, inunderstanding a dynamical system, the challenge lies in finding the right (differential) equa-tions that describe how the system evolves with time. The specific point at which one startssolving the evolution equations is of no particular importance, especially because the intentis to apply the same equations to understanding the behavior over a wide range of systemsand initial conditions. In cosmology, on the other hand, the relevant dynamical equationsare well-known. Measurements of the cosmic microwave background and other astrophys-ical experiments confirm that Einstein’s classical theory of general relativity describes thetime evolution of our large-scale universe since primordial nucleosynthesis to an astonishingaccuracy, if we specify the initial geometry and radiation-matter content. The very same ob-servations, though, also teach us that the cosmic initial conditions cannot be explained by the– 1 –ame physics that underlies the evolution equations because, in that context, they are expo-nentially rare or finely tuned. Some mechanism operating before primordial nucleosynthesisis needed to explain how these initial conditions arise.The goal of this paper is to show that slow contraction is a ‘robust’ smoothing mechanismthat can naturally generate the cosmic initial conditions. The slow contraction phase isinduced by a canonical scalar field evolving down a steep negative potential coupled to Einsteingeneral relativity, as postulated in many bouncing and cyclic cosmological models [17, 18, 20,26]. We adapt the tools of numerical general relativity to solve the full set of coupled Einstein-scalar field equations beginning from a wide range of highly inhomogeneous and anisotropicstates that would not be compatible with observations of the cosmic microwave backgroundand large-scale structure if there were not a robust smoothing and flattening mechanism. Themeasure of ‘robustness’ is how sensitive the outcome is to the initial state and whether theoutcome converges closely and rapidly to the initial conditions needed to explain cosmologicalobservations.More precisely, the large-scale universe evolves as predicted by the Einstein equations if the cosmic initial conditions at primordial nucleosynthesis are very precisely set as follows:i. the background spacetime geometry is smooth and spatially flat described by a sin-gle dynamical quantity, the Friedmann-Robertson-Walker (FRW) scale factor a ( τ ) , orequivalently, the associated Hubble radius Θ = | H − | . (Here and throughout, the Hub-ble parameter H is defined as the logarithmic time derivative of the scale factor H = a (cid:48) /a and prime denotes differentiation with respect to the proper FRW time coordinate τ .);ii. to first sub-leading order, there is only a single kind of deviation from this geometry,namely a nearly scale-invariant and gaussian spectrum of adiabatic curvature fluctua-tions;iii. there are no measurable first-order tensor, vector, or isocurvature fluctuations in theinitial space-time geometry;iv. the initial background geometry as well as the initial spectrum of curvature fluctuationsare correlated over (or, equivalently, e ) Hubble-sized patches.Notably, under these initial conditions, the Einstein equations dramatically simplify and thesubsequent large-scale evolution of 14 billion years is given by the Friedmann equation, d H − d τ = ε eff , (1.1)where ε eff ≡ (3 / p/(cid:37) ) is the equation of state of the dominant stress-energy component,relating its pressure p to its energy density (cid:37) . (Throughout, we work in reduced Planckunits.) Furthermore, as the background slowly expands, the initial curvature fluctuations inthe space-time geometry grow and become the seeds from which galaxies, stars and planetsform.However, in a classical general relativistic space-time, this combination of initial con-ditions is highly non-generic if the stress-energy is sourced by ordinary (baryonic or dark)matter and radiation. Accordingly, the cosmic initial conditions problem consists in identify-ing a mechanism that generically (i.e., for most initial data) leads to the conditions as giventhrough (i.)-(iv.). – 2 –he traditional approach to resolving the cosmic initial conditions problem has beento propose a ‘ classical ’ smoothing mechanism [8]. For this purpose, it is necessary that themechanism possesses a dynamical attractor solution in which the relative contribution of small inhomogeneities and anisotropies to the total energy density shrinks with time. In all knownexamples, classical smoothing is achieved by introducing a novel stress-energy component,typically sourced by scalar field φ , that evolves to dominate all other contributions to thegeneralized Friedmann constraint, H = 13 (cid:16) (cid:37) m a + (cid:37) r a + (cid:37) φ a ε (cid:17) − ka + σ a . (1.2)Here, (cid:37) m , (cid:37) r , (cid:37) φ represent the energy density of matter, radiation and the scalar field φ ,respectively, at some initial time t , and the scale factor is normalized such that a ( t ) = 1 .The last two terms correspond to spatial curvature and anisotropy.From Eq. (1.2), it is immediately clear that there can only be two kinds of classicalsmoothing mechanisms: • in an expanding universe ( H > ): ≤ ε < ( inflation ); • in a contracting universe ( H < ): ε > ( slow contraction ).The respective equation of state ε is achieved by choosing a suitable potential energy density V ( φ ) . For example, assuming a scalar field with canonical kinetic energy density and anegative exponential potential V ( φ ) = V exp( − φ/M ) , (1.3)where V < , the contracting FRW space-time is stable to small perturbations and hence adynamical attractor solution of the Einstein-scalar field equations, if the characteristic energyscale M < √ . In this case, M − = √ ε .In addition, both inflationary expansion and slow contraction transform quantum fluc-tuations generated on scales smaller than the Hubble radius into ‘squeezed modes’ on scaleslarger than the Hubble radius, as required by condition (ii). During slow contraction, for ex-ample, the characteristic wavelength of a fluctuation, λ , decreases in proportion to the scalefactor, λ → λ × a ( t ) ; but the Hubble radius decreases more rapidly, as Θ ∼ a ε where ε > .Consequently, a mode that originates on scales much smaller than a Hubble patch evolves toa wavelength extending over scales exponentially larger than the Hubble radius.Yet, classical smoothing is necessary but not sufficient by far to explain our large-scaleuniverse. Satisfying this criterion alone would mean that only an infinitesimal set of initialconditions, namely small classical perturbations around FRW, lead to a universe as we observeit. Another, obvious requirement is that the classical evolution remains stable to quantumfluctuations. Notably, inflation is a classical but famously not a quantum smoother [16, 27, 30];but slow contraction is both a classical and a quantum smoother [8].More importantly, a classical smoothing phase, even if it is stable to quantum fluctu-ations, does not solve the cosmic initial conditions problem if it is only stable under smallperturbations around a smooth and flat FRW background. After all, the set of space-timegeometries that represent small deviations from a FRW space-time only represent a measure-zero set of all permitted and (physically) plausible initial conditions. To generically drivethe universe towards a flat, homogeneous and anisotropic space-time, a classical smoothing– 3 –echanism must be a ‘ robust ’ smoother, i.e. , insensitive to a wide range of arbitrary initialconditions including those well outside the perturbative regime of the attractor FRW solution.In this paper, we examine quantitatively the robustness of slow contraction using anumerical scheme that enables the variation of all freely specifiable physical quantities thatcharacterize the initial spatial hypersurface, such as the initial shear and spatial curvaturecontributions as well as the initial field and velocity distributions of the scalar that drivesthe cosmological evolution. In fact, the only restrictions we have on the initial data resultfrom imposing periodic spatial boundary conditions and choosing an initial three-metric thatis conformally flat.In particular, we ‘empirically’ confirm the well-known ultra-locality conjecture [6] anddemonstrate that, generically, all gradients rapidly become negligible as the evolution pro-ceeds. Finally, we show that the homogenous end states we identified numerically are theonly stable attractor solutions of the underlying relativistic system of evolution and con-straint equations. To study the robustness of slow contraction in smoothing and flattening the universe, we solvethe full Einstein-scalar field equations for a wide range of highly non-perturbative initial con-ditions using the techniques of numerical general relativity. The study entails numericallyevolving a system of coupled non-linear, second-order partial differential equations. Per-forming such a computation necessitates finding a ‘good’ formulation of the field equationssatisfying two criteria: • the formulation is ‘well-suited’ to the physical situation; and • the formulation is well-posed.The term ‘formulation’ is given multiple definitions in the literature. Here, we follow theconvention established in mathematical relativity ( e.g. , see [15]): A ‘formulation’ (sometimescalled an initial value formulation) of a given theory is the representation of the underlyingsystem of differential equations obtained by choosing a particular coordinate system. Fur-thermore, we define a formulation to be ‘well-posed’ if the underlying system of differentialequations can be put into a form such that, for given initial conditions, there exists a uniquesolution that depends continuously on the initial conditions. Note that a formulation is notequivalent to a gauge choice. Different gauge conditions can be implemented in the sameformulation, but not all gauge choices lead to a well-posed formulation. Typical cosmological studies in the literature are perturbative and only consider thesecond criterion. The diffeomorphism invariance of the field equations is exploited to simplifythe theoretical analysis and straightforwardly relate the predictions of different cosmologicalmodels to observations.For example, in studying the evolution of perturbations away from FRW in the veryearly universe, a coordinate basis using the well-known (3 + 1) or Arnowitt-Deser-Misner(ADM) decomposition [2] proves to be well-suited. The ADM formulation naturally restson the homogeneity and isotropy of the FRW background solution, enabling the separation We note that some fields use the phrase ‘well-posed problem’ as we use ‘well-posed formulation’ here. Thelatter is a common usage in the modern relativity literature and implicitly includes that appropriate boundaryconditions and initial data are specified. – 4 –f linear perturbations around this background into decoupled scalar, vector, and tensordegrees of freedom that each evolve independently mode by mode [3, 14, 22]. Finally, asuitable gauge choice is made for the physical situation at hand. In a gauge such as unitary,the scalar part of the linearized spatial metric component can be identified with the co-moving curvature perturbation, an invariant of the linear theory that directly determines thetemperature anisotropies of the cosmic microwave background (CMB) [5, 23].However, for the non-perturbative numerical analyses presented in this paper, the com-mon procedure adapted in the cosmology literature is insufficient. Without a well-posedformulation of the full (non-perturbative) Einstein-scalar field equations, the existence of aunique solution is not guaranteed. For a given set of initial conditions, the system of equationsmight admit no solution at all or it might admit multiple solutions. That is, in an ill-posedformulation, no predictions can be derived. A common manifestation is the ‘blow up’ of thenumerical code within finite time even if there is no fundamental instability in the underlyingtheory.Notably, many formulations of the Einstein-scalar field equations, including the ADMdecomposition popular with cosmologists, are ill-posed . In particular, for the case where thelapse and shift are given by algebraic relations, it is straightforward to show that the resultingpartial differential equations for the evolution subsystem of the Einstein-scalar field systemis only weakly hyperbolic rather than well-posed (strongly hyperbolic). Strictly speaking,cosmologists focusing on perturbative cosmological analyses are only able to get by with anADM decomposition because there exist other formulations of Einstein scalar-field equationsthat are well-posed and that have been shown to uniquely admit the FRW solution assumedin ADM.In general, we must not only choose a well-posed formulation of the field equations but, insome cases, one that also admits an effective constraint damping scheme. Constraint dampingis a common tool used in numerical relativity. In principle, by choosing initial conditions thatsatisfy the Hamiltonian and momentum constraints, i.e. , the Einstein equations projectedorthogonally onto a spacelike hypersurface at some initial time t , the field equations ofgeneral relativity propagate and preserve the constraints going forward in time. In practice,though, the constraints can only be satisfied initially up to some numerical error, and, insome cases, the error grows exponentially. Then, even a well-posed formulation can resultin a numerical blow-up if these constraint violations are not addressed. For this reason, insome numerical setups it is often prudent to add terms to the evolution equations that dampthe growth of constraint violating numerical errors, while leaving the underlying solutionunaffected in the continuum limit [24, 25]. As it turns out, in the numerical scheme describedin this paper, constraint damping is not required for stability of the numerical evolution.In the cases considered in this paper, we need a ‘suitable’ formulation of the equationsthat address issues specifically related to slowly contracting spacetimes: a stiffness problemand accurately tracking contraction over many e-folds before reaching a big crunch. Thestiffness problem arises because the equation of state parameter ε is greater than three (withtypical models having ε (cid:29) ) which means that Θ ∝ a ε is very rapidly decreasing comparedto a ( t ) . For example, in models discussed in the literature [18], during a period when theaveraged scale factor decreases a factor of two or three, the Hubble radius shrinks by a factorof e or more. To handle this stiffness problem, we choose a Hubble-normalized formulationin which the Hubble radius Θ does not enter explicitly.To handle the crunch problem, we choose a time coordinate such that it follows themean curvature growth, i.e. a time slicing where e t ≡ . This way, reaching the curvature– 5 –ingularity (or vanishing of the Hubble radius Θ ) takes an infinite (coordinate) time. We findthat our Hubble-normalized formulation using this time slicing is sufficiently suitable (ourthird criterion above) for analyzing slow-contraction.We find that an effective way of incorporating these methods of handling the stiffnessand crunch problems within a well-posed initial value formulation is to adapt the orthonormaltetrad formulation of the Einstein-scalar field equations in our numerical scheme. Originally,the formulation was used by Schücking to find all exact vacuum solutions describing spatiallyhomogeneous spacetimes [11, 21]. Later, it was successfully implemented to numericallystudying contracting vacuum spacetimes [4, 7, 12] as well as spacetimes with a canonicalscalar field [13]. In Sec. 2.1, we first introduce the basic tetrad variables. Then, we derive theEinstein scalar field equations in tetrad form in Sec. 2.2. Finally, in Sec. 2.3, we convert thetetrad equations to partial differential equations using local coordinates, making the systemreadily usable as a numerical evolution scheme. Tetrad formulations of the Einstein-scalar field equations assign each spacetime point a familyof unit basis 4-vectors (or vierbeins ) { e , e , e , e } (as opposed to coordinates { x , x , x , x } )that describe local Lorentz frames with the spacetime metric being given by the dot productof the basis vectors. The starting point is a timelike vector field e that defines a future-directed timelike reference congruence, to which it is tangent. It is supplemented by a triadof spacelike unit 4-vectors { e , e , e } that lie in the rest 3-spaces of e .In an orthonormal tetrad formulation that we shall employ, the spacetime metric iseverywhere given by e α · e β = η αβ , (2.1)where η αβ = diag( − , , , the Minkowski metric and " · " is the spacetime inner product.Throughout, spacetime indices (0 − are Greek and spatial indices (1-3) are Latin. Thebeginning of the alphabet ( α, β, γ or a, b, c ) is used for tetrad indices and the middle of thealphabet ( µ, ν, ρ or i, j, k ) is used for coordinate indices. Tetrad frame indices are raised andlowered with η αβ .The geometric variables of the formulation are the sixteen tetrad vector components andthe twenty-four Ricci rotation coefficients γ αβγ ≡ e α · ∇ γ e β = − γ βαγ , (2.2)which define the deformation of the tetrad when moving from point to point. Here, ∇ γ is thespacetime covariant derivative projected onto a tetrad e γλ ∇ λ . Note that the γ αβγ are the‘tetrad components’ of the Christoffel symbols Γ µνρ .Similar to coordinate-based formulations of the field equations, the tetrad formulationgreatly simplifies when making a space-time split. Unlike in the 3+1 (coordinate based) ADMformulation, where the split is defined by the constant-time spacelike hypersurface, here thesplit is relative to the timelike congruence defined by e . Note, though, that the timelikecongruence does not uniquely define the auxiliary spatial congruence. Rather, the spatialtriad vectors are fixed by imposing gauge conditions.To perform the spacetime split, we first write the fifteen connection coefficients that haveat least one timelike index in terms of 3-dimensional quantities b a , Ω a , and K ab , reflecting theantisymmetry of γ αβγ in its first two indices: γ a = − γ a = b a , (2.3)– 6 – ab = − γ ba = (cid:15) abc Ω c , (2.4) γ ab = − γ a b = − K ba ; (2.5)where (cid:15) abc is the Levi-Civita symbol. As shown by Ehlers et al. [9, 19], the fifteen quantitiesdescribe kinematic fields associated with the timelike congruence tangent to e : the 3-vector b a is the local proper acceleration; the 3-vector Ω a is the local angular velocity of the space-like triads { e , e , e } relative to Fermi-propagated axes; and K ba is the local rate-of-strain(or shear) tensor. (A spatial triad is ‘Fermi-propagated’ if Ω a ≡ , i.e. , it is a local, inertiallynon-rotating reference frame.)Similarly, we express the remaining nine purely spatial connection coefficients γ abc thatdescribe the induced curvature of the auxiliary 3-congruence using a 3-tensor, N ab = (cid:15) bcd γ cda . (2.6)The spatial connection coefficients N ab and the components of the shear tensor K ab are theeighteen dynamical variables. The acceleration and angular velocity vectors, b a , Ω a , are thesix (tetrad) gauge source functions.For completeness, we note that some authors use as basic variables the commutation (orstructure) coefficients C αβγ rather than the Ricci rotation coefficients γ αβγ [10, 28, 29]. Here,the C αβγ are given through [ e α , e β ] ≡ ∇ α e β − ∇ β e α ≡ C αβγ e γ . (2.7)It is straightforward to relate the two conventions: With the definitions in Eqs. (2.2) and(2.7), γ αβγ = 12 (cid:16) C γβα + C αγβ − C βαγ (cid:17) , (2.8)or, alternatively, C αβγ = γ γβα − γ γαβ , (2.9) i.e. , expressed in terms of the twenty-four 3D quantities, the full set of commutation coeffi-cients takes the following form, C a = −C a = b a , (2.10) C ab = −C ba = 2 (cid:15) abc ω c , (2.11) C ab = −C a b = − K ( ab ) + (cid:15) abc ( ω c − Ω c ) , (2.12) C abc = −C bac = (cid:15) cbd N ad . (2.13)Here, the 3-vector ω a is the antisymmetric part of K ab , ω a ≡ (cid:15) abc K bc ; (2.14)it measures the vorticity (or twist) vector of the e -congruence.Finally, the geometric variables have to be supplemented by the dynamical quantitiesthat describe the matter source. In our case, this is the canonical scalar field φ which isspecified in the Einstein-scalar field equations through its potential energy density V ( φ ) .– 7 – .2 Tetrad equations Next we present the tetrad evolution and constraint equations for the (eighteen) dynamicalvariables K ab , N ab . The remaining (six) gauge variables are fixed by our tetrad frame gaugechoice.A natural gauge choice is a frame withi. Fermi-propagated axes ( Ω a ≡ ); andii. hypersurface orthogonal timelike congruence ( ω a ≡ or, equivalently, K ab ≡ K ( ab ) ).Here and throughout, parentheses denote symmetrization, i.e. , K ( ab ) ≡ ( K ab + K ba ) . Inthis (frame) gauge, the time-like vierbein e is the future-directed unit normal to the space-like hypersurfaces Σ t of constant time, and the spatial tetrad vectors are tangent to { Σ t } .Furthermore, K ab is the extrinsic curvature of Σ t and the N ab are the nine (intrinsic) spatialcurvature variables. All Ricci rotation coefficients, K ab , N ab act as scalars on Σ t .Employing these gauge conditions, the tetrad evolution and constraint equations takethe following form: D K ab = (cid:15) acd D c N db + D a b b + b a b b − (cid:15) bcd N ac b d + N cc N ab − K ac K cb − N ca N cb (2.15) + (cid:15) adf (cid:15) bce ( K dc K fe − N dc N fe ) + s ab − δ ab (cid:0) (cid:37) + 3 p (cid:1) ,D N ab = − (cid:15) acd D c K db + (cid:15) bcd K ac b d − N cc K ab + 2 N c [ a K b ] c + (cid:15) adf (cid:15) bce N dc K fe (2.16) − (cid:15) abc j c , D b A b = N ab N ab + (cid:16) K ab K ab − N ab N ab − ( K aa ) − ( N aa ) (cid:17) + (cid:37), (2.17) D b K ab − D a K cc = (cid:15) abc K bd N dc + 2 K ac A c − j a , (2.18) D b N ab − D a N cc = − (cid:15) abc N bd N dc , (2.19)where A b ≡ (cid:15) bcd N cd . (2.20)is the antisymmetric part of N ab ; D is the Lie derivative along the timelike vierbein e ;and D a denotes the directional derivative along the spatial vierbein e a . The matter variablesassociated with the stress-energy T αβ , such as the energy density (cid:37) , pressure p , 3-momentumflux j a , and (spatial) stress tensor s ab , are defined as follows: (cid:37) ≡ e α e β T αβ , (2.21) j a ≡ − e α e aβ T αβ , (2.22) s ab ≡ e aα e bβ T αβ , (2.23) p ≡ s cc . (2.24)Notably, there is no evolution equation for the acceleration 3-vector b a , which reflects the factthat it is a (frame) gauge source function.The tetrad evolution and constraint equations (2.15-2.19) were first obtained in Ref. [7].Their derivation is straightforward when using the tetrad form of the Riemann tensor, R αβγδ = D γ γ αβδ − D δ γ αβγ + γ α(cid:15)γ γ (cid:15)βδ − γ α(cid:15)δ γ (cid:15)βγ + γ αβ(cid:15) ( γ (cid:15)γδ − γ (cid:15)δγ ) . (2.25)More exactly, substituting Eq. (2.25) into the (trace-reversed) Einstein-scalar field equations, R αβ = T αβ − η αβ η γδ T γδ , (2.26)– 8 –here R αβ ≡ R γαγβ is the Ricci tensor, yields the evolution equation (2.15) for K ab as wellas the Hamiltonian and momentum constraints, Eqs. (2.17) and (2.18), respectively. Theevolution and constraint equations for the spatial curvature variables N ab , Eqs. (2.16, 2.19),follow from the Riemann identities, R αβγδ = R γδαβ , (2.27) R αβγδ + R αδβγ + R αγδβ = 0 . (2.28)For a single scalar field with canonical kinetic energy density and (non-zero) potential V ( φ ) , the hydrodynamical (macroscopic) matter variables (cid:37), p, j a and s ab are given by (cid:37) = D φD φ + D a φD a φ + V ( φ ) , (2.29) s ab = D a φD b φ + (cid:0) D φD φ − D c φD c φ − V ( φ ) (cid:1) δ ab , (2.30) p = D φD φ − D a φD a φ − V ( φ ) , (2.31) j a = − D φD a φ. (2.32)Note that, in general, j a is non-zero and s ab is non-diagonal, which reflects the fact thatchoosing a hypersurface-orthogonal tetrad frame gauge generally does not coincide with theco-moving frame of the (scalar field) matter source. In fact, a hypersurface-orthogonal frame-gauge is co-moving only in the homogeneous (ultra-local) limit.The system (2.15-2.19) is completed by adding the scalar-field evolution and constraintequations: D φ = W, (2.33) D a φ = S a , (2.34) D W = − δ ab K ab W + D a S a + ( b a − A a ) S a − V, φ , (2.35) D S a = D a W + b a W − K ( ab ) S b . (2.36)Here, Eqs. (2.33) and (2.34) are the defining relations for the auxiliary variables W and S a , which denote the velocity and gradient of the scalar field φ , respectively. Eq. (2.35) isobtained by expanding the Laplacian of the Klein-Gordon equation ( (cid:3) φ = V, φ ) using theRicci rotation coefficients; and Eq. (2.36) is obtained by evaluating the commutation relation [ e , e a ] φ = −C a W + C ab S b using Eqs. (2.10-2.13). In order to evolve the tetrad equations (2.15-2.19, 2.33-2.36) numerically, we must write themas a system of partial differential equations. That means, we must give a representation ofthe tetrad vector components { e α } using a particular set of local coordinates { x µ } and thenconvert the directional derivatives D α in the tetrad equations to partial derivatives alongthese coordinates.To this end, we introduce the transformation matrix { λ µα } between coordinate and tetradbasis vectors defined thru e α = λ µα e µ . (2.37)The coordinate metric components are then given by g µν = η αβ λ µα λ νβ ; (2.38)– 9 –nd directional derivatives along tetrads can now be written as partial derivatives alongcoordinate directions, D = N − (cid:0) ∂ t − N i ∂ i (cid:1) , D a = E ai ∂ i , (2.39)where N is the tetrad lapse function and the N i are the three coordinate components of thetetrad shift vector. Both the tetrad lapse function and the tetrad shift vector describe theevolution of the coordinates relative to the tetrad congruence (as opposed to the ADM lapseand shift that describe the evolution of the proper time and co-moving spatial coordinatesrelative to the coordinates of a particular foliation). The nine coordinate components E ai describe projections of the spatial tetrads tangent to the constant-time hypersurface Σ t . Notethat the spatial triad vectors have zero time component, since we have chosen our tetrad framegauge to be hypersurface-orthogonal. (For arbitrary tetrad frame gauge choices, this is notthe case.)The E ai are dynamical variables determined by the evolution and constraint equations N − ∂ t E ai = − K ac E ci , (2.40) (cid:15) cab E ai ∂ i E bj = N dc E dj − N dd E cj , (2.41)which are derived from applying the commutators of the basis vectors to the spatial coordi-nates { x i } .The lapse function N and the shift vector N i are gauge variables that we can freely spec-ify. A natural coordinate gauge choice when studying the evolution of cosmological spacetimesis to have co-moving coordinates, such that the x i are constant along the congruence and N i ≡ . (2.42)For the lapse function, we will impose the condition that hypersurfaces Σ t of constanttime be constant mean curvature (CMC) hypersurfaces, i.e. , the trace of the extrinsic curva-ture K ab is spatially uniform for each Σ t , Θ − ≡ − K aa (cid:12)(cid:12) Σ t = const > . (2.43)Imposing CMC slicing, the trace of Eq. (2.15) reduces to a linear, elliptic equation forthe lapse N , (cid:16) − (cid:0) D a − A a (cid:1) D a + Σ ab Σ ab + 3Θ − + W − V ( φ ) (cid:17) N = 0 , (2.44)where Σ ab ≡ K ab − K cc δ ab (2.45)is the trace-free part of the extrinsic curvature. Few works that rigorously prove well-posednesstreat elliptic gauge conditions, and we are not aware of any for the particular tetrad formula-tion we use; however, see Ref. [1] for a proof in a closely related coordinate based formulation.We note, though, that if a system of partial differential equations is not well posed it would bechallenging, if not impossible, to obtain stable, convergent numerical solutions with any dis-cretization scheme. That our code is stable and convergent can thus be viewed as ‘empiricalevidence’ that the underlying equations are well posed.Note that the (elliptic) equation for the (tetrad) lapse also determines the tetrad gaugesource function b a , which denotes the acceleration of the tetrad congruence worldlines. Thiscan be seen, e.g. , by computing the commutator [ e , e a ] as applied to the time coordinate x , b a = N − E ai ∂ i N. (2.46)– 10 –or the time coordinate t , we choose a particular (re-)scaling, e t = 3Θ , (2.47)that is consistent with the CMC slicing condition. If Eq. (2.43) is satisfied, the inverse trace ofthe extrinsic curvature (here and throughout denoted by Θ ) is the well-known Hubble radius(as measured by the proper time coordinate τ ), i.e. , e − t = d ln a ( τ ) dτ . (2.48)During contraction, the Hubble radius decreases. Accordingly, the time coordinate t ≤ ,running from small negative to large negative values. Yet, due to CMC slicing, all curvaturevariables remain finite and non-zero for any evolution of finite duration. This is a necessarycondition to stably evolve contracting phases that last several hundreds of e -folds, which isrequired to study the robustness of slow contraction.In addition, for a sufficiently long evolution, we must ensure that no two dynamicalvariables grow (or shrink) at significantly different rates to avoid a stiffness problem . Since inthe case of slow contraction, some spatial metric variables, such as the spatially averaged scalefactor, decrease at a significantly lower rate than some curvature variables, such as the Hubbleradius, we eliminate the latter by introducing dimensionless Hubble-normalized variables, N → N ≡ N/ Θ , (2.49) (cid:110) E ai , Σ ab , A a , n ab , W, S a , V (cid:111) → (cid:110) ¯ E ai , ¯Σ ab , ¯ A a , ¯ n ab , ¯ W , ¯ S a , ¯ V (cid:111) , (2.50)where N is the Hubble-normalized lapse function; bar denotes multiplication by the Hubbleradius Θ ; and n ab ≡ N ( ab ) (2.51)is the symmetric part of N ab . Substituting into Eq. (2.44), yields an elliptic equation for theHubble-normalized lapse N− ¯ E ai ∂ i (cid:0) ¯ E aj ∂ j N (cid:1) + 2 ¯ A a ¯ E ai ∂ i N + N (cid:16) ab ¯Σ ab + ¯ W − ¯ V (cid:17) = 3 . (2.52)Imposing CMC slicing as defined in Eqs. (2.43, 2.47) and using Hubble-normalized vari-ables in Eqs. (2.15-2.16, 2.33-2.36, and 2.40), the gravitational quantities ¯ E ai , ¯Σ ab , ¯ n ab , and ¯ A a as well as the scalar field matter variables φ, ¯ W , ¯ S a satisfy the hyperbolic evolution equations ∂ t ¯ E ai = − (cid:16) N − (cid:17) ¯ E ai − N ¯Σ ab ¯ E bi , (2.53) ∂ t ¯Σ ab = − (cid:16) N − (cid:17) ¯Σ ab − N (cid:16) n (cid:104) ac ¯ n b (cid:105) c − ¯ n cc ¯ n (cid:104) ab (cid:105) − ¯ S (cid:104) a ¯ S b (cid:105) (cid:17) + ¯ E (cid:104) ai ∂ i (cid:16) ¯ E b (cid:105) i ∂ i N (cid:17) (2.54) − N (cid:16) ¯ E (cid:104) ai ∂ i ¯ A b (cid:105) − (cid:15) cd ( a (cid:16) ¯ E ci ∂ i ¯ n b ) d − A c ¯ n b ) d (cid:17)(cid:17) + (cid:15) cd ( a ¯ n b ) d ¯ E ci ∂ i N + ¯ A (cid:104) a ¯ E b (cid:105) i ∂ i N ,∂ t ¯ n ab = − (cid:16) N − (cid:17) ¯ n ab + N (cid:16) n ( ac ¯Σ b ) c − (cid:15) cd ( a ¯ E ci ∂ i ¯Σ b ) d (cid:17) − (cid:15) cd ( a ¯Σ b ) d ¯ E ci ∂ i N , (2.55) ∂ t ¯ A a = − (cid:16) N − (cid:17) ¯ A a − N (cid:16) ¯Σ ab ¯ A b − ¯ E bi ∂ i ¯Σ ab (cid:17) − ¯ E ai ∂ i N + ¯Σ ab ¯ E bi ∂ i N , (2.56) ∂ t φ = N ¯ W , (2.57) ∂ t ¯ W = − (cid:16) N − (cid:17) ¯ W − N (cid:16) ¯ V ,φ + 2 ¯ A a ¯ S a − ¯ E ai ∂ i ¯ S a (cid:17) + ¯ S a ¯ E ai ∂ i N , (2.58)– 11 – t ¯ S a = − (cid:16) N − (cid:17) ¯ S a − N (cid:16) ¯Σ ab ¯ S b − ¯ E ai ∂ i ¯ W (cid:17) + ¯ W ¯ E ai ∂ i N . (2.59)(Here angle brackets denote traceless symmetrization defined as X (cid:104) ab (cid:105) ≡ X ( ab ) − X cc δ ab .)The system of Eqs. (2.52-2.59) will serve as our numerical scheme. Notably, this systemis well-posed, as was shown, e.g. , in [4], and, by construction, it admits stable numericalevolution that involves several hundreds of e -folds of contraction.In addition, the same variables satisfy the constraint equations C G ≡ E ai ∂ a ¯ A a − A a ¯ A a − ¯ n ab ¯ n ab + (¯ n cc ) − ¯Σ ab ¯Σ ab (2.60) − ¯ W − ¯ S a ¯ S a − ¯ V = 0 , ( C C ) a ≡ ¯ E bi ∂ i ¯Σ ab − ab ¯ A b − (cid:15) abc ¯ n bd ¯Σ cd − ¯ W ¯ S a = 0 , (2.61) ( C J ) a ≡ ¯ E bi ∂ i ¯ n ba + (cid:15) bca ¯ E bi ∂ i ¯ A c − A b ¯ n ba = 0 , (2.62) ( C S ) a ≡ ¯ S a − ¯ E ai ∂ i φ = 0 , (2.63) ( C E ) ia ≡ (cid:15) bca (cid:16) ¯ E bj ∂ j ¯ E ci − ¯ A b ¯ E ci (cid:17) − ¯ n ad ¯ E di = 0 . (2.64)The constraints are the Hubble-normalized version of Eqs. (2.17-2.19, 2.34, 2.41), again im-posing CMC slicing conditions as in Eqs. (2.43, 2.47). (The subscripts G, C and J stand forGauss, Codazzi, and Jacobi, respectively, referring to the commonly used terminology.) Asdetailed in the following two sections, we shall use the constraint equations to specify theinitial conditions and to ensure constraint damping and numerical convergence. In testing the robustness of slow contraction, the set of initial conditions under study playsa key role. Analytic-perturbative analyses of smoothing mechanisms can only establish ‘clas-sical smoothing,’ i.e. stability of the attractor FRW solution to small inhomogeneities andanisotropies [8]. To establish ‘robustness,’ the evolution must be studied under a wide setof initial conditions including those that are far outside the perturbative regime of the FRWattractor solution. As we shall describe in this section, our scheme enables the variation ofall freely specifiable variables, such as the initial shear and spatial curvature contributions, { ¯Σ ab , ¯ A a } , as well as the initial field and velocity distributions of the scalar, { φ, ¯ W } . To specify the initial conditions, we first choose a particular time t . With the tetrad andcoordinate gauge conditions remaining the same as specified above for the evolution scheme,the t -hypersurface has constant mean curvature Θ − , the value of which we can freely choose,and zero shift. We note that, using Hubble-normalized variables in our evolution scheme, thenatural units are set by this initial Hubble radius Θ rather than the Planck scale.Next, we must fix the variables (cid:8) ¯ E ai , ¯ n ab , ¯ A a , ¯Σ ab (cid:9) (3.1)that describe the geometry of the t -hypersurface. Not all of these variables are freely specifi-able, though, as they must satisfy the constraint equations (2.60-2.64). (Notably, the evolutionequations (2.53-2.59) propagate the constraints, i.e., ensure that the constraints are satisfiedat later times provided they are satisfied on the initial time slice.)– 12 –dapting the York method [31] commonly used in numerical relativity computations,we choose the initial metric to be conformally flat, g ij ≡ ψ ( x, t ) δ ij , (3.2)where ψ is the conformal factor. The conformal factor is not a free function but fixed by theHamiltonian (or Gauss) constraint (2.60), as we will detail below.With Eq. (2.38), this choice for g ij simultaneously fixes the coordinate components ofthe spatial triad: ¯ E ai = ψ − Θ − δ ai . (3.3)Substituting into the spatial commutator, [¯ e a , ¯ e b ] ≡ ¯ C abc ¯ e c = (cid:16) A [ a δ b ] c + (cid:15) abd ¯ n dc (cid:17) ¯ e c , (3.4)as defined in Eq. (2.7), we find the components of the intrinsic curvature tensor, ¯ n ab = 0 , (3.5) ¯ A a = − ψ − ¯ E ai ∂ i ψ ; (3.6)and the constraint equations (2.62) and (2.64) are trivially satisfied by this combination of ¯ E ai , ¯ n ab , ¯ A a .We stress that, in general, ¯ A a is non-zero, reflecting the fact that the anti-symmetricpart of the intrinsic curvature tensor does not transform trivially under conformal rescaling,as pointed out in Ref. [4]. That means, most especially, assuming the initial metric to beconformally flat does not imply uniform spatial curvature on the initial slice and, hence, doesnot impose a real restriction on the initial data set.Having set the geometric variables { ¯ E ai , ¯ n ab , ¯ A a } through specifying the spatial metricfor the initial slice, it remains to determine the components of the Hubble-normalized sheartensor ¯Σ ab , as defined in Eq. (2.45), to close the set of variables describing the geometry ofthe initial t -hypersurface.It is a particular advantage of the conformal rescaling as suggested by the York methodthat the constraint equations significantly simplify. For this reason, we will determine theinitial shear contribution ¯Σ ab ( x, t ) by first specifying its conformally rescaled counterpart Z ab ( x, t ) ≡ ψ ¯Σ ab ( x, t ) , (3.7)using the momentum constraint (2.61) as evaluated for Z ab ( x, t ) , ¯ E ai ( x, t ) ∂ i Z ab ( x, t ) = Q ( x, t ) ¯ E bi ( x, t ) ∂ i φ ( x, t ) . (3.8)Here, Q ( x, t ) is the conformally rescaled scalar field kinetic energy defined by Q ( x, t ) ≡ ψ ( x, t ) ¯ W ( x, t ) . (3.9)Then, solving the Hamiltonian constraint (2.60) for the conformal factor ψ ( x, t ) yields ¯Σ ab ( x, t ) .It is immediately obvious that we can follow two strategies in specifying Z ab ( x, t ) : eitherwe freely specify the initial shear contribution Z ab ( x, t ) that then fixes parts of the initialfield and velocity distribution { φ ( x, t ) , Q ( x, t ) } , or we freely specify only parts of the initial– 13 –hear contribution Z ab ( x, t ) so we can freely choose { φ ( x, t ) , Q ( x, t ) } that together fix therest of Z ab ( x, t ) using the momentum constraint. In this analysis, we have chosen the latteroption.By definition, the vacuum contribution Z ab ( x, t ) of the conformally rescaled shear ten-sor, i.e. , the divergence-free part of Z ab ( x, t ) , is independent of any matter source. Accord-ingly, we will freely specify only Z ab ( x, t ) and constrain the rest of Z ab ( x, t ) by the choiceof the initial scalar field and velocity distribution { φ ( x, t ) , Q ( x, t ) } .For the numerical simulations, we use periodic boundary conditions ≤ x ≤ π with and π identified. In addition, we restrict ourselves to deviations from homogeneity along asingle spatial direction x so that the spacetimes have two Killing fields. Since the variablesdepend only on x and since x is periodically identified, we specify their spatial variation as asum of Fourier modes with period π .A simple yet generic divergence-free and trace-free choice for Z ab ( x, t ) is given by Z ab ( x, t ) = b ξ ξ a cos x + b a cos x a cos x − b − b − a cos x , (3.10)where ξ, a , a , b and b are constants. Note that since ¯Σ ab must be trace-free, so must Z ab .The rest of the initial shear distribution, Z ab − Z ab is obtained by solving the momentumconstraint (3.8), which turns into an algebraic relation for the Fourier coefficients of thisnon-zero divergence piece of Z ab upon specifying the initial scalar field matter variables. In setting the initial velocity and field distribution, we freely specify the Fourier coefficientsof Q ( x, t ) and φ ( x, t ) via φ ( x, t ) = f cos (cid:0) m x + d (cid:1) , (3.11) Q ( x, t ) = Θ (cid:16) f cos (cid:0) m x + d (cid:1) + Q (cid:17) . (3.12)where f , m , d , f , m , d , and Q are constants.To compute the initial value of the (Hubble-normalized) scalar field gradient ¯ S a ( x, t ) ,we substitute φ ( x, t ) into the constraint equation (2.63). Finally, imposing the Hamiltonian constraint (2.60), our ansatz for the initial set of geometricand matter variables yields an elliptic equation for the conformal factor ψ : ∂ i ∂ i ψ = (cid:0) − − V (cid:1) ψ − (cid:0) ∂ i φ∂ i φ (cid:1) ψ − (cid:16) Q + Z ik Z ik (cid:17) Θ ψ − . (3.13)We numerically solve these equations in the following section. Using the numerical general relativity scheme described in Secs. 2 and 3, we have conductedextensive series of studies to gauge the robustness of slow contraction in smoothing andflattening the universe beginning from a broad spectrum of highly non-perturbative, inhomo-geneous initial matter, spatial curvature and shear distributions. To date, there has been no– 14 –nalogous test of robustness for an expanding cosmology, including inflation; see discussionin Ref. [8].As described in Sec. 3, we solve the full -dimensional Einstein-scalar field equationsbeginning from arbitrary combinations of shear ( Ω s ), spatial curvature ( Ω k ) and matter ( i.e. ,scalar field) density ( Ω m ), where Ω s ≡ ¯Σ ab ¯Σ ab , (4.1) Ω m ≡ ¯ W + ¯ S a ¯ S a + ¯ V , (4.2) Ω k ≡ − ¯ E ai ∂ i ¯ A a + ¯ A a ¯ A a + ¯ n ab ¯ n ab − (¯ n cc ) , (4.3)and (cid:80) i Ω i = 1 . For simplicity, all deviations from homogeneity are along a single spatialdirection x , as in Ref. [8].The scalar field potential energy takes the form V ( φ ) = V exp (cid:0) − √ ε φ (cid:1) (4.4)where V < and ε (cid:29) . The dynamical Ω m = 1 FRW attractor solution corresponds tothe effective equation-of-state ε eff → ε . The initial spatial inhomogeneities in the matterdensity are set by f in the expression for Q ( x, t ) in Eq. (3.12). More precisely, Q ( x, isthe initial velocity distribution, and f is the magnitude of a velocity fluctuation mode withwavenumber m about the mean initial scalar field velocity Q at t = 0 . The sign conventionis that positive Q ( x, t ) corresponds to rolling down the potential (towards more negative valuesof V ( φ ) ). Spatial variations in the initial shear Z ab are set by the parameters a and a inEq. (3.10). Once these parameters are set, the initial conditions for the spatial curvature,shear and matter density are completed by computing the conformal factor ψ in Eq. (3.13)and using Eqs. (3.7) and (3.9), as described in Sec. 3. The studies consisted of a seriesof simulations in which each of these parameters that set the initial conditions was variedindependently from zero to O (1) .Parameters that were found not to affect significantly the robustness tests were heldfixed. For the purposes of the simulations, the coefficient of the potential in Eq. (4.4) wasset to V = − . (in units of the initial Θ ); the initial scalar field was set to φ ( x, t ) = 0 ; theperiod and shift of the sinusoidal spatial variation of scalar field velocity Q ( x, t ) in Eq. (3.12)were set to m = 1 ; d = 0 ; and the remaining constant coefficients in Z ab in Eq. (3.10) wereset to ξ = 0 . , b = − . and b = − . .The result of the numerical studies can be summarized in a series of ‘phase diagrams’indicating the final state as a function of the initial mean scalar field velocity Q and potentialenergy density parameter ε . The different phase diagrams correspond to different types ofinitial conditions as described in the text below. The phase diagrams were made by firstmaking runs for a coarse sampling of Q and ε and characterizing the final states at theend of a representative run lasting N H = 200 e-folds, where N H is equal to the number ofe-folds of contraction of the Hubble radius or, equivalently, Θ = | H | − . Then further runswere made, holding ε fixed and performing a bisection search in Q to identify the boundariesbetween different end states to 1 decimal place precision. By comparing results at differentresolutions, we conclude that the dominant numerical error in our results arises from theprecision used in the bisection search. We note also that smoothing occurs very slowly for thelowest value of ε considered, ε = 4 . , and longer lasting runs suggest that the simulations thathave not smoothed, or only partially smoothed by N H = 200 , will continue to smooth further.– 15 – m = 1 FRW KL completelysmoothed Figure 1 : Phase Diagram I shows the two possible final states ( Ω m = 1 FRW and KL)reached for the case of homogeneous initial conditions as a function of the mean initial scalarfield velocity Q and the scalar field potential parameter ε . The curve divides the diagraminto two regimes.Nevertheless, the diagrams below are based on the end state at N H = 200 for consistency.For ε ≈ . , our results should be viewed as conservative upper bounds for smoothing. Phase Diagram I.
We first consider the case of fully homogeneous initial conditions: f = a = a = 0 . In this case, the simulations converge to one of two homogeneous fixed-pointoutcomes: • an Ω m = 1 FRW universe with Ω k = Ω s = 0 , in which the matter density consists of auniform combination of scalar field kinetic and potential energy (no scalar field gradientenergy density) that has converged to the dynamical attractor solution for all x : namely, ε eff → ˙ φ / ( ˙ φ + V ) → ε (cid:29) ; or, • a “Kasner-like” (KL) universe that is homogeneous, spatially flat ( Ω k = 0 ) and comprisedof some uniform mixture of Ω m and Ω s such that Ω m + Ω s = 1 ; the matter density consistspurely of scalar field kinetic energy density (with zero gradient or potential energy density)corresponding to ε eff → ; the ratio of Ω m to Ω s at the fixed point depends on the initialvalues of the Ω i .The phase diagram in Fig. 1 shows a curve that divides the plot into two regions. (Thecurve is identical to the one marked ∆ f = 0 curve in Fig. 3 of Ref. [8].). All initial conditionsabove the curve converge to the Ω m = 1 FRW fixed point and all initial conditions belowconverge to the KL fixed point. Snapshots of the final distributions of the Ω i ’s for the twotypes of fixed points are shown in Fig. 2.Note that the curve falls below Q = 0 for ε (cid:38) . As discussed in Ref. [8], in cyclic andmost bouncing cosmological models, the natural initial condition at the onset of a contractingphase is that the scalar field is either at rest or evolving down the potential, though withspeeds that can vary with x . This corresponds to mean initial scalar field velocity Q ≥ or‘down the potential.’ Also, as shown in Ref. [8], practical bouncing models require ε (cid:38) inorder that the smoothing and flattening (beginning from non-perturbative initial conditions)is completed rapidly enough that a nearly scale-invariant spectrum of density perturbations– 16 – (cid:127)€(cid:129)€(cid:127)‚‚ Ω (cid:127)€(cid:129)€(cid:127)‚‚ (a) (b) KL Ω m = 1 FRW N (cid:127) N (cid:127) = 120 n (cid:127) N (cid:127) N (cid:127) = 120 Figure 2 : Plots of the final states corresponding to the two phases shown in Phase DiagramI: (a) the dynamical attractor Ω m = 1 FRW state; and (b) the homogeneous Kasner-like (KL)state with non-zero Ω s and Ω m . both final states have Ω k = 0 . Each plot shows normalizedenergy density in matter Ω m (blue), curvature Ω k (red) and shear Ω s (green) as a function of ≤ x ≤ π . N H is equal to the number of e-folds of contraction of Θ = | H | − .consistent with observations can be generated from quantum fluctuations before the bounce.Consequently, the remainder of the discussion will be confined to Q ≥ and, while we willcomment on some interesting behaviors with ε < , the focus, as far as realistic applicationsto cyclic and bouncing cosmology are concerned, should be on the results when ε (cid:38) . InPhase Diagram I, the fixed point in this entire region converges to the desired Ω m = 1 FRWdynamical attractor.
Phase Diagram II:
The second phase diagram (Fig. 3) corresponds to initial conditions inwhich Z ab is homogeneous ( a = a = 0 ) but the initial scalar field velocity distribution is Ω m = 1 FRW mixed: Ω m = FR W + KL un s m oo t h completelysmoothed Δ f = . Figure 3 : Phase Diagram II shows the final states reached for the cases in which the initialscalar field velocity is highly inhomogeneous ( ∆ f = 0 . ) but the Z ab is homogeneous. Overmost of the phase diagram, spacetime is completely smoothed ( Ω m = 1 FRW) or in a mixedstate that is smoothed to an exponential degree (as measured by proper volume).– 17 – Ω x x x x N (cid:127) €(cid:129)€‚ N (cid:127) €(cid:129)€ N (cid:127) €(cid:129)€!" N (cid:127) €(cid:129)€ ‚ (cid:127)€(cid:127)(cid:129)‚(cid:129)(cid:127)(cid:129)‚(cid:129)(cid:127)€€€ (cid:127)€(cid:127)(cid:129)‚(cid:129)(cid:127)(cid:129)‚(cid:129)(cid:127)€€€ Figure 4 : Series of snapshots showing the evolution of the normalized energy density inmatter Ω m (blue), curvature Ω k (red) and shear Ω s (green) for two cases with ∆ f = 0 . drawn from Phase Diagram II: (a) for Q = 0 . , convergence to the Ω m = 1 FRW dynamicalattractor with (cid:15) eff → (cid:15) = 8 ; and (b) for Q = 0 . , convergence to a mixed state in whichspace is divided into segments that are Ω m = 1 FRW separated by segments that are KL butspatially varying. N H is equal to the number of e-folds of contraction of Θ = | H | − .not ( f (cid:54) = 0 ). Here the division between phases depends on ∆ f = f /Q attr , where Q attr is thevalue of Q ( x, t ) for the dynamical attractor solution. ∆ f = O (1) corresponds to large initialvelocity fluctuations in which the initial velocity deviates far from the attractor as a functionof x . We show two bounding curves in the plot (which correspond to two of the examplesshown in Fig. 3 of Ref. [8]). For the practical reasons described above, we restrict the plot tothe regime relevant for cyclic and bouncing cosmology, Q > .In this case, initial conditions corresponding to points above the curve converge to the Ω m = 1 FRW dynamical attractor. Initial conditions below the curve evolve into “mixedstates” in which the volume appears to be divided into segments that converge to the dynam-ical attractor with ε eff = ε (cid:29) separated by segments that are ‘KL but spatially varying,’as shown in Fig. 4. The latter refers to segments in which Ω k = 0 and ε eff = 3 (as in caseof homogeneous KL), but the ratio of Ω m to Ω s varies continuously across the CMC hyper-surfaces. Because the physical volumes corresponding to each segment of length ∆ x contractas a ( τ )∆ x = | τ | /ε eff ∆ x, (4.5)where τ is the proper FRW time coordinate, the slowly contracting Ω m = 1 FRW segmentswith ε eff (cid:29) become exponentially larger than the KL segments with ε eff = 3 as contractionproceeds. The situation is illustrated in Fig. 5 where the main plot is the physical distancedivided by the Hubble radius; the spatially varying KL segment shown in the inset is exponen-tially small by comparison. Similar remarks apply to all the examples of mixed states below.( n.b. This more intuitive argument agrees with the more formal proper volume analysis inRef. [13].) – 18 – (cid:127)€(cid:129)€(cid:127)‚‚
D / | H | -1 x e N H (cid:129) e N H Figure 5 : Sketch of the final (mixed) state shown in Fig. 4(b) in physical distance coordinate D compared to the Hubble radius | H | − showing that the spatially varying KL segment isexponentially small compared to the homogeneous Ω m = 1 FRW region that occupies mostof the spacetime volume.In other words, except for the sliver of small (cid:15) and Q marked ‘unsmooth,’ all initialconditions in the phase diagram either converge to the dynamical attractor solution for all x or into mixed states in which an exponentially large fraction of space time converges to thedynamical attractor but there are also (cosmologically irrelevant) infinitesimal regions withspatially varying KL behavior.Fig. 6 shows state space orbit plots associated with two examples in Fig. 4. The orbit (a) (b) Σ - Σ + - - Σ - Σ + - - Figure 6 : The state space orbits comparing worldlines at x = π for the two models in Fig. 4:(a) converges to the Ω m = 1 FRW dynamical attractor, and (b) evolves into a mixed statethat includes a segment corresponding to spatially varying KL. The point x = π lies in thespatially varying KL region in the second example; the corresponding orbit never reaches thecenter corresponding to the Ω m = 1 FRW dynamical attractor.– 19 – m = 1 FRW completelysmoothed Ω m = 1 FRW completelysmoothed m i xe d un s m oo t h e d se r i es o f n ea r l y K L Ω m = F R W a = 0.8, a = f = 0 a = 0.8, a = f = 0 m i xe d un s m oo t h e d Figure 7 : Phase Diagram III (left) shows the final states reached for the cases in which theinitial scalar field velocity is homogeneous f = 0 but the off-diagonal components of Z ab areinhomogeneous ( a = 0 . but a = 0 ). For Phase Diagram IV (right), only the on-diagonalcomponents are highly inhomogeneous ( a = 0 . but a = 0 ).plots enable the visualization of the evolution of the shear at a chosen point x in the ( ¯Σ + , ¯Σ − ) plane, where ¯Σ + = (cid:16) ¯Σ + ¯Σ (cid:17) , ¯Σ − = √ (cid:16) ¯Σ − ¯Σ (cid:17) . (4.6)The ¯Σ ± are normalized so that the unit circle ( ¯Σ + ¯Σ − = 1 ) corresponds to the vacuumKasner solution; trajectories that approach Ω m = 1 FRW converge to the center. See Refs. [8,13] for other examples. In the case of Fig. 6, the orbit converges to the center ( Ω m = 1 FRW)for the case in Fig. 4a for all x ; for the mixed state case in Fig. 4b, the point x = π lies inthe spatially varying KL region between the Kasner circle and the origin. Phase Diagrams III and IV:
The third phase diagram (Fig. 7, left) corresponds to initialconditions in which the only initial inhomogeneity is due to the spatially-varying off-diagonalcomponents of Z ab in Eq. (3.10) (the terms proportional to a ). The initial scalar field velocitydistribution is homogeneous ( f = 0 ) and the diagonal components of Z ab proportional to a are also zero. Here we find that, as in Phase Diagram I (Fig. 1), the entire range of Q ≥ and ε (cid:38) converges directly to the Ω m = 1 FRW dynamical attractor. That is, unlikethe previous example, introducing inhomogeneity by making a non-zero does not reduce therobust smoothing and flattening behavior over the range of parameters relevant for cyclic andbouncing models.In fact, a curious feature occurs for examples in the slivers of the phase diagram where ε (cid:46) and Q is near zero. This is the range that leads to homogeneous KL behavior in PhaseDiagram I. When a is non-zero and (cid:46) O (0 . , we find instead that, beginning from a nearlyhomogeneous state, the evolution first approaches a nearly homogeneous KL fixed point butthen veers away and goes through a series of further nearly KL stages before converging on a Ω m = 1 FRW dynamical attractor solution. This is illustrated in the upper panels shown inFig. 8 and by the state space orbit diagram in Fig. 9 (left). Effectively, this means that non-zero a actually enlarges the phase space region that converges to the attractor solution. As a is increased by another order of magnitude, the initial state is significantly inhomogeneous– 20 – Ω x x x x x (cid:127)€(cid:127)(cid:129)‚(cid:129)(cid:127)(cid:129)‚(cid:129)(cid:127)€€€ (cid:127)€(cid:127)(cid:129)‚(cid:129)(cid:127)(cid:129)‚(cid:129)(cid:127)€€€ (cid:127)€(cid:129)(cid:127)‚(cid:129) N !"! N !"!$% N !"!& N !"!' N !"!%% N !"! N !"!$% N !"!(( N !"!)** N !"!)( Figure 8 : Series of snapshots showing the evolution of the normalized energy density inmatter Ω m (blue), curvature Ω k (red) and shear Ω s (green) for ≤ x ≤ π for two caseswith (cid:15) = 13 and Q = 0 in which the inhomogeneity is solely an off-diagonal inhomogeneouscontribution to Z ab : (a) a = 0 . and f = a = 0 ; and (b) a = 1 . and f = a = 0 . Thefirst evolves through a series of nearly homogeneous KL stages, then a series of inhomogeneousand mixed stages before converging to a final homogenous Ω m = 1 FRW dynamical attractor;the second case begins significantly inhomogeneous (large a ) and passes through a longsequence of different mixed states that include spikes (such as the one visible near the centerof the lower panel labeled N H = 55 ) before reaching the final homogenous Ω m = 1 FRWdynamical attractor.compared to the small a case; nevertheless, after evolving through a different sequence ofstages, it also converges to the homogeneous Ω m = 1 FRW dynamical attractor, as illustratedin the lower panels shown in Fig. 8 space orbit diagram in Fig. 9 (left). (In the slivers of thephase diagram labeled as undergoing a series of nearly KL stages or as ending with mixedregions, the numerical evolution includes some spikes of negligible extent in x , as describedin Ref. [13]. An example is visible in the middle of the snapshot labeled N H = 55 in thelower row of Fig. 8; although not apparent in the subsequent snaphots, its presence can beidentified in higher resolution runs.)The fourth phase diagram (Fig. 7, right) is the complementary case where the non-zero inhomogeneous components of Z ab are on the diagonal ( a (cid:54) = 0 ) and the off-diagonalcomponents are zero ( a = 0 ). As in Phase Diagram II, the entire range of Q ≥ and ε (cid:38) converges directly to the Ω m = 1 FRW dynamical attractor. For the corner of thephase diagram where a = O (1) , as in Fig. 7 (right) and ε (cid:46) , the result is a mixed state inwhich an exponentially large fraction of space time converges to the FRW dynamical attractorbut there are also (cosmologically irrelevant) infinitesimal regions with spatially varying KLbehavior. For yet smaller ε , the spacetime is not smoothed.– 21 – a) (b) Σ - Σ + Σ + - - Σ - - - Figure 9 : The state space orbits comparing worldlines at x = π for the two models in Fig. 8:(a) one that passes through a series of nearly homogeneous KL stages but keeps veering awayuntil finally converging to the Ω m = 1 FRW dynamical attractor; and (b) one that evolvesthrough a different sequences of mixed states before converging to the attractor. The point x = π lies in a region that goes through multiple stages in both cases. Phase Diagrams V:
Finally, we present a simplified phase diagram corresponding to casesin which any combination of inhomogeneities (non-zero f , a and a up to O (1) ) is con-sidered (Fig. 10). This diagram encapsulates the take-away lesson regarding the remarkablerobustness of slow contraction: The entire regime of practical relevance to cyclic and bouncingcosmology models – that is, having ε (cid:38) (the condition for sufficiently rapid smoothing) a = 0.5 & a = 0.5 Ω m = 1 FRW completelysmoothed mixed: Ω m = F R W + K L un s m oo t h e d Δ f = . completelysmoothed (modulo spikes) Figure 10 : Phase Diagram V shows the final states reached when a combination of inhomo-geneities is considered (that is, is f , a and a are all non-zero). The entire region relevantto cyclic and bouncing cosmology models – ε (cid:38) and Q ≥ converges completely or to anexponential degree (as measured by proper volume) to the desired smooth, anisotropic andflat dynamical attractor solution. – 22 – igure 11 : Two types of potentials V ( φ ) (solid curves) used to test the robustness of slowcontraction in smoothing and flattening spacetime. The scalar field evolves from right to leftin each case (see arrow), first encountering a finite range of V ( φ ) that closely matches thenegative exponential potential (dotted lines) considered in the phase diagrams above with ε = 50 . From that point onward (after just a few e -folds of contraction), the potentialsdeviate significantly: (a) the negative potential reaches a minimum and then approaches zerofrom below; and (b) the negative potential reaches a minimum and rises rapidly above zerowith a shape that is super-exponential ( e.g. , φ exp ( αφ ) for some constant α ).and initial mean scalar field velocity at rest or evolving down the potential V ( φ ) for all x )– converges completely or to an exponential degree (as measured by proper volume) to the desiredsmooth and flat FRW dynamical attractor solution with Ω m = 1 and Ω s = Ω k = 0 . (In theAppendix A, we describe our tests for numerical convergence and demonstrate that the testsare well-satisfied in those regions of the phase diagram that completely smooth. The sameis found to be true in the regions of spacetime that directly smooth without spikes in thecase of mixed final states. More complex results are found in regions with spikey behavior.)In this and in all the diagrams above, we allow highly non-perturbative initial conditions,but we restrict them to be less than or O (1) , which is at a level at or beyond what wouldbe naturally expected as plausible initial conditions. Pushing parameters beyond the valuesconsidered here will move boundaries to slightly higher values of ε , say; but there alwaysremains a substantial range of the phase diagram that converges to the desired smooth andflat FRW dynamical attractor solution with Ω m = 1 and Ω s = Ω k = 0 . In particular, modelswith ε > (or M < m Pl / in Eq. (1.3) where m Pl is the reduced Planck mass) are plausi-ble in microphysical models, and they are exponentially more powerful and rapid in drawingmuch broader ranges of initial conditions to the dynamical attractor solution compared tothe already-robust cases considered here.As a further test of robustness, a series of simulations were performed in which the scalarfield potential V ( φ ) significantly deviates away from the simple negative potential (Eq. (4.4)that was assumed in the phase diagrams above. The plots of V ( φ ) in Fig. 11 illustrate twovariations that were explored. In each case the scalar field evolves from right to left beginningwith a segment that closely matches the negative exponential potential Eq. (4.4). During thisfirst stage, there is a short period of slow contraction with effective equation-of-state ε eff → ε ,during which the large initial inhomogeneities are substantially smoothed and flattened asabove. However, as φ → −∞ (that is, evolving to the left), the potentials diverge and ε eff decreases rapidly and significantly. – 23 – Ω x x x x N (cid:127) €(cid:129)€‚ N (cid:127) €(cid:129)€ !‚ N (cid:127) €(cid:129)€"! N (cid:127) €(cid:129)€$‚‚ (cid:127)€(cid:129)€(cid:127)€(cid:129)€(cid:127)€(cid:129)€(cid:127)€(cid:129)€ %&'%(' ε €(cid:129)€)!" ε €(cid:129)€ !$ Figure 12 : Series of snapshots showing the evolution of the normalized energy density inmatter Ω m (blue), curvature Ω k (red) and shear Ω s (green) for ≤ x ≤ π for the two typesof potentials described in Fig. 11. In the top example, ε eff → as φ → −∞ ; in the bottomexample, ε eff falls below three as the φ climbs the steeply rising positive part of the potential.In either case, the smoothing and flattening created by the initial slow contraction phase ismaintained. Σ + Σ - - - Figure 13 : The state space orbits comparing worldlines at x = 3 π/ for the two models inFig. 12: (a) as φ → −∞ , V ( φ ) approaches zero from below (red, solid); and (b) as φ → −∞ , V ( φ ) rises above zero and increases super-exponentially (blue, dotted). In both cases, theorbits converge rapidly to the center during the early slow contraction phase and remain therewhen the potential sharply deviates from negative exponential and the slow contraction phaseends. – 24 – igure 14 : A plot showing the effective equation-of-state ε eff ( x ) at a sequence of times forthe case shown in Fig. 11b in which V ( φ ) rises above zero and increases super-exponentially.For early times when V < and exponentially decreasing ε eff rises to a large positive value(red dashed curves for which the number of e-folds of contraction are n H = 3 . and . ).At later times when V > and super-exponentially decreasing, the spacetime and ε becomenearly uniform and ε eff decreases and falls below ε = 3 (sequence of green thin solid curveswith progressively increasing values of n H ). After reaching a minimum at n H = 7 . , ε eff increases and approaches the homogeneous FRW fixed point with ε = 3 (blue solid thickcurve).In the first case, the negative potential reaches a minimum and then rises to approachzero as φ → −∞ . The equation-of-state ε eff → . Despite this significant change, the earlystage of slow contraction is rapid and powerful enough that the spacetime remains smoothand flat, as shown in the top panel of Fig. 12 and by the red solid trajectory in the statespace orbit plot shown in Fig. 13.In the second case, the negative potential reaches a minimum and then rises super-exponentially rapidly above zero; e.g., V ( φ ) ∝ φ exp ( αφ ) as φ → −∞ . The super-exponentialbehavior is artificially introduced to force the equation-of-state to reach ε eff < , as shownin Fig. 14. Note that the potential is not well-motivated physically and does not correspondto any bouncing or cyclic picture; rather, it is specifically introduced as an extreme test ofrobustness. Even in this case, the early stage of slow contraction is rapid and powerful enoughthat the spacetime remains smooth and flat, as shown in the bottom panel of Fig. 12 and bythe blue dotted trajectory in the state space orbit plot shown in Fig. 13.We have constructed examples of even more steeply rising, super-exponential potentialsin which ε eff not only falls below three, but also well below zero. In this case, some smalldeviations from perfect homogeneity appear for some ranges of x (a kind of mixed state)as the field climbs up the steep positive wall, though even then smoothness and flatness isretained for most of the range of x . – 25 – Analytic approximation
The key result of our numerical analysis is that slow contraction is an astonishingly robustsmoothing mechanism: for most initial conditions, the scalar field energy density rapidlyhomogenizes and quickly becomes the dominant component ( Ω m ≡ ), driving the geometryto a spatially-flat, homogeneous, and isotropic (FRW) spacetime everywhere and leading to complete smoothing for all x . For extreme initial conditions that do not result in a completelysmoothed FRW spacetime, the end state is either a mixed that is smooth almost everywhere,as measured by co-moving volume, or a homogeneous but anisotropic spacetime described by a‘Kasner-like’ (KL) metric. To conclude our study, we complement our numerical computationwith an analysis in which we demonstrate that the homogeneous end states found numericallycorrespond to the only possible stable critical points for a contracting universe. Notably, in all the cases we studied, the evolution towards the homogeneous end states foundnumerically approach the so-called ultra-local limit, in which all terms that involve gradientsbecome negligible in the evolution and constraint equations, i.e. , ¯ E ai → , ¯ A a → , ¯ S a → , (5.1)as the universe contracts ( t → −∞ ). Most remarkably, ultra-locality is generically reachedrapidly even in those cases where the initial conditions include large gradients, as is apparentfrom the plot of Ω k in the N H = 0 panels shown in Figs. 4, 8 (lower panel) and 12. Notethat, by choosing the spatial metric on the initial t -hypersurface to be conformally flat, Ω k is a direct measure of spacetime gradients.In the ultra-local limit, the evolution and constraint equations (2.53, 2.56, 2.59 and2.62-2.64) for ¯ E ai , ¯ A a and ¯ S a , respectively, are automatically satisfied, while the remainingsystem of non-trivial evolution and constraint equations dramatically simplifies, leaving asystem of first-order ordinary differential equations (ODEs) for twelve metric and two mattervariables constrained by three algebraic relations, as opposed to the original set of coupledpartial differential equations with twenty-four metric and five matter variables.Most importantly, in the ultra-local limit the momentum constraint (2.61), (cid:15) abc ¯ n bd ¯Σ cd = 0 , (5.2)is equivalent to demanding that the shear tensor ¯Σ ab and (the symmetric part of) the intrinsiccurvature tensor ¯ n ab commute, which means that ¯Σ ab and ¯ n ab share the same eigenvectors .Furthermore, as shown in the Appendix B, the eigenvectors at time t remain eigenvectorsat later times. The dynamics in the ultra-local limit, therefore, is completely determined bythe evolution of the eigenvalues of ¯Σ ab and ¯ n ab .As a consequence, the twelve coupled ODEs for the metric variables ¯Σ ab and ¯ n ab , ˙¯Σ ab = − (cid:16) N − (cid:17) ¯Σ ab − N (cid:16) n c (cid:104) a ¯ n b (cid:105) c − ¯ n cc ¯ n (cid:104) ab (cid:105) (cid:17) , (5.3) ˙¯ n ab = − (cid:16) N − (cid:17) ¯ n ab + 2 N ¯ n c ( a ¯Σ b ) c , (5.4)can be replaced by six ODEs for the corresponding eigenvalues σ i , ν i ( i = 1 , , ), such thatthe dynamical system is fully described by the following set of six evolution equations ˙ σ = (cid:16) − N (cid:17) σ − N (cid:16)(cid:0) ν − ν − ν (cid:1) ν − (cid:0) ν − ν (cid:1) (cid:17) , (5.5)– 26 – σ = (cid:16) − N (cid:17) σ − N (cid:16)(cid:0) ν − ν − ν (cid:1) ν − (cid:0) ν − ν (cid:1) (cid:17) , (5.6) ˙ ν = (cid:16) N (cid:0) σ − (cid:1)(cid:17) ν , (5.7) ˙ ν = (cid:16) N (cid:0) σ − (cid:1)(cid:17) ν , (5.8) ˙ ν = (cid:16) − N (cid:0) σ + 2 σ + 1 (cid:1)(cid:17) ν , (5.9) ˙¯ W = − (cid:16) N − (cid:17) ¯ W + √ ε N ¯ V , (5.10)subject to the Hamiltonian constraint N − − (cid:16) ν + ν + ν (cid:17) + (cid:16) ν + ν + ν (cid:17) − (cid:16) σ + σ σ + σ (cid:17) − ¯ W = 0 , (5.11)and the ultra-local limit of the Hubble-normalized lapse equation (2.52), N − = 1 + (cid:16) σ + σ σ + σ (cid:17) + (cid:16) ¯ W − ¯ V ( φ ) (cid:17) . (5.12)As before, dot denotes differentiation with respect to (coordinate) time t , as defined inEq. (2.47). The third shear eigenvalue σ could be eliminated since the trace-freeness of theshear tensor implies that the three eigenvalues must sum to zero. In addition, we substitutedin Eq. (5.10) ¯ V ,φ ¯ V = −√ ε, (5.13)to eliminate the term proportional to ¯ V ,φ in Eq. (2.58) (and, hence, eliminate any explicit φ -dependence from the equations above). As detailed above, our goal is to show that the homogenous end states we identified numer-ically are the stable attractor solutions of the underlying system of evolution and constraintequations above. Since any attractor is a (stationary) critical point, we begin with identifyingall critical points of the system as given by Eqs. (5.5-5.10). We list the complete set of (seven)critical point solutions of Eqs. (5.5-5.10) in Table 1.If the potential energy density V ( φ ) is negative, as needed to drive slow contraction (seeSec. 1 and Refs. [17, 18, 20, 26]), there exists only the three distinct critical point solutionslisted in the first three lines of Table 1:- flat, homogeneous, and isotropic (FRW) scaling solution:- flat, homogeneous, and anisotropic (Kasner-like with Ω m ≥ ) solution; and- curved, homogeneous, and anisotropic ( Ω m = 0 ) solution.Since the remaining four critical points listed in rows 4-7 of Table 1 only exists for positivepotentials, they cannot trigger a phase of slow contraction and lie outside of the scope of thispaper. For this reason, we do not further consider them here.– 27 – σ σ ν ν ν ¯ W ¯ V ε ε > − ε (cid:54) = 0 (cid:54) = 0 − ( σ + σ σ + σ ) − − (cid:54) = 0 ν ± (cid:113) ε − ν ν ε > W (cid:54) = 0 ν ε = ¯ W ε +89 ε − εε +8 2 ε − ε +8 ± √ (1 − ε )( ε − ε +8 < ε ( ε +8) < · (4 − ε )( ε +8) ε +23 ε < − ε ε ε − ε ± √ ε − ε +2 − ν ε (2+ ε ) ε (2+ ε ) Table 1 : Critical point solutions corresponding to the autonomous system of ordinary dif-ferential equations (5.5-5.10) describing the evolution in the ultra-local limit.
A critical point is an attractor solution if it is stable to small perturbations. Otherwise, it isa repeller.
Accordingly, to decide which of the three critical points are attractor solutions,we linearize the system around each critical point and solve the resulting constant-coefficientsystem for the complete set of perturbation variables. (Here, we only show the three sets ofequations corresponding to the critical points. For the perturbed evolution and constraintequations around an arbitrary background, see Appendix C.) ε > First, we linearize the system around the homogeneous and isotropic FRW critical point σ i = ν i = 0 (for all i ); N − = ¯ W = ε ; ¯ V = 3 − ε. (5.14)The perturbed shear, spatial curvature and scalar-field matter decouple and obey thefollowing evolution equations: δ ˙ σ i = (cid:18) − ε (cid:19) δσ i , i = 1 , (5.15) δ ˙ ν j = (cid:18) − ε (cid:19) δν j , j = 1 , , (5.16) δ ˙¯ W = (cid:18) − ε (cid:19) δ ¯ W . (5.17)Solutions to the system admit a simple exponential scaling behavior, δσ i , δ ¯ W ∝ e (cid:0) − ε (cid:1) t , δν j ∝ e (cid:0) − ε (cid:1) t . (5.18)– 28 –t is immediately apparent that, for a contracting spacetime ( t → −∞ ), the FRW scalingsolution is a stable attractor for ε > and a repeller otherwise.Notably, the FRW critical point solution recovers the well-known scaling attractor solu-tion, a = ( τ /τ ) /ε , φ = (cid:112) /ε ln( τ /τ ) , (5.19)(where τ is the proper FRW time), as often cited in the cosmology literature. In particular,the eigenvalues correspond to the so-called ‘Friedmann variables’ commonly used to identifythe scaling solution while assuming an FRW background. (For details, see the Appendix D.) Linearizing Eqs. (5.5-5.10) for the homogeneous and anisotropic Kasner-like critical pointsolution, N = , σ , σ (cid:54) = 0 , ν , ν , ν = 0 , ¯ W = 3 − ( σ + σ σ + σ ) , ¯ V = 0 , (5.20)the system reduces to two simple decoupled constant-coefficient matrix equations: one for thethree spatial curvature eigenvalue variables, δ ˙ ν δ ˙ ν δ ˙ ν = 23 σ + 1 0 00 σ + 1 00 0 1 − σ − σ δν δν δν , (5.21)and one for the shear eigenvalue and scalar-field matter variables, δ ˙ σ δ ˙ σ δ ˙¯ W = J − p J δσ δσ δ ¯ W , (5.22)where p ≡ (cid:18) − √ ε · (cid:113) − ( σ + σ σ + σ ) (cid:19) , (5.23)and the matrix J is given by J = − ¯ W σ + σ − σ +2 σ σ + σ σ ¯ W −√ ε σ ¯ W −√ ε . (5.24)As before, solutions to the system admit a simple exponential scaling behavior: δν , ∝ e (cid:0) σ , (cid:1) t , δν ∝ e (cid:0) − σ − σ (cid:1) t , (5.25) δσ , , δ ¯ W ∝ e (cid:0) −√ ε · √ − ( σ + σ σ + σ ) (cid:1) t . (5.26)Intriguingly, though, the Kasner-like solution can behave both as an attractor or a repellerin a contracting universe ( t → −∞ ): if ε = 0 and σ i > − for i = 1 , , ; or if ε > , σ i > − and ( σ + σ + σ ) > − /ε , the Kasner-like solution is an attractor. Otherwise, it is arepeller. Examples for both behaviors are shown in Figure 15.– 29 – a) (b) Σ - Σ + Σ + - - Σ - - - Figure 15 : State orbit plots for two cases: (a) starting with initial conditions ¯ W = 0 , { σ , σ , σ } = {− . , − . , . } and { ν , ν , ν } = { . , . , } and ε = 0 in Eq. (5.10),spacetime undergoes a series of bounces before converging to a homogeneous Kasner at-tractor solution with ¯ W = 0 , { σ , σ , σ } = { . , − . , − . } and { ν , ν , ν } = { , , } .and (b) starting with initial conditions ¯ W = 0 , { σ , σ , σ } = {− . , − . , } and { ν , ν , ν } = { . , . , } and ε = 6 , spacetime undergoes a different sequence of bouncesbefore converging to a homogeneous and isotropic (spatially-flat) FRW attractor solutionwith ¯ W = √ , { σ , σ , σ } = { , , } and { ν , ν , ν } = { , , } .In Case (a) on the left, the system converges to a Kasner-like state, while in Case (b)on the right, the system is driven to the FRW scaling attractor solution. This second caseis especially important since it makes apparent the difference to the well-known vacuum case(that is, pure gravity with no scalar field). In the pure vacuum case, the Kasner solution isthe only stable attractor. In the presence of a scalar-field, though, reaching the Kasner-likeattractor is only possible under very special initial conditions, namely, when the initial scalarfield velocity is uphill (that is, in the direction that V ( φ ) increases). In this case, the scalarfield’s relative contribution to the total energy density ( Ω m ) becomes negligible and there isthe same dynamical behavior as in the vacuum case. But for cosmologically motivated initialconditions, as discussed in Ref. [8] and the previous section, the initial scalar field velocityis at rest or downhill; then one finds that the scalar field energy density increases relative tothe other components and the FRW scaling solution becomes the only attractor.An illustrative example for how large the set of initial data is that belongs to the basinof attraction of the FRW solution was given above in Figure 8. There the initial scalar fieldmatter energy density was chosen so small that the system first approached a Kasner-likecritical point, just as it would in the absence of matter. Then the system was dynamicallydriven away by the growing homogeneous spatial curvature to another Kasner-like criticalpoints through a series of (mixmaster) bounces. Yet, as the scalar field’s energy density con-tinued to grow relative to other contributions, all Kasner-like critical points became repellersand the system eventually settled in the FRW attractor solution.Finally, we note that this analysis also complements and generalizes our numerical studyin that it enables us to follow the evolution under homogeneous but spatially curved ( ν j (cid:54) = 0 for at least one j ∈ { , , } ) initial data. This is the one type of initial data that is excluded– 30 –n our numerical study by assuming a conformally-flat spatial metric on the initial spacelike t -hypersurface. Thirdly and lastly, we linearize the system (5.5-5.10) around the curved homogeneous andanisotropic critical point solution, N = , σ = σ = − , ν = ν ≡ ν (cid:54) = 0 , ν = 0 , ¯ W = 0 , ¯ V = 0 , (5.27)leading to a simple set of evolution equations: δ ˙ σ δ ˙ σ δ ˙ ν δ ˙ ν δ ˙ ν δ ˙¯ W = − ν ν ν − ν − ν − ν − ν − ν − ν − ν
00 0 0 0 2 0 √ ε √ ε √ ε ν · δσ δσ δν δν δν δ ¯ W . (5.28)For some initial perturbations δσ i , δν j , δ ¯ W , the system admits the following solutions: δσ = 12 (cid:16) δσ − δσ + (cid:16) δσ + δσ (cid:17) e t (cid:17) , (5.29) δσ = 12 (cid:16) δσ − δσ + (cid:16) δσ + δσ (cid:17) e t (cid:17) , (5.30) δν = δν + ν (cid:16) δν − δν (cid:17) t + ν (cid:16) δν + δν + ν δν (cid:17) (cid:0) − e t (cid:1) , (5.31) δν = δν − ν (cid:16) δν − δν (cid:17) t + ν (cid:16) δν + δν + ν δν (cid:17) (cid:0) − e t (cid:1) , (5.32) δν = δν e t , (5.33) δ ¯ W = δ ¯ W − (cid:114) ε (cid:16) δν + δν + ν δν (cid:17) (cid:0) − e t (cid:1) . (5.34)Unlike before, as contraction proceeds, the dominant scaling behavior is not the expo-nential but the constant or power-law terms in the solutions. More exactly, as t → −∞ , theshear and scalar field fluctuations converge to constants δσ → (cid:0) δσ − δσ (cid:1) , δσ → (cid:0) δσ − δσ (cid:1) , (5.35) δ ¯ W → δ ¯ W − (cid:113) ε (cid:16) δν + δν + ν δν (cid:17) , (5.36)while two of the spatial curvature perturbations grow δν → ν (cid:16) δν − δν (cid:17) t, δν → − ν (cid:16) δν − δν (cid:17) t, (5.37)driving the system away from the critical point, which is therefore unstable .– 31 – Summary and discussion
The combination of numerical relativity simulations in Sec. 4 and the critical-point analysisin the ultra-local limit in Sec. 5 provide complementary information about the power of slowcontraction to smooth and flatten spacetime.The numerical studies show the robustness of slow contraction in transforming space-times with wildly non-perturbative inhomogeneous initial conditions into homogeneous space-times that approach the ultra-local limit. The analytic studies prove that the only possibleultra-local end states are either Ω m = 1 FRW or a Kasner-like universe with a combinationof anisotropy and matter energy density.Which end point is reached depends on the initial conditions. In a series of phasediagrams exploring different combinations of shear and intrinsic curvature inhomogeneities, wehave shown that for negative exponential potentials with ( V ,φ /V ) ≡ ε (cid:38) and physicallyplausible initial scalar field velocity distributions (that is, at rest or downhill for all x ), theuniverse is driven to an Ω m = 1 FRW spacetime with ε eff = ε , as required in bouncing andcyclic models of the universe. Similar results were found for more complicated potentials forwhich the condition on V ,φ /V is only maintained for a short interval of time.The phase diagrams also show that, in some cases, an initial inhomogeneity can favoreventual convergence to the Ω m = 1 FRW spacetime. For example, there are regions of PhaseDiagram III that would converge to Kasner if there is no initial inhomogeneity ( f = a = a = 0 ), but that are instead driven, after a few bounces, to FRW when a is set to an evenrelatively small value, such as a = 0 . in Fig. 8(a). This suggests that the basin of attractionfor the KL critical point is small. In other cases, the result is a mixed state that is almostentirely FRW (as measured by proper volume) interspersed with exponentially tiny Kasner-like regions over which the ratio of matter energy density and shear vary. The ultimate fateof these comparatively infinitesimal Kasner-like regions is an interesting academic questionthat is not yet resolved, but one that is not of practical relevance to cosmology because oftheir insignificant volume weight.The bottom-line of the complementary studies is that slow contraction is an even morerobust smoothing and flattening mechanism than imagined based on earlier heuristic argu-ments or than has been shown for any other proposed cosmological smoothing and flatteningprocess. Acknowledgments
We thank David Garfinkle for helpful comments and discussions. The work of A.I. is sup-ported by the Lise Meitner Excellence Program of the Max Planck Society and by the SimonsFoundation grant number 663083. W.G.C. is partially supported by the Simons Foundationgrant number 654561. F.P. acknowledges support from NSF grant PHY-1912171, the SimonsFoundation, and the Canadian Institute For Advanced Research (CIFAR). P.J.S. is supportedin part by the DOE grant number DEFG02-91ER40671 and by the Simons Foundation grantnumber 654561. – 32 –
Numerical methods and convergence tests
In this Appendix, we describe our tests for numerical convergence. The key result is thatcases of interest to bouncing cosmology – those regions of the phase diagram that completelysmooth and converge to the Ω m = 1 FRW dynamical attractor solution – strongly satisfyall tests. The same is found to hold in cases ending with mixed states for those regions ofspacetime that smooth without encountering spikes; these regions occupy almost the entirespacetime, as measured by proper volume. The effects on convergence are also shown for theexponentially small non-smooth regions that undergo spikey behavior.To numerically solve the system of equations detailed in Eqs. (2.52-2.59), we use 2ndorder accurate spatial derivatives, and a 3 step method for time integration given by the Iter-ated Crank-Nicolson method. The evolution equations consist of a coupled elliptic-hyperbolicsystem of equations, so at each sub-step of the time integration, we first solve the elliptic equa-tion for the Hubble-normalized lapse N through a relaxation method, and then update thehyperbolic equations to the next Iterated Crank-Nicolson sub-step. In the simulations demon-strated above we use a grid of 1024 points, with ∆ x = 2 π/ , with a Courant factor of0.5. To demonstrate the convergence of our code, the error and convergence was analyzedfor a broad range of examples drawn from Phase Diagram V (Fig. 10). Here we present theresults for two representative cases: c = 5 , a = 0 . , a = 0 . , f /c = 0 . , Q = 2 . (A.1)corresponding to a simulation that smooths everywhere to FRW without encountering spikes,and c = 5 , a = 0 . , a = 0 . , f /c = 0 . , Q = 1 . (A.2)corresponding to a simulation that ends in a mixed state with spikes.Fig. 16 shows the L2 norm of the Hamiltonian constraint (Eq. 2.60) integrated over thespatial domain as a function of time. We see that, in the case of a spacetime that smoothsdirectly to FRW everywhere, after an initial period of second order convergence, the constraintconverges faster than expected, at third order (Fig. 16 left). In the case where a region ofthe spacetime does not smooth, we see a reduction in convergence at these points, droppingto second order or worse when spikes form. Between the times of spike formation we retain3rd order convergence (Fig. 16 right).For the second case, Fig. 17 shows the modulus of the Hamiltonian constraint as afunction of the spatial coordinate x at a fixed time (in this case, n H = 99 ). Where spikes form,we see that the Hamiltonian constraint in the regions where the spikes form ( . (cid:47) x (cid:47) . )do not converge at the required rate; but in the outer regions, where the spacetime hassmoothed to FRW without encountering spikes, the convergence properties remain unaffected.For the same mixed state example, Fig. 18 tracks the evolution of Σ + in Eq. (4.6) attwo fixed spatial points. The left panel focuses on a point that smooths to FRW, and theright focuses on a point that remains unsmoothed and encounters spikes. In the region thatsmooths to FRW, we see second order convergence for this variable as expected. Assumingsecond order convergence and using the data with the production resolution of 1024 gridpoints (medium resolution subscript m ) and data with double the resolution (fine resolutionsubscript f ), we perform a Richardson extrapolation, which is defined as R = 4(Σ + ,f − Σ + ,m / . (A.3)– 33 – ixed final state with spiking S ca l e I n t e g r a t e d H a m il t on i a n C on s t r a i n t V i o l a t i on Complete smoothing to FRWwith no spiking coarse/64medium/8fine n H n H coarse/64medium/8fine
Figure 16 : The Hamiltonian constraint integrated over the spatial domain as a functionof time, for 3 resolutions, coarse ∆ x c = 2 π/ , medium ∆ x m = 2 π/ and fine ∆ x f =2 π/ . Left , parameters of Eq. (A.1): the evolution encounters no spikes while smoothingto FRW. We plot the integrated Hamiltonian constraint rescaled by the appropriate factorfor 3rd order convergence ( || H || c / , || H || m / , || H || f ) . Right , parameters of Eq. (A.2): theevolution does form spikes at this point. We plot the integrated Hamiltonian constraintrescaled by the appropriate factor for 2nd order convergence ( || H || c / , || H || m / , || H || f ) . H a m il t on i a n C on s t r a i n t V i o l a t i on -8 -10 -12 -14 -16 x / 2 π n H = 99 Figure 17 : The Hamiltonian constraint as a function of space at fixed times for the same 3resolutions (with the same color coding) described in Fig. 16 (right), rescaled by factors todemonstrate 3rd order convergence, for parameters of Eq. (A.2) that lead to a mixed final statecontaining an ‘inner region’ that does not smooth to FRW surrounded by a smooth regionthat does. Recall that the inner region is exponentially small compared to the smoothedregion when measured by proper volume. In the regions that smooth to FRW, the thirdorder convergence is unaffected by the presence of the spikes in the inner region. By contrast,convergence is lost in the inner region. – 34 – moothing with spikes Smoothing to Ω m =1 FRW (no spikes) n H n H x -7 For x if complete smoothing (no spikes)or for x in smooth regions away from spikes Example of x that undergoes spikes D i ff e r e n ce b e t w ee n R i c h a r d s on ex t r a po l a t i on & M e d i u m R es o l u t i on Figure 18 : The difference between the Richardson extrapolation of the trajectory of Σ + and the calculated value at resolution ∆ x = 2 π/ for spatial points x = 3 π/ (left) and x = π (right), with parameters of Eq. (A.2). For a spatial point that smooths directly toFRW the error is consistently small, ∼ − (left). For a point that does not smooth andwhich remains in a Kasner-like region with spikes, the error remains small until spikes form,at which point the error grows in magnitude (right). The growth in error corresponds to atransition between Kasner states at x = π .We take the error as the difference between this extrapolation and the medium resolutionresults. We estimate the absolute size of the error as approximately − in the trajectoryof Σ + that smooths to FRW (Fig. 18 left). For the trajectory not smoothing to FRW we seethat the variable still converges at 2nd order until the formation of spikes, at which pointerrors can grow as large as − (Fig. 18 right). B Evolution of the eigensystem in the ultra-local limit
Here, we show that, in the ultra-local limit ( ¯ E ai → , ¯ A a → , ¯ S a → ), the eigenvectors ofthe shear and intrinsic curvature tensors, ¯Σ ab and ¯ n ab , at some time t remain eigenvectorsat later times. As a result, all the dynamics is incapsulated in the evolution of the shear andspatial curvature eigenvalues σ i , ν i ( i = 1 , , ). Of course, this is a trivial statement if theshear and intrinsic curvature tensors are diagonal. Here, we generalize to the case that ¯Σ ab and ¯ n ab have non-zero off-diagonal components.Since both ¯Σ ab and ¯ n ab are symmetric and real ( i.e. , Hermitian) and commute in theultra-local limit (see Eq. 5.2), the two tensors share a common orthogonal system of eigen-vectors .Let { ψ , ψ , ψ } be an orthogonal eigensystem for ¯Σ ab and ¯ n ab and let σ i and ν i be theeigenvalue corresponding to the eigenvector ψ i , i.e. , ¯Σ ab · ψ i = σ i ψ i , ¯ n ab · ψ i = ν i ψ i , ( i = 1 , , , (B.1)where there is no summation over i . Then, from the evolution equations (5.3-5.4), it is im-mediately apparent that all time derivatives can be eliminated and replaced by combinations– 35 –f the original tensors ¯Σ ab and ¯ n ab . It follows that { ψ i } are also eigenvectors of the timederivatives ˙¯Σ ab and ˙¯ n ab : ˙¯Σ ab · ψ i = − (cid:16)(cid:0) N − (cid:1) ¯Σ ab + N (cid:0) n c (cid:104) a ¯ n b (cid:105) c − ¯ n cc ¯ n (cid:104) ab (cid:105) (cid:1)(cid:17) · ψ i ≡ λ i ψ i , (B.2) ˙¯ n ab · ψ i = − (cid:16)(cid:0) N − (cid:1) ¯ n ab − N ¯ n c ( a ¯Σ b ) c (cid:17) · ψ i ≡ η i ψ i . (B.3)Using the orthogonality of the eigensystem, (cid:0) ψ j · ψ i (cid:1) = 0 for i (cid:54) = j , and the expressionsabove, one can compute ψ j · d t (cid:16) ¯Σ ab · ψ i (cid:17) in two different ways that must necessarily beequivalent: ψ j · d t (cid:16) ¯Σ ab · ψ i (cid:17) = ˙ σ i (cid:8)(cid:8)(cid:8)(cid:8)(cid:8)(cid:42) (cid:16) ψ j · ψ i (cid:17) + σ i (cid:16) ψ j · ˙ ψ i (cid:17) (B.4) = λ i (cid:8)(cid:8)(cid:8)(cid:8)(cid:8)(cid:42) (cid:16) ψ j · ψ i (cid:17) + ψ j · ¯Σ ab · ˙ ψ i = σ j (cid:0) ψ j · ˙ ψ i (cid:1) , meaning that for i (cid:54) = j and σ i (cid:54) = σ j , ψ j and ˙ ψ i are orthogonal to one another, ψ j · ˙ ψ i ≡ , (B.5)such that ˙ ψ i is identical to ψ i up to stretching; if all eigenvectors are simply stretched overtime, then the subspace spanned by an individual eigenvector is unchanged or, equivalently,every eigenvector at time t remains an eigenvector. Note that each eigenvector ψ i keepsits own direction and the triad { ψ , ψ , ψ } does not undergo an overall rotation. The sameargument can be repeated for ¯ n ab and its eigenvalues ν i .It remains to discuss the degenerate case. We note that the shear tensor is defined tobe trace-free ( σ + σ + σ = 0 ) and must therefore have at least two different eigenvalues,unless all σ i ( i = 1 , , are zero. So the only non-trivial degenerate case is one in whichtwo eigenvalues are the same and one is different. Without loss of generality, let us assumethat σ = σ and σ (cid:54) = σ , and that ψ and ψ are two linearly independent eigenvectorswith eigenvalue σ at time t . Together, ψ and ψ span the two-dimensional subspace ofeigenvectors with eigenvalue σ . In this case, the same sort of argument as above shows thatthe time-evolution of ψ and ψ maps them into the same two-dimensional subspace, althoughin general they could be stretched and rotated. Furthermore, the time-evolved ψ and ψ andevery linear combination thereof have the same eigenvalue. This can be seen from Eqs. (B.2-B.3) by Taylor-expanding the shear and intrinsic curvature tensors around the initial time t and operating on an arbitrary linear combination of the t = t eigenvectors ψ and ψ : ¯Σ ab ( t + ∆ t ) · ( aψ + bψ ) (cid:39) ¯Σ ab ( t ) · ( aψ + bψ ) + ˙¯Σ ab ( t )∆ t · ( aψ + bψ ) (B.6) = (cid:0) σ + λ ∆ t (cid:1) · ( aψ + bψ ) , where we have used the fact that Eqs. (B.2) imply λ must equal λ . As in the non-degeneratecase, every eigenvector at time t remains an eigenvector. A similar analysis applies to ¯ n ab .Accordingly, the dynamics in the ultra-local limit is completely determined by the time-evolution of the eigenvalues. – 36 – Linearized evolution equations in the ultra-local limit
In Sec. 5, we identified the end states that we found numerically as the only dynamical attrac-tors, i.e. critical points of the evolution scheme in the ultra-local limit, given by Eqs. (5.5-5.10). There, we only listed the equations as linearized around the three critical points thatwe previously identified and listed in Table 1. Here, we present the perturbed system aslinearized for an arbitrary background solution, which underlied our calculations. (We shalldenote perturbed quantities by δ , e.g. , the i th linearized shear eigenvalue is given by δσ i . Allother variables are background quantities.)The linearly perturbed system of ordinary differential equations describing the dynamicsin the ultra-local limit around an arbitrary background solution takes the following form: δ ˙ σ = (cid:16) − N (cid:17) δσ − (cid:16) σ + (cid:0) ν − ν − ν (cid:1) ν − (cid:0) ν − ν (cid:1) (cid:17) δ N (C.1) − N (cid:16) ν − ν − ν (cid:17) δν + N (cid:16) ν + 2 (cid:0) ν − ν (cid:1)(cid:17) δν + N (cid:16) ν − (cid:0) ν − ν (cid:1)(cid:17) δν ,δ ˙ σ = (cid:16) − N (cid:17) δσ − (cid:16) σ + (cid:0) ν − ν − ν (cid:1) ν − (cid:0) ν − ν (cid:1) (cid:17) δ N (C.2) − N (cid:16) ν − ν − ν (cid:17) δν + N (cid:16) ν + 2 (cid:0) ν − ν (cid:1)(cid:17) δν + N (cid:16) ν − (cid:0) ν − ν (cid:1)(cid:17) δν ,δ ˙ ν = (cid:16) N (cid:0) σ − (cid:1)(cid:17) δν + (cid:0) σ − (cid:1) ν δ N + 2 N ν δσ , (C.3) δ ˙ ν = (cid:16) N (cid:0) σ − (cid:1)(cid:17) δν + (cid:0) σ − (cid:1) ν δ N + 2 N ν δσ , (C.4) δ ˙ ν = (cid:16) − N (cid:0) σ + 2 σ + 1 (cid:1)(cid:17) δν − (cid:16) σ + σ + (cid:17) ν δ N − N ν (cid:16) δσ + δσ (cid:17) , (C.5) δ ˙¯ W = − (cid:16) N − (cid:17) δ ¯ W − (cid:16) W − √ ε ¯ V (cid:17) δ N + √ ε N δ ¯ V , (C.6)where δ N ≡ −N δ ( N − ) and δ ¯ V are given by δ ( N − ) = (cid:16) ν − ν − ν (cid:17) δν + (cid:16) ν − ν − ν (cid:17) δν + (cid:16) ν − ν − ν (cid:17) δν (C.7) + (cid:16) σ + σ (cid:17) δσ + (cid:16) σ + σ (cid:17) δσ + ¯ W δ ¯ W ,δ ¯ V = − δ ( N − ) + 2 (cid:16) σ + σ (cid:17) δσ + 2 (cid:16) σ + σ (cid:17) δσ + 2 ¯ W δ ¯ W . (C.8)
D Conventional critical-point analysis using Friedmann variables
In Sec. 5, we identified the FRW scaling attractor solution as the stable end state for allphysically plausible initial conditions (assuming V ,φ /V = −√ ε = constant ). This solution iswidely known in cosmology as the exact solution of the Einstein-scalar field equations whenassuming an exponential potential. Here, we present the conventional derivation from thecosmology literature and relate the commonly used quantities from there with our variables.The starting point is a gravitational action involving a single scalar field φ with canonicalkinetic energy and and and exponential potential that is minimally coupled to gravity: S = (cid:90) d x √− g (cid:18) R − ∇ α φ ∇ α φ − V ( φ ) (cid:19) . (D.1)Evaluating the action for the flat homogeneous but anisotropic Kasner-like metric ds = − dτ + a ( τ ) (cid:88) i e β i ( τ ) (D.2a)– 37 – here β ( τ ) + β ( τ ) + β ( τ ) = 0 (D.2b)and τ is the proper FRW time coordinate, we find the following system of evolution andconstraint equations H − (cid:0) β (cid:48) + β (cid:48) + β (cid:48) (cid:1) = φ (cid:48) + V ( φ ) , (D.3a) H (cid:48) + (cid:0) β (cid:48) + β (cid:48) + β (cid:48) (cid:1) = − φ (cid:48) , (D.3b) β (cid:48)(cid:48) i + 3 Hβ (cid:48) i = 0 , ( i = 1 , , , (D.3c) φ (cid:48)(cid:48) + 3 Hφ (cid:48) = − V ,φ . (D.3d)As before, prime denotes differentiation with respect to proper FRW time τ . In the following,we shall eliminate β using the identity (D.2b).Next, we define the dimensionless variables that are often quoted as ‘Friedmann vari-ables,’ x = φ (cid:48) H , y = (cid:112) | V | H , u = β (cid:48) H , v = β (cid:48) H ; (D.4)and the dimensionless time variable N a = ln (cid:0) a/a (cid:1) , (D.5)measuring the number of e -folds of contraction starting from a = a . (By definition, N isnegative if the universe contracts and positive if it expands.) Note that, in the ultra-local limit,the Friedmann variables u, v are identical to the shear eigenvalues σ , σ and the remainingFriedmann variables x and y can be identified with ¯ W and ¯ V , respectively.Using the dimensionless Friedmann variables, the Hamiltonian constraint (D.3a) reducesto y = − x + 2 (cid:16) u + v + uv (cid:17) , (D.6)and the homogeneous system of evolution equations (D.3b-D.3d) takes the simple form x , N a = (cid:16) x + V ,φ V (cid:17) y , (D.7a) y , N a = (cid:16) V ,φ V x + y + 3 (cid:17) y, (D.7b) u , N a = u y , (D.7c) v , N a = v y . (D.7d)The system admits the following critical point solutions for V ,φ /V = −√ (cid:15) < : ( √ , , , ε = 3 , FRW) (D.8a) ( − V ,φ /V, (cid:113) ( V ,φ /V ) − , , ε > , FRW) (D.8b) ( (cid:113) − (cid:0) u + v + uv (cid:1) , , u, v ); (KL) . (D.8c)Obviously, these are the same critical points we found and listed in the first two rows ofTable 1. However, since no (homogeneous) spatial curvature is included, the third criticalpoint as listed in the third row of Table 1 cannot be recovered by means of this analysisIt is straightforward to linearize this simple autonomous system around each of thecritical points and conclude that the FRW scaling solution is a dynamical attractor if | V ,φ /V | – 38 –s sufficiently large, in agreement with our analysis in Sec 5. Yet, again, this stability analysisis limited by assuming a flat Kasner-like metric as given in Eq. (D.2a). For example, it is not possible to recover the role of the (homogeneous) spatial curvature in driving the systemaway from a Kasner-like and towards the FRW solution.Notably, for an exponential potential as given in Eq. (4.4), finding the stable criticalpoint, x = 2 ε, V ( τ ) = (3 − ε ) H ( τ ) , u = v = 0 , immediately yields the well-known FRWscaling attractor solution given in Eq. (5.19) after a series of simple integrations, using x / d ( H − ) (cid:48) . References [1] Lars Andersson and Vincent Moncrief. Elliptic hyperbolic systems and the Einstein equations.
Annales Henri Poincare , 4:1–34, 2003.[2] Richard L. Arnowitt, Stanley Deser, and Charles W. Misner. Dynamical Structure andDefinition of Energy in General Relativity.
Phys. Rev. , 116:1322–1330, 1959.[3] James M. Bardeen. Gauge Invariant Cosmological Perturbations.
Phys. Rev. , D22:1882–1905,1980.[4] James M. Bardeen, Olivier Sarbach, and Luisa T. Buchman. Tetrad formalism for numericalrelativity on conformally compactified constant mean curvature hypersurfaces.
Phys. Rev. D ,83:104045, 2011.[5] James M. Bardeen, Paul J. Steinhardt, and Michael S. Turner. Spontaneous Creation of AlmostScale - Free Density Perturbations in an Inflationary Universe.
Phys.Rev. , D28:679, 1983.[6] V.A. Belinsky, I.M. Khalatnikov, and E.M. Lifshitz. Oscillatory approach to a singular point inthe relativistic cosmology.
Adv. Phys. , 19:525–573, 1970.[7] L.T. Buchman and James M. Bardeen. A Hyperbolic tetrad formulation of the Einsteinequations for numerical relativity.
Phys. Rev. D , 67:084017, 2003. [Erratum: Phys.Rev.D 72,049903 (2005)].[8] William G. Cook, Iryna A. Glushchenko, Anna Ijjas, Frans Pretorius, and Paul J. Steinhardt.Supersmoothing through Slow Contraction, 6 2020.[9] Jürgen Ehlers. Beiträge zur relativistischen Mechanik kontinuierlicher Medien.
Mainz AkademieWissenschaften Mathematisch Naturwissenschaftliche Klasse , 11:792–837, January 1961.[10] G.F.R. Ellis. Dynamics of pressure free matter in general relativity.
J. Math. Phys. ,8:1171–1194, 1967.[11] F.B. Estabrook and H.D. Wahlquist. Dyadic analysis of space-time congruences.
J. Math.Phys. , 5:1629–1644, 1964.[12] Frank B. Estabrook, R.Steve Robinson, and Hugo D. Wahlquist. Hyperbolic equations forvacuum gravity using special orthonormal frames.
Class. Quant. Grav. , 14:1237–1247, 1997.[13] David Garfinkle, Woei Chet Lim, Frans Pretorius, and Paul J. Steinhardt. Evolution to asmooth universe in an ekpyrotic contracting phase with w > 1.
Phys. Rev. , D78:083537, 2008.[14] U. H. Gerlach and U. K. Sengupta. Gauge invariant perturbations on most general sphericallysymmetric space-times.
Phys. Rev. , D19:2268–2272, 1979.[15] Robert Geroch. Gauge, diffeomorphisms, initial-value formulation, etc. In Piotr T. Chruścieland Helmut Friedrich, editors,
The Einstein Equations and the Large Scale Behavior ofGravitational Fields , pages 441–477, Basel, 2004. Birkhäuser Basel.[16] Alan H. Guth. Eternal inflation and its implications.
J. Phys. , A40:6811–6826, 2007. – 39 –
17] Anna Ijjas and Paul J. Steinhardt. Bouncing Cosmology made simple.
Class. Quant. Grav. ,35(13):135004, 2018.[18] Anna Ijjas and Paul J. Steinhardt. A new kind of cyclic universe.
Phys.Lett. , B795:666–672,2019.[19] Pascual Jordan, Jürgen Ehlers, and Wolfgang Kundt. Republication of: Exact solutions of thefield equations of the general theory of relativity.
General Relativity and Gravitation ,41(9):2191–2280, September 2009.[20] Justin Khoury, Burt A. Ovrut, Paul J. Steinhardt, and Neil Turok. The ekpyrotic universe:Colliding branes and the origin of the hot big bang.
Phys. Rev. , D64:123522, 2001.[21] Wolfgang Kundt. “Golden Oldie”: The Spatially Homogeneous Cosmological Models.
GeneralRelativity and Gravitation , 35(3):491–498, March 2003.[22] E. Lifshitz. Republication of: On the gravitational stability of the expanding universe.
J.Phys.(USSR) , 10:116, 1946. [Gen. Rel. Grav.49,no.2,18(2017)].[23] Viatcheslav F. Mukhanov. Quantum Theory of Gauge Invariant Cosmological Perturbations.
Sov. Phys. JETP , 67:1297–1302, 1988. [Zh. Eksp. Teor. Fiz.94N7,1(1988)].[24] Frans Pretorius. Numerical relativity using a generalized harmonic decomposition.
Class.Quant. Grav. , 22:425–452, 2005.[25] Frans Pretorius. Simulation of binary black hole spacetimes with a harmonic evolution scheme.
Class. Quant. Grav. , 23:S529–S552, 2006.[26] Paul J. Steinhardt and Neil Turok. Cosmic evolution in a cyclic universe.
Phys.Rev. ,D65:126003, 2002.[27] Paul Joseph Steinhardt. Natural inflation. In G.W. Gibbons, Hawking. S., and S. Siklos,editors,
The Very Early Universe , pages 251–66. Cambridge University Press, 1983.[28] Henk van Elst and Claes Uggla. General relativistic (1+3) orthonormal frame approachrevisited.
Class. Quant. Grav. , 14:2673–2695, 1997.[29] Henk van Elst, Claes Uggla, and John Wainwright. Dynamical systems approach to G(2)cosmology.
Class. Quant. Grav. , 19:51–82, 2002.[30] Alexander Vilenkin. The Birth of Inflationary Universes.
Phys.Rev. , D27:2848, 1983.[31] Jr. York, James W. Gravitational degrees of freedom and the initial-value problem.
Phys. Rev.Lett. , 26:1656–1658, 1971., 26:1656–1658, 1971.