Dynamical attractors in contracting spacetimes dominated by kinetically coupled scalar fields

We present non-perturbative numerical relativity simulations of slowly contracting spacetimes in which the scalar field driving slow contraction is coupled to a second scalar field through an exponential non-linear sigma model-type kinetic interaction. These models are important because they can generate a nearly scale-invariant spectrum of super-Hubble density fluctuations fully consistent with cosmic microwave background observations. We show that the non-linear evolution rapidly approaches a homogeneous, isotropic and flat Friedmann- Robertson-Walker (FRW) geometry for a wide range of inhomogeneous and anisotropic initial conditions. Ultimately, we find, the kinetic coupling causes the evolution to deflect away from flat FRW and towards a novel Kasner-like stationary point, but in general this occurs on time scales that are too long to be observationally relevant.


Introduction
Cosmic microwave background observations have shown that the spectrum of temperature anisotropies is nearly scale invariant and gaussian with an average amplitude of one part in hundred thousand and no detectable B-mode polarization thus far. According to our current understanding, the temperature anisotropies are an imprint of primordial curvature fluctuations sourced by quantum excitations of one or more scalar fields, the energy density of which dominated the early-universe [2,14].
In order for the curvature fluctuations to match the observed spectrum, by the time the modes are generated, the cosmological background must be smooth and flat as described by a Friedmann-Robertson-Walker (FRW) geometry with line element ds 2 = −dτ 2 + a(τ ) 2 dx i dx i , where a(τ ) is the FRW scale factor. Then, once the modes are generated, the smoothing mechanism must continue for at least 60 additional e-folds. Neither of these two conditions is trivial to satisfy. For example, if the FRW solution is a robust attractor for a wide range of initial conditions, most of the volume is eventually smoothed and flattened by the time the smoothing phase ends; but, if the smoothing is not sufficiently rapid, most volume converges to a smooth and flat FRW spacetime with only a few e-foldings remaining before the end of smoothing. In this case, modes generated around the 60 e-fold mark would carry an imprint of an unsmoothed inhomogeneous and/or anisotropic geometry and would thus not match the observed spectrum.
Slow contraction [3], a primordial phase that connects to the hot expanding phase through a gentle classical bounce, has been shown to be both robust and rapid. The phase can be reached via a canonical scalar field φ that is minimally coupled to Einstein gravity and has a negative potential V (φ). The scalar field energy density naturally evolves to dominate the total stress-energy while driving the geometry to a smooth and flat FRW space-time. The robustness of slow contraction as a dynamical attractor -its insensitivity to initial conditions including those that lie outside the perturbative regime of FRW spacetimes -was recently shown in Refs. [6,9,11]; the remarkable rapidity, with smoothing and flattening typically occurring within less than 10 e-folds of contraction of the Hubble radius, was demonstrated in Ref. [9].
The goal of this paper is to examine whether the powerful smoothing property remains in cases in which the smoothing scalar field φ is coupled to a second scalar field χ through an exponential non-linear σ model-type kinetic interaction. These models are important because it has been shown that they can lead to the generation of a nearly scale-invariant spectrum of super-Hubble density fluctuations fully consistent with cosmic microwave background observations [12,13]. More generally, this study is important for exploring whether slow contraction remains a powerful, robust and rapid dynamical attractor even when there is an exponentially strong kinetic interaction with a secondary scalar field and whether any distinctive features result compared to the case of scalar fields with canonical kinetic energy density.
As demonstrated in Ref. [11], the key to rapid and robust smoothing is ultralocality: Starting with arbitrary initial data that lies outside the perturbative regime of the FRW state, contracting spacetimes rapidly evolve to an anisotropic and spatially curved state where gradients, i.e., spatial derivatives, are suppressed relative to the other so-called velocity contributions. These causally separated regions then each independently evolve to the homogeneous, isotropic and spatially flat FRW state which is the only stable stationary point of the ultralocal limit.
Once the flat FRW state is reached, physical distances between objects shrink proportional to the scale factor a(τ ) and exponentially slower than the Hubble radius, where φ 3; hence the name: slow contraction. The rate at which |H −1 | contracts is determined by the effective equation of state associated with the scalar field on an FRW background: where p φ = 1 2 φ 2 −V (φ) is the co-moving 'pressure' and ρ φ = 1 2 φ 2 +V (φ) the co-moving energy density of the scalar field φ in the homogeneous FRW limit and prime denotes differentiation with respect to τ .
For a negative exponential potential, as will be considered throughout this paper, the scaling attractor solution of the Einsteinscalar system in the flat FRW limit is given by and τ is running from large negative to small negative values. (Here and throughout the paper, reduced Planck units, M Pl −2 ≡ 8πG N with G N being Newton's constant, are used.) For a potential with the characteristic mass scale M ∼ 0.1 M Pl , the effective equation of state φ = 50 is such that the scale factor contracts only by a factor of 2 while the Hubble radius shrinks by a factor of 2 50 . A particularly important feature of slow contraction is the fact that, because the Hubble radius shrinks much faster than the scale factor, the wavelengths of fluctuations (which are proportional to a(τ )) necessarily end up on super-Hubble scales by the end of slow contraction. Unlike inflation, though, slow contraction is a 'supersmoother,' meaning that adiabatic (a.k.a. curvature) modes decay, whether they are of classical or quantum origin [4]. The decay is the opposite of what occurs in expanding universes and is associated with the fact that the evolution of the adiabatic modes during slow contraction is subject to an antifriction-like term (H 0) that has the opposite sign than the friction-like term (H 0) in the expanding case. This is an important and appealing feature of slow contraction because it suppresses quantum runaway effects. At the same time, a mechanism is needed to generate the spectrum of temperature anisotropies observed in the cosmic microwave background and the seeds for structure formation.
In Refs. [12,13], it has been shown that entropy modes (i.e. pressure fluctuations on hypersurfaces of constant energy density) generated by quantum fluctuations of a second scalar field during slow contraction can fulfill this role. First, like adiabatic modes, the entropic modes are generated by quantum fluctuations and their wavelengths also end up on super-Hubble scales. Second, unlike the adiabatic modes, the entropy modes can experience a net red shift effect. This can occur, for example, if the modes are sourced by a light scalar field χ that is kinetically coupled to the background field φ through a non-linear σ-type interaction, e.g., where g is the four-metric determinant and R the Ricci scalar.
Assuming an exponential coupling with a characteristic mass scale m M as will be considered throughout this paper, the stable attractor solution of the Friedmann-scalar system of equations with U (χ) 0, is that χ is constant (χ ≡ 0) while a(τ ) and φ(τ ) evolve with time according to the scaling solution given by Eq. (1.4). The constant χ solution is stable [12] because, on an FRW background, the kinetic coupling 1 2 κ(φ)∇ µ χ∇ µ χ enters the evolution equation (2.25) of the χ-field in a way that adds to the Hubble anti-friction a friction-like term, Since the combination is positive, the χ field's kinetic energy experiences a net damping as if it were in a de Sitter-like background and freezes out at some constant value. Notably, quantum fluctuations of the χ field also experience a net de Sitter-like damping term 3H + φ /m 0 such that their amplitude grows, leading to a nearly scale-invariant and gaussian spectrum of entropy modes on super-Hubble wavelengths. Finally, it has been shown that the entropy modes can source curvature modes, e.g. when exiting slow contraction and entering the bounce stage, as shown in Refs. [7,8].
The remaining question is whether adding a χ field with the exponential non-linear σ model-type kinetic interaction with φ needed to generate a nearly scale-invariant and gaussian spectrum of curvature perturbations preserves the robustness and rapidity of slow contraction as established for the single field scenario. In this paper, we address this issue by adapting the mathematical and numerical techniques developed for the single-field case in Refs. [3,5,6,9,11].
Our non-perturbative analysis yields some surprising results that could not be obtained using the conventional methods of cosmological perturbation theory. First, in the special case where χ is precisely massless, we find that the evolution can be rapidly deflected away from the flat FRW stationary point if the characteristic scale M associated with V (φ) is too close to the Planck scale (M = 1), or equivalently, if ε φ is not sufficiently large; instead, the evolution is driven towards a Kasner-like solution in which the gradient of χ,S χ x (τ, x), is non-zero and time-independent ∂ τSχ x ≡ 0. (1.11) Second, in generic models which have moderately smaller values of M 0.1 or weakly broken shift symmetry in χ, e.g., a small mass for χ, we find that the deflection effect is strongly suppressed such that the evolution leads to a long-lived state with negligible S χ x (τ, x) and the FRW scaling solution in Eq. (1.4), similar to the impressive single-field result. This state persists long enough and has just the conditions needed to generate a nearly scale-invariant spectrum of curvature perturbations that can account for temperature fluctuations observed in the cosmic microwave background. The kinetic coupling ultimately causes the evolution to deflect from flat FRW, but on a time scale that is too long to be relevant for cosmologies that connect to the hot expanding phase through a smooth (non-singular) bounce.

Numerical Scheme
For the non-perturbative, numerical analysis, we shall adapt the orthonormal tetrad form of the Einstein-scalar field equations corresponding to the action given in Eq. (1.5), as developed for the single-field case. As per convention, g µν is the spacetime metric, and G µν is the Einstein tensor. (Here and for the remainder of the paper, we express dimensional quantities in reduced Planck units where M Pl = 1.) Throughout, spacetime indices (0−3) are Greek and spatial indices (1 − 3) are Latin. The beginning of the alphabet (α, β, γ or a, b, c) denotes tetrad indices and the middle of the alphabet (µ, ν, ρ or i, j, k) denotes coordinate indices. In the following, we only give a brief overview with the goal of making the paper selfcontained. A comprehensive description of the formulation as well as a complete derivation of the evolution and constraint equations are provided in Refs. [6,11].

Variables
As with any tetrad formulation of the field equations, spacetime points are represented through a family of unit basis four-vectors, or vierbein, {e 0 , e 1 , e 2 , e 3 }, where e 0 is the timelike four-vector and the spacelike four-vectors of the triad {e 1 , e 2 , e 3 } each lie in the rest three-space of e 0 . The local Lorentz frame set by the tetrad basis is flat, i.e., with · denoting the inner product of the tetrad and η αβ = diag(−1, 1, 1, 1) being the Minkowski metric.
The geometric variables of the formulation are the sixteen tetrad vector components {e α µ } and the twenty-four Ricci rotation coefficients where ∇ λ is the projection of the spacetime covariant derivative ∇ µ onto the tetrad e λ , defines the nine components of the shear tensor; and is the induced curvature tensor associated with the spatial triad where abc denotes the Levi-Civita-symbol. The eighteen components of K ab and N ab are dynamical variables. The three-vectors b a ≡ γ a00 , Ω a ≡ 1 2 a bc γ cb0 (2.8) are frame gauge quantities with b a defining the proper local acceleration of the congruence and Ω a defining the local angular velocity of the spatial triad relative to Fermi-propagated axes.
The geometric variables must be supplemented by the variables describing the two scalars. These are the field distributions φ, χ, their velocities and gradients, respectively, as detailed below in Section 2.3.

Gauge fixing
For our numerical scheme, we fix the six gauge degrees of freedom of the tetrad frame in a way that makes the connection of the geometric variables to physical quantities straightforward: -we fix the spatial triad {e 1 , e 2 , e 3 } to be inertially non-rotating a.k.a. Fermi propagated (Ω a ≡ 0); -we require the shear tensor to be symmetric (K ab ≡ K (ab) ), meaning that the tetrad congruence is hypersurface-orthogonal such that it defines a particular foliation of spacetime into spacelike hypersurfaces of constant time {Σ t } with e 0 being the future directed timelike unit normal to {Σ t } and the spatial triad vectors being tangent to {Σ t }.
With this choice of tetrad frame gauge, K ab denotes the extrinsic curvature of {Σ t } and the components of N ab are the spatial (or intrinsic) curvature variables. Note that the acceleration vector b a is implicitly fixed by this gauge choice through where x 0 is the time coordinate of {Σ t } and × denotes scalar multiplication. Furthermore, we must write the tetrad evolution and constraint equations in the form of partial differential equations (PDEs) by way of which we shall numerically evolve the geometric and scalar field variables specified on an initial spacelike hypersurface. To do so is particularly straightforward given our choice of a hypersurface-orthogonal tetrad because, in this frame gauge, elements of the transformation matrix {λ α µ } between tetrad and coordinate basis four-vectors, are easily identified with quantities of the 3+1 Arnowitt-Deser-Misner (ADM) formalism, namely: 11) where N is the ADM lapse, N i the ADM shift and the coordinate metric is being given by g µν = η αβ λ α µ λ β ν . Finally, the directional derivatives along the tetrads can be written as To fix the four coordinate gauge degrees of freedom, we specify the lapse function N and the shift vector N i by requiring that The CMC slicing condition fixes the lapse function in that it leads to an elliptic equation (2.17) for N ; and -the spatial coordinates are co-moving (N i = 0), i.e. constant along both the congruence and the foliation.
As emphasized previously in Refs. [6,11], the particular choice of our coordinate gauge has several advantages: since the trace of the extrinsic curvature Θ −1 is spatially uniform on each {Σ t }, we can define the time coordinate t to track Θ, i.e., such that, in the homogeneous limit, Θ is the Hubble radius |H −1 |. In addition, we can use Θ to rewrite our equations in terms of dimensionless Hubble-normalized variables , where N is the Hubble-normalized lapse and bar denotes normalization by the mean curvature Θ −1 on constant time hypersurfaces. The combination of the time coordinate t tracking the three-curvature with Hubblenormalized variables is particularly useful for our purposes to study the rapidity and robustness of slow contraction because it enables us to run the simulation for any finite period without encountering singular behavior or stiffness issues: first, the putative singularity is at t → −∞ and, second, the rapidly changing Hubble radius is not entering the numerical calculation as a dynamical variable, leaving the slowly changing scale factor as the only relevant dynamical factor.

Evolution and constraint equations
Taking everything together, the Einstein-scalar field evolution equations (2.1-2.3) in Hubblenormalized, orthonormal tetrad form yield an elliptic-hyperbolic system: The lapse N is defined at each time step through an elliptic equation, while the evolution of the remaining geometric as well as the scalar field variables is given by a hyperbolic system of PDEs: where curved brackets denote symmetrization X (ab) ≡ 1 2 (X ab +X ba ) and angle brackets denote traceless symmetrization defined as X ab ≡ X (ab) − 1 3 X c c δ ab . The geometric variables are the symmetric and antisymmetric components, respectively of the Hubble-normalized, spatial curvature tensorN ab ;Σ ab is the trace-free extrinsic curvature tensor, The variablesW φ ,W χ denote the Hubble-normalized scalar field velocities. In addition, the evolution equations are supplemented by a set of constraints which we shall use to specify the initial data as well as to verify convergence of the numerical computation: The variablesS φa ,S χa denote the Hubble-normalized scalar field gradients. Note that, for simplicity, we rescaled the velocity and gradient terms of the χ field with the non-linear σ-type kinetic interaction κ(φ), as can be seen, e.g., in Eqs. (2.25) and (2.35).

Initial data
To study the robustness of the kinetically-coupled two-field model as given through Eq. (1.5) to cosmic initial conditions, it is essential to perform a large number of numerical relativity simulations corresponding to a wide range of initial conditions including those that lie far outside the perturbative regime of homogeneous and isotropic FRW spacetimes.
As detailed in Refs. [6,11], our numerical scheme allows for the variation of all freely specifiable geometric and scalar field variables, {Ē a i ,n ab ,Ā b ,Σ ab } and {φ, χ,W φ ,W χ }, respectively. To ensure that the initial data satisfy the constraint equations (2.30-2.35), in particular, energy and momentum conservation, we adapt, as in our earlier work, the York method [15] that is commonly used in numerical relativity studies. (If the constraint equations are satisfied at the initial time, the Einstein equations propagate them such that they are satisfied at all later times, a condition that is checked numerically.) Employing the same tetrad frame and coordinate gauge conditions that we detailed above in Sec. 2.2, we first fix the value of the inverse mean curvature Θ 0 of the spatial hypersurface Σ t 0 at some initial time t 0 and then define the three-metric of Σ t 0 to be conformally-flat, i.e., where the conformal factor ψ is not a free function but is determined by an elliptic equation (given in Eq. 3.11 below) upon setting all other variables. Note that the conformally-flat metric choice in Eq. (3.1) does not impose a limitation on cases that can be studied since it does not propagate; in fact, it is immediately violated after the first evolution step.
Our choice of Θ 0 and g ij (t 0 , x) fixes the coordinate components of the spatial tetrad basis vectors,Ē and the intrinsic curvature variables, , it is straightforward to verify that the constraints (2.32 and 2.33) are trivially satisfied. Furthermore, with a conformally-flat spatial metric, the momentum constraint (2.31) reduces to the following simple relation, denote the conformally rescaled Hubble-normalized shear and scalar field velocity variables, respectively. For simplicity and without loss of generality, we choose the initial value of φ and χ to be zero, turning the momentum constraint (3.4) into the condition that the initial shear component Z 0 ab be divergence free and giving us full freedom to choose the initial scalar field velocities Q φ and Q χ . We fix these quantities as follows: and where the parameters ξ, a 1 , a 2 , b 1 , b 2 , f φ , µ φ , d φ , f χ , µ χ , d χ , and Q 0 are constants. The sinusoidal form reflects the fact that, for the numerical simulation, we choose periodic boundary conditions 0 ≤ x ≤ 2π with 0 and 2π identified. For simplicity, all deviations from homogeneity are along a single spatial direction x, as in Ref. [6]. Finally, with the initial value of all freely specifiable geometric and scalar field variables Θ 0 ,Ē a i ,Ā b ,n ab , Z 0 ab , φ, χ, Q φ , Q χ fixed and satisfying the constraint equations (2.31-2.35), we impose the remaining Hamiltonian constraint (2.30) on these variables. This yields an elliptic equation for the conformal factor ψ, which we solve numerically.

Numerical results
To numerically solve the Einstein-scalar field equations, we discretize the elliptic-hyperbolic system (2.17-2.27) using second order accurate spatial derivatives and a three-step method for time integration employing the Iterated Crank-Nicolson algorithm. At each sub-step, we first solve the elliptic equation (2.17) for the Hubble-normalized lapse N through a relaxation method and then update the hyperbolic equations (2.18-2.27) to the next Iterated Crank-Nicolson sub-step. In the simulations presented below, we use a grid of 4096 points with ∆ x = 2π/4096 and a Courant factor of 0.5. To demonstrate the convergence of our code, the error and convergence was analyzed for a broad range of examples using the same methods as detailed in the Appendices of Refs. [6,9]. To summarize those tests, our code shows no signs of numerical instability and exhibits clear second order convergence at early times. At later times when a smooth, ultralocal spacetime develops, we empirically see the convergence improve to third order.
In this section, we present three representative examples from our extensive numerical studies of slow contraction with kinetically coupled scalar fields φ and χ given through the action in Eq. (1.5) with kinetic coupling and potential energy densities respectively. We begin with the special case that m χ = 0 and demonstrate a subtle instability compared to the single-field case in Refs. [6,9,11]. We then show that, for the generic case with m χ = 0, the instability is suppressed and robust and rapid convergence to a long-lived flat FRW state occurs.
For each example, we show the evolution of the total scalar field energy density (Ω φ−χ ), shear (Ω s ) and curvature (Ω k ) contributions to the normalized energy density, defined as: where Ω φ−χ + Ω s + Ω k = 1. We also formally define the separate contributions of the two fields, ) Recall that the gradientS χ a includes dependence on φ through the coupling κ(φ); see Eq. (2.35).
The specify the parameters of the initial shear components Z 0 ab and scalar field velocities Q φ , Q χ given in Eqs. (3.8) and (3.9-3.10), respectively, such that they correspond to initial conditions with highly non-perturbative deviations from a flat FRW spacetime: In addition, for the three representative cases, we fix the model parameters where the ratio M/m is chosen such that the predicted tilt of the temperature fluctuation spectrum matches current observations [1]. Note that for ξ, a 1 , a 2 , b 1 , b 2 , f φ , µ φ , d φ , Q 0 , and V 0 , we chose the same values as in the single field case studied Ref. [6], to facilitate comparison. The time coordinate runs from initial time t = 0 (or Θ 0 = 3) towards −∞ (or Θ → 0). Equivalently, the time can be characterized by n H ≡ −t, the number of e-folds of contraction of the inverse mean curvature Θ, where n H runs from zero towards +∞. In practice, models with a classical non-singular bounce undergo slow contraction until n H ≈ 120 before the bounce occurs [10]. For each example studied, the same initial data was evolved with several resolutions to confirm second order convergence; the highest resolution has 4096 points on the base level.
The three representative cases correspond to different choices for M and m χ : The choice m χ = 0 is a special case where the action (1.5) has a shift symmetry (χ → χ+ const.). The value of M = 0.2 ( φ = 13) was shown in the single field studies to be just above the minimum required for robust and rapid smoothing and flattening [6,9]. In this case with two fields and the special shift symmetry, though, we find a different outcome. Fig. 1 shows four snapshots of the evolution of the normalized scalar field energy density Ω φ−χ (green), spatial curvature Ω k (red) and shear Ω s (blue). The first snapshot (n H = 0) shows the initial non-perturbative deviations from flat FRW. By the second snapshot, a short time later (n H = 9), the spacetime approaches flat FRW with Ω φ−χ ≈ 1 and Ω s ≈ Ω k ≈ 0, seemingly similar to the single-field case. But instead of remaining smooth and flat as found in the single-field case, something different occurs: the flat FRW phase is not stable, and the shear (blue curve) begins to grow, as illustrated in the last two snapshots. At n H = 25, the geometry is spatially flat but with a mix of field energy plus shear. This combination was termed as 'Kasner-like' in the single field case [6] where it was found in certain examples with values M bigger than 0.2 M Pl . As we will see, though, something quite different is happening in the case of two kinetically coupled fields. Note that spikes will form at locations where the χ-gradient vanishes; i.e., due to the ultralocal nature of the evolution, such points get stuck at the FRW state, whereas adjacent regions, destabilized by a non-zero χ-gradient, transition to the final Kasner-like states. Thus, using the terminology of Ref. [6], the state should be referred to as 'Kasner-like (modulo spikes).' Fig. 2, which tracks the evolution of Ω φ and Ω χ separately, reveals more details. As the sequence progresses, it can be seen that Ω φ rapidly comes to dominate, but then energy is transferred to χ through their non-linear coupling such that Ω χ begins to grow. More precisely,W χ → 0 almost immediately, by n H = 2, and remains negligible; however, the rescaled χ-gradientS χ x = κ(φ)D x χ grows rapidly proportional to the coupling factor, κ(φ), that is growing rapidly due to the growth of φ which is quickly rolling downhill its potential V (φ). (N.B. This growth of the gradient is real, i.e. not a frame effect, as explained in the Appendix A.) By the last snapshot n H = 25, which is far short of the bounce, Ω χ grows to dominate over Ω φ ; and, then, looking back to the last snapshot in Figure 1, we see that, at the same time, the gradient sources a growing shear component Ω s . We stop the code at n H = 25 for reasons described in the next section where we also use analytics to determine how the evolution continues. Here again we consider the special case with m χ = 0 in which the action (1.5) has a shift symmetry. This example illustrates that slightly decreasing M (or, equivalently, increasing the equation of state during slow contraction φ ) significantly delays the onset of the instability but does not eliminate it, as illustrated by this example. Fig. 3 shows that the universe is rapidly flattened as in Case I, but in this case the shear Ω s only begins to dominate for some ranges of x at n H = 120. If there were no bounce and the simulation were running for a yet longer range of n H , the behavior would be similar to the Case I but shifted in time. As a practical matter, though, by decreasing M modestly further, M < 1/15, the instability could be delayed beyond n H = 120, which is more than sufficient for bouncing cosmologies with a smooth (non-singular) bounce (where the bounce occurs at n H < 120). Now we turn to the generic case of m χ = 0, thus breaking the shift symmetry, which we show also acts to suppress the instability. (This case is generic because there is no reason to expect an exact shift symmetry in χ since there is no shift symmetry in φ.) Bearing in mind that the initial value of the inverse extrinsic curvature is Θ −1 0 10 −42 GeV in bouncing cosmologies, Figs. 4 and 5 show that it suffices to break the shift symmetry in χ with even a small mass m χ in order to obtain a qualitatively different result. The reason will be explained in Section 5 below. Fig. 4 shows that, beginning from the same initial conditions as in Cases I and II, the geometry robustly and rapidly (by n H = 4) converges to flat FRW, just as in the single field case. (A similar result is found for M = 0.2 for a slightly greater value of m φ .) In particular, the gradientsS x χ and the shear Ω s remain small over the n H = 120 e-folds of the simulation, sufficient for practical applications to bouncing cosmologies with a gentle (non-singular) bounce. Fig. 5 shows a subtle difference from the single-field case, though. Although the total normalized scalar field energy density Ω φ−χ is close to unity and the shear is negligible, the kinetic coupling between the two fields through the gradientS x χ result in an exchange of energy Ω x x x x Figure 4: For the case with M = 0.1 and m χ = 300Θ −1 0 , snapshots of the normalized energy density following the same color code as in Figs. 1 and 3. Compared to the case in Fig. 3, one observed that introducing a small m χ suffices to obtain robust and rapid smoothing. between the two components, as shown by the compensating oscillations at n H = 120, albeit an exchange that generates negligible shear. These can be viewed as classically generated entropic fluctuations; in this example and for a wide range of M and m χ , these fluctuations have an amplitude that is irrelevant compared to the quantum-generated entropic fluctuations on the length scales measured by cosmic microwave background observations. Their presence is a sign that the evolution is beginning to deviate away from flat FRW; in principle, given additional time, the system would evolve to a Kasner-like fixed point as in the cases above (but this has no practical relevance for cosmologies involving a slow contraction phase that connect to the hot expanding phase through a smooth non-singular bounce.)

Analytic approximation
The numerical studies in the previous section show that slow contraction in the two-field models considered in this paper is in general a rapid and robust smoother over the time scales and length scales of interest for bouncing cosmologies. However, the three numerical results also show that there is a subtlety that does not arise in the case of a single canonical scalar field with a steep negative potential. Namely, a flat FRW state is not the ultimate fixed point attractor; although the evolution initially approaches a flat FRW fixed point, it is subsequently deflected towards a homogeneous but anisotropic 'Kasner-like' state.
To understand this novel phenomenon, we study different stages in the evolution using analytic perturbative analyses. These explain how the deflection from flat FRW arises and show that that the characteristic time scale for the deflection is in general of O(100) or more e-foldings of contraction of the inverse mean curvature, which is of no practical relevance in cosmologies where slow contraction connects to the hot expansion phase through a smooth classical (non-singular) bounce.
The evolution beginning from highly non-perturbative deviations from flat FRW involves four stages: Stage 1: During the first stage, the system evolves from the freely specified initial state (as detailed above in Sec. 3) to one that is ultralocal, meaning that terms involving spatial derivatives (gradients) are quickly dominated by velocity terms (i.e., terms that do not involve spatial derivatives) as the evolution proceeds. In particular, within only a few e-folds of contraction of the inverse mean curvature Θ, Note that ultralocal does not mean flat FRW, as extensively detailed in Ref. [11].
Stage 2: For M ≤ 0.2 and Q 0 > 0, the same range considered in Ref. [6], the Einstein-scalar system (2.18-2.27) rapidly and non-linearly evolves from an ultralocal (but not FRW) state towards the flat FRW stationary point with just as we observe by n H = 9 in all of our simulations as illustrated in the second panels of Figures 1, 3 and 4.
The first two stages are very rapid (complete by n H 10 ) and are determined by the slow contraction sourced by the φ-field. As result, the evolution during these two stages is similar to what is found for the single-field case.
Stage 3: For the case of a single canonical scalar, flat FRW is a stable fixed point of the evolution. For the kinetically-coupled two-field models considered here, though, flat FRW is not a stable fixed point. Rather, there remain small deviations from flat FRW that eventually grow large enough to deflect the evolution from the flat FRW fixed point. This effect can be seen by perturbing the Einstein-scalar system of equations (2.18-2.27) around the flat FRW fixed point given by Eqs. (5.2-5.3). Inspecting the linearized system, where δ denotes linear perturbations around the background solutions, it is immediately apparent that the perturbations of all geometric variables as well as the linearized scalar field variables δW φ , δS φa form a closed system and decay at the same rate as in the single field case as |t| grows, i.e., 14) where N = 2M 2 ≤ 2/5. 3), it is straightforward to verify that the kinetic interaction κ(φ) = e −φ/m makes theW χ = 0 solution stable. In matrix form, the closed system can be written as The two eigenvalues λ ± corresponding to the coefficient matrix are given by Although the analysis to this point might suggest that the flat FRW fixed point is stable, there remains the (rescaled) spatial gradient of χ to consider, which turns out to be the source of the instability. More precisely, the same factor κ(φ) that stabilizes the time derivative of χ destabilizes theS χ x 0 gradient contribution: Evaluating the evolution equation (5.12) of S χ x for the flat FRW stationary point solution in Eq. (5.3), with M/m = 1.015, it becomes apparent that the gradient perturbation slowly grows, where δS χ x (t FRW ) is the value of δS χ x at the beginning of Stage 3, when the flat FRW stage begins. This expression is valid at the linear level during Stage 3, the small entropic fluctuations shown in the last panel of Fig. 5 arise as non-linear gradient effects on the evolution grow to become non-negligible.

Discussion
Using the tools of numerical general relativity, we studied the cosmological evolution during slow contraction in models where the stress-energy is sourced by two kinetically-coupled scalar fields with an exponential non-linear σ model-type interaction beginning from inhomogeneous and anisotropic initial conditions that deviate far from flat FRW spacetimes. Our main finding was that slow contraction led to a robust and rapid convergence to a flat FRW geometry, similar to previous results obtained in Refs. [3,6,11] for the case where slow contraction is sourced by a single canonical scalar field. However, our study revealed a subtle difference as well. Whereas the flat FRW solution is a stable attractor in the single-field case, it is not in the case of the two kinetically coupled fields considered here. Instead, the evolution in the two-field case rapidly becomes ultralocal and approaches close to the flat FRW fixed point where it remains for a considerable period, but ultimately it is deflected away from that fixed point and towards a Kasner-like fixed point. We investigated the instability analytically and showed that the characteristic time for the instability to develop enough to cause the deflection is O(100) or more e-folds of contraction of the inverse mean curvature Θ.
In addition, we showed that slightly decreasing the characteristic scale M of the negative exponential potentialV = −V 0 e −φ/M associated with the field φ that drives slow contraction as well as increasing the massm χ of the light χ field which couples to φ through a non-linear σ type kinetic interaction further increases the duration of the flat FRW period and delays the deflection to the Kasner-like fixed point.
The result is reminiscent of the single-field case discussed in Ref. [9] where, for some initial data sets, we found that the system first evolved close to a Kasner-like fixed point before it was deflected to the stable flat FRW attractor solution. There are two key differences between the two cases, though: -in the single field case, the system approaches close to the (unstable) Kasner-like fixed point for only a very small subset of initial conditions and it remains there for only a short period before being deflected to the stable attractor fixed point (the flat FRW state), -in the the kinetically-coupled two-field case, the system approaches close to the (unstable) flat FRW fixed point for a very wide range of initial conditions and remains there for a long period before being deflected to the stable attractor Kasner-like fixed point.
This result is surprising and was not anticipated in previous studies based on the conventional perturbative techniques commonly applied in cosmology. Hence, this study is a fine demonstration of the power of non-perturbative, numerical relativity to reveal novel cosmological dynamics that one could not anticipate through conventional perturbative methods but that are critically important to understand in developing cosmological models.
In the present example, the difference between the single-field and kinetically-coupled two-field case turns out not to be observationally relevant in bouncing cosmologies because the end of slow contraction and the transition to the hot expanding phase begins around 100 e-foldings of contraction of Θ -before the flat FRW state destabilizes. This is significant because the kinetically-coupled two-field models studied here are examples that can generate a nearly scale-invariant spectrum of super-Hubble density fluctuations fully consistent with cosmic microwave background observations.

A Tetrad frame transformation rules
In this Appendix, we consider the relationship between geometric variables as represented in the hypersurface-orthogonal tetrad frame (used in our tetrad codes to study smoothing of initial conditions) and the the co-moving tetrad frame. To this end, we will first derive the transformation rules for the tetrad vector components, the Ricci rotation coefficients and the (effective) 'fluid' variables describing the stress-energy under the Lorentz boost {Λ α β } that connects two arbitrary tetrad frames. Then, we apply the rules to the case of scalar fields to transform geometric and 'fluid' variables from the hypersurface-orthogonal tetrad frame to the co-moving tetrad frame.
We are interested in tetrad frame transformations under Lorentz boosts {Λ α β } that transform a time-like tetrad e 0 to anotherẽ 0 . The boost is defined through the rapidity β between the two 4-vectors e 0 andẽ 0 and through the projection {w a } ofẽ 0 into the 3-surfaces spanned by the spatial triad e a (a = 1, 2, 3): where w a w a = 1. Note that Γ ≡ coshβ = 1/ √ 1 − v 2 is the Lorentz factor with v a ≡ tanhβ w a . In particular, Vierbein. The tetrad frame vectors {e α } and {ẽ α } (with α = 0, ..., 3) are related as follows: Ricci Rotation Coefficients. The Ricci Rotation Coefficients γ αβλ which are defined as with ∇ λ ≡ e λ µ ∇ µ transform as follows: Components of the Stress-Energy Tensor. Characterizing an arbitrary stress-energy tensor through effective 'fluid' variables, j a ≡ −e 0 µ e a ν T µν , (A.8) where is the energy density, j a the three-momentum flux, s ab the spatial stress tensor and p denotes the pressure and π ab denotes the trace-free anisotropic stress (with π ab ≡ s (ab) − 1 3 s a a δ ab ), under a Lorentz boost as defined in Eq. (A.1), the 'fluid' variables transform as follows:T i.e.,˜ For the numerical relativity codes that we use to study the robustness and rapidity of slow contraction, we fixed the tetrad frame gauge to be hypersurface-orthogonal and Fermi propagated, as described in Sec. 2.2. In particular, the timelike vierbein e 0 is normal to spacelike hypersurfaces. In general, j a = 0 in this frame, which means that, typically, e 0 and the effective fluid's 4-velocity (which is the timelike congruenceẽ 0 of the co-moving tetrad) do not coincide.
For example, in the case of a single scalar field, where the 'fluid' variables , j a , s ab , and p in the hypersurface-orthogonal tetrad frame take the following form: Here, D 0 denotes the Lie derivative along e 0 and D a is the directional derivative along e a .
Manifestly, e 0 = −D a φ/D 0 φ ad hence the hypersurface-orthogonal and co-moving tetrad frames do not coincide in general. The Lorentz boost, v a = −D a φ/D 0 φ, transforms the tetrad to the co-moving frame, where the stress-energy takes the form of a perfect fluid: j a = π ab = 0. (A.21) Note, though, that the hypersurface-orthogonal tetrad frame converges to the co-moving frame if v a = (D a φ/φ) → 0, as is the case in our simulations that involve a single scalar which is minimally-coupled to gravity and has a negative potential. By contrast, our initial data as specified above in Sec. 3 is not represented in the co-moving frame. This does not affect the conclusion that the spacetime converges to a homogeneous, isotropic and spatially flat FRW universe. However, it does mean that the initial conditions are not as a Eulerian co-moving observer would measure them to be. More generally, the Lorentz boost v a = j a /( + p c ) with p c being the co-moving pressure transforms an arbitrary tetrad frame to the co-moving frame. If, in addition, s ab = j a j b +p c δ ab , it is straightforward to show that the 'fluid' takes the perfect fluid form in the co-moving tetrad frame with˜ = 1 Γ 2 + p c − p c ,s ab = p c δ ab , j a = π ab = 0. (A.22) In the case of two kinetically-interacting scalars like those considered in this paper with T µν = ∇ µ φ∇ ν φ + κ(φ)∇ µ χ∇ ν χ − 1 2 ∇ µ φ∇ µ φ + 1 2 κ(φ)∇ µ χ∇ µ χ + V (φ) + U (χ) g µν , (A.23) the 'fluid' variables , j a , s ab , and p in the hypersurface-orthogonal tetrad frame take the following form: and the co-moving 'fluid' variables take the following form: In general, interacting scalar fields act as 'imperfect fluids' with non-zero momentum flux and non-diagonal spatial stress tensor.