Stable cosmological Kaluza-Klein Spacetimes

We consider the Einstein flow on a product manifold with one factor being a compact quotient of 3-dimensional hyperbolic space without boundary and the other factor being a flat torus of fixed arbitrary dimension. We consider initial data symmetric with respect to the toroidal directions. We obtain effective Einsteinian field equations coupled to a wave map type and a Maxwell type equation by the Kaluza-Klein reduction. The Milne universe solves those field equations when the additional parts arising from the toroidal dimensions are chosen constant. We prove future stability of the Milne universe within this class of spacetimes, which establishes stability of a large class of cosmological Kaluza-Klein vacua. A crucial part of the proof is the implementation of a new gauge for Maxwell-type equations in the cosmological context, which we refer to as slice-adapted gauge.

where M (4) corresponds to the macroscopic spacetime and B is a compact q-dimensional Riemannian manifold referred to as internal space. The latter models compactified dimensions practically invisible to observers. Identifying the ground state of Kaluza-Klein theory has been a long standing open problem, which may be considered in different contexts. Original works show semiclassical instabilities in the case M (4) is equipped with the Minkowski metric [Wi]. If results on the Einstein-Maxwell system [BZ, Sp], which relate to the special case B = S 1 and M (4) being equipped with the Minkowski metric, are excluded, then mathematically rigorous nonlinear stability or instability of Kaluza-Klein spacetimes in the context of classical general relativity was unknown until recently. In a recent work Wyatt established stability of Kaluza-Klein spacetimes for the class of models, where M (4) carries the Minkowski metric and B is a symmetric q-dimensional torus [Wy], (1.2) g K K = η M (4) + g flat,T q .
With this imposed symmetry the vacuum Einstein equations on M (4+q) reduce to an Einstein-wave map-Maxwell type system on M (4) , which is shown to have the Minkowski metric as its stable ground state. From the perspective of classical general relativity this result justifies the interpretation of the corresponding higher-dimensional Kaluza-Klein background spacetime (M (4+q) , g K K ) as the ground state of the generalized higher dimensional field equations.

Cosmological spacetimes.
A prerequisite for the stability analysis of the class of Kaluza-Klein spacetimes with Minkowski space as their macroscopic part is the corresponding nonlinear stability result for the classical 4-dimensional vacuum Einstein equations [CK, LR]. In the class of asymptotically flat spacetimes Minkowski spacetime is the only solution known to be stable. The analogous problem for the Kerr family is still open. There is only one other spacetime in the class of solutions of the Einstein equations with vanishing cosmological constant that is known to be stable, which is the Milne model. This solution belongs to the class of cosmological spacetimes, i.e. it has spatial slices with compact topology that carry a negative Einstein metric γ. The Milne model is future complete and past incomplete and its future nonlinear stability problem has been resolved in the vacuum setting by . This result covers also the higher-dimensional case, however, not in the sense of compactified dimensions. In analogy to the asymptotically flat case we ask for the natural ground state for Kaluza-Klein theory in the class of cosmological spacetimes. It will be shown in this paper that the generalized Kaluza-Klein spacetime arising from the Milne model, reading is future nonlinearly stable for perturbations that obey the symmetry of g flat,T q .
1.3. Main theorem. The stability theorem formulated in terms of terminology introduced in Section 2 is given in the following.
1.3.1. Result. We formulate the main result in terms of rescaled variables, adapted to the evolution. To obtain the final reduced equations two rescalings are performed. The first is a conformal rescaling necessary to avoid regularity problems arising from the Kaluza-Klein reduction. The second rescaling uses the CMC-time function to obtain variables which are scale free and independent of the expansion. The detailed reduced system is given in (3.13)-(3.18). We consider initial data sets consisting of a Riemannian metric g , the trace-free part of the second fundamental form Σ, a set of one-forms A and its time deriv-ativeȦ as well as set of wave-type maps Φ and their time derivativesΦ fulfilling the reduced constraint equations (3.13) -(3.14).
(1.5) (g , Σ, A, Φ) → (γ, 0, 0, Φ ∞ ) for some set Φ ∞ of constant functions. In particular, the Milne model is an attractor for the macroscopic geometry of product spacetimes with a torus as an internal space within the class of perturbations that preserve the full symmetry group of the torus.
REMARK 1.2. To obtain a suitable solution theory for (3.9), we need to impose a gauge condition, e.g. the Lorentz gauge, which turns (3.9) into a hyperbolic equation. In the main theorem, we have imposed the Lorentz gauge and the initial data (A,Ȧ) is meant as initial data with respect to this hyperbolic equation. However, to control the long-time behaviour of (3.9), a different gauge turned out to be more suitable, see Subsection 1.3.2 below. REMARK 1.3. The Kaluza-Klein reduced Einstein equations restrict all possible perturbations of the background to those where the torus remains symmetric. This, however, still allows for each point in the macroscopic space that the torus (the internal space at this point) evolves within the class of flat tori. Moreover, this evolution within the moduli space of the internal space is not uniform in the macroscopic space but can vary in general in a sufficiently regular way. While all fields do not depend on the point in internal space, the dynamics of the geometry of the internal space is captured within this approach.
REMARK 1.4. By the conformal rescaling we perform in Section 3.2 of this paper, the physical Riemannian metricsḡ of the spatial hypersurfaces satisfy after a natural rescaling that τ 2ḡ → det(Φ ∞ ) − 1 2 γ as τ → 0, where τ is the mean curvature function. As Φ ∞ is constant, this limit is also negative Einstein but with a possible different constant. It is interesting to note that the macroscopic geometry encoded inḡ gets a contribution from the geometry of the internal space via the above rescaling. REMARK 1.5. Our main result is also interesting for the stability analysis of classical vacua in string theory since toroidal compactifications are used as toy models in string theory [Po05,Chapter 8].
We comment in the following on some technical aspects of the stability proof.

Gauging Kaluza-Klein fields:
The slice-adapted gauge. A standard gauge for a vector potential that is used to consider Maxwell-type equations is the Lorentz gauge ∇ µ A µ = 0. One obtains a nonlinear wave equation of second order on A. However, it turns out to be surprisingly difficult to analyze this equation in the present context and to construct a natural energy which yields optimal bounds for the decay of the perturbation. A source of this difficulty may arise from the fact that the vector potential A is not determined by the Lorentz gauge as this gauge is preserved by transformations A → A + d f if f = 0 and thus has infinitely many degrees of freedom. To overcome this problem, we choose a gauge which is adapted to a foliation of the spacetime by spacelike hypersurfaces and which uniquely determines A: We demand that the spatial components of A, ω, associated to this foliation are orthogonal to the kernel of the Hodge Laplacian and that the time component of A, Ψ, regarded as a function of the spacetime has vanishing integral on each spatial slice div g ω = 0, ω ⊥ ker(∆ H ),ˆM ΨdV g = 0. (1.6) In this gauge, the Maxwell equation is a wave equation on the spatial part of A coupled to an elliptic equation for its time component. Details are provided in Lemma 6.4 and Proposition 6.5. To the best of our knowledge, such a gauge has not been used in the context of related problems so far. However, the slice-adapted gauge seems to be useful for Maxwell-type equations on other spacetimes with compact spatial hypersurfaces (e.g. on the de-Sitter space).
1.3.3. Regularity aspects and the momentum constraint. Another interesting aspect of the Kaluza-Klein reduced system is the fact that the momentum constraint, which is not explicitly used in controlling the perturbation in the pure 3+1-dimensional vacuum stability proof, does play an important role in the present problem in the following sense. Below, we will use energies that control the H 4 -norm of an evolving metric g (in terms of a fixed background metric) and the H 3 -norm of the tracefree part Σ of the second fundamental form. However, when differentiating the energies for the perturbation of the fields generated by the internal space, one obtains s derivatives of Σ and 3 derivatives of its time derivative. Those in turn can not be controlled by the H 3 -norm of Σ and the H 2 -norm of its time derivative, respectively. A closer analysis however reveals that these terms only appear as third derivatives of div g Σ and second derivatives of ∂ T div g Σ. Replacing those terms using the momentum constraint improves the regularity by one order and closes the estimate.
1.4. Related systems. Theorem 1 has some immediate consequences for related systems and in particular automatically implies the following results.
1.4.1. Einstein-Maxwell-Dilaton system. In the special case of B = T q = S 1 the 5-dimensional U (1)symmetric vacuum field equations with S 1 being the symmetry direction reduce to the 4-dimensional Einstein-Maxwell-Dilaton system [OW]. This implies in particular the following corollary. COROLLARY 1.6. The Milne model is future stable as a trivial solution to the Einstein-Maxwell-Dilaton system.
We use the terminology trivial solution in the sense that it is actually a solution to the Einstein vacuum equations where the matter fields vanish.
Setting the field Φ to zero and removing the corresponding equation from the system, we obtain a new system that is formally equivalent to the classical Einstein-Maxwell system. This implies COROLLARY 1.7. The Milne model is future stable as a trivial solution to the Einstein-Maxwell system.

1.4.2.
Brans-Dicke theory. Another well-known system that is captured by our main result is the Brans-Dicke model of general relativity. This system is obtained by setting the one-forms A to zero. The Brans-Dicke model couples pure gravitation with a scalar field in which the value of the scalar field can be interpreted as a dynamical version of Newton's gravitational constant, see [OW] for more details.
COROLLARY 1.8. The Milne model is future stable as a trivial solution to the Brans-Dicke model.

U(1)-symmetric spacetimes.
There is a third relation of Theorem 1.1 with previously considered models, where in this specific case the present result can be considered as a higher-dimensional analog. In their work on the stability of certain Bianchi type-III models Choquet-Bruhat and Moncrief consider spatial topologies of the form Σ×S 1 , where Σ is a closed two-dimensional higher genus surface [CM]. The background solution being investigated is −4d t 2 +2t 2 σ Σ +d x 2 , where σ Σ is a metric of constant negative scalar curvature on Σ. They prove future stability of this solution considered within the set of solutions to the 4-dimensional vacuum Einstein equations obeying a U (1)-symmetry in the S 1 direction. By a Kaluza-Klein reduction this symmetric system is equivalent to the 2+1-dimensional Einstein equations on R × Σ with a source term given by a massless scalar field. In a way this can be seen as an analogue to the problem considered in the present work, where the Kaluza-Klein fields are replaced by a single massless scalar field. However, the approach of Choquet-Bruhat and Moncrief does not carry over to higher dimensions as it relies on the particular features of the 2+1-dimensional geometric setting. Those are for instance the existence of a monotone L 2 -energy and the usability of the momentum constraint to control the trace-free part of the second fundamental form. In 3+1-dimensions these methods are not available and need to be replaced by the energies provided by . Nevertheless, the structure of a torus bundle over a negatively curved compact Riemannian manifold is present in both cases. The result in this paper implies that those geometries are stable under the Einstein flow irrespective of the low dimensional features used in [CM].
1.4.4. Higher-dimensional backgrounds. Finally, we mention that by the methods used in this paper, one can also prove nonlinear stability of a higher-dimensional Kaluza-Klein Milne model under the same class of perturbations. Here, γ is a metric of constant sectional curvature −m 2 on a compact m-dimensional manifold. In higher dimension, the conformal behaviour of the Maxwell-type equation yields a faster decay of F µν and improves the energy estimates.
1.5. Organization of the paper. This paper is organized as follows. In Section 2 notations are introduced as well as the rescaling of the macroscopic geometry and several auxiliary quantities. In Section 3 we perform the Kaluza-Klein reduction and derive the reduced Einstein-wave map-Maxwell system. In Section 4 we compute the energy-momentum tensor in the reduced Einstein equations in terms of the fields generated by the internal space and introduce norms to estimate them. Sections 5 presents the elliptic estimates for the macroscopic lapse function and the shift vector field. Section 6 derives energy estimates for all evolution equations individually and thereby constitutes the core step of the stability analysis. Section 7 presents the proof of the main theorem and Section 8 presents all related systems listed above for which our stability analysis of the Milne model applies.
Acknowledgements. D.F. has been supported by the Austrian Science Fund (FWF) project P29900-N27 Geometric Transport equations and the non-vacuum Einstein flow. V.B. gratefully acknowledges the support of the Austrian Science Fund (FWF) through the project P30749-N35 Geometric variational problems from string theory.

PRELIMINARIES
2.1. Notation. Throughout this paper, M is a compact manifold eventually equipped with different Riemannian metrics and I ⊂ R is an open interval. In this paper, the appearing Lorentzian metrics on M = I × M will be denoted by h, and the associated covariant derivative will be denoted by ∇. The wave operator associated to h is defined with the sign convention such that = tr h ∇ 2 . In this paper, we will sometimes also denote Lorentzian metrics byh,h,ĥ and the associated covariant derivatives and wave operators will be denoted by ∇, ∇,∇ and , andˆ , respectively. Riemannian metrics on M will be denoted by g ,g and the associated covariant derivatives will be denoted by D, D, respectively. The Laplacian of g is defined as ∆ = tr g D 2 and the volume form will be denoted by dV g . The exterior derivative acting on differential forms on M is denoted by d and the formal adjoint with respect to g is d * . The Hodge-Laplacian acting on differential forms is then ∆ H = d * d + d d * . The Lie-Derivative of a tensor T in the direction of a vector field X will be denoted by L X T . Throughout this paper, Greek indices α, β, γ, . . . will denote spacetime coordinates on I × M and Latin indices i , j , k, . . . will denote coordinates on M . The coordinates on the torus T q will be denoted by m, n, p, . . .. The index 0 will either refer to a time coordinate or to a timelike vector field. Its meaning will be clarified in the subsection where it is used.

The macroscopic spatial background geometry.
In what follows we consider M equipped with a negative Riemannian Einstein metric γ with Ric[γ] = − 2 9 γ fixed once and for all. The Einstein operator ∆ E associated with γ acting on symmetric 2-tensors, ∆ E ≡ −∆ − 2R, has trivial kernel, i.e. ker ∆ E = {0}. This fact is relevant for the features of the natural energy associated with ∆ E . This has been discussed in [AF17] and is mentioned here for the sake of completeness.
2.3. Geometric formalism for the evolving spacetime. In the following sections, we will study the evolution of a 3+1-dimensional Lorentzian metrich (more precisely of its rescaled version h introduced below). For this purpose, we will now introduce some geometric quantities that will be used throughout the paper. In the ADM formalism,h is written as and the tracefree part of the second fundamental form of the hypersurfaces {τ = const } is denoted by Σ.
Here we assume that these hypersurfaces all have constant mean curvature and that the mean curvature of {τ = const } is τ. We define rescaled quantities g , N , Σ, X by It is easily seen that with respect to this new time coordinate, the above Lorentzian metric is given bỹ Let Π be the second fundamental form of the slice {T ≡ const } with respect to the Lorentzian metric h. Then one can show that The future-directed timelike unit normal of the hypersurfaces {T ≡ const } with respect to h is We use e 0 to split 1-forms on M described in the following. For A ∈ Ω 1 ( M), we define a function Ψ ∈ C ∞ ( M) and a time-dependent family of one-forms ω ∈ C ∞ (I , Ω 1 (M )) by Ψ := A(e 0 ) and ω(∂ i ) = A(∂ i ). Throughout the paper, we will view any ω ∈ C ∞ (I , Ω 1 (M )) as an element in Ω 1 ( M ) by demanding ω(e 0 ) = 0. This allows us to write the above splitting as A = ω + Ψe * 0 where e * 0 ∈ Ω 1 ( M) is the dual of e 0 . We compute the connection coefficients for the rescaled Lorentzian metric h. Using the Koszul formula, one shows k are coordinates on M and the index 0 refers to the vector field e 0 given in (2.6). The following lemma is technically relevant for computations performed further below.
Proof. At first, we compute Let {∂ 1 , ∂ 2 , ∂ 3 } be local coordinate fields on M such that D ∂ i ∂ j = 0 at some fixed point p and with respect to a fixed metric g t 0 . We then extend these local vector fields to elements in C ∞ (I , X(M)) by defining Then at the point (t , p), we compute (2.10) Therefore, we obtain (2.11) The third formula follows from the second and the fact that [L e 0 , D] f = 0.

KALUZA-KLEIN REDUCTION
In this section we perform the Kaluza-Klein reduction beginning with the physical Lorentzian metric on the full spacetime.
3.1. Kaluza-Klein metrics. Following [CH09,p. 653] we consider the Kaluza-Klein ansatz for the Lorentzian metricĥ on R × M × T q , where q ≥ 1 denotes the dimension of internal space, Here,h αβ denotes a Lorentzian metric on R× M , {Φ mn } is a set of functions on M and A i α is an R q -valued 1-form on M . Moreover, θ α and θ m are suitable co-frames on M and T q , respectively. We obtain for the macroscopic part of the Ricci tensor (cf. [CH09,p. 659 Here, Φ mn 1≤m,n≤q is the inverse of the matrix {Φ mn } 1≤m,n≤q and F is the curvature of A given by In particular, T [Φ, F ] determines the matter source terms in the effective 3+1-dimensional Einstein equations. The remaining parts of the Einstein vacuum equation yield the equations of motion for the fields F and Φ, cf. [CH09,p. 659,eq. (5.3,5.4)]. Those equations read 3.2. Conformal rescaling. From an analytical point of view the second order terms of Φ on the righthand side of (3.3) are problematic. In this section, we therefore perform a standard conformal rescaling of the metric that yields an equivalent system that has a better analytic structure. Let us recall some standard transformation formulas. Ifh = e 2uh , we havē where n is the dimension of the spacetime. As a conformal factor, we set with a constant c whose value is to be determined. We have . We now put c = −1/4, n = 4 and the relation between the physical metrich and the conformal metrich ish = e 2uh = 1 detΦh . Then forh, F and Φ, the equations read ( 3.11) 3.3. Macroscopic Einstein equations. (3.11) are the effective macroscopic Einstein equations. Their energymomentum tensor, arising from the geometry of the internal space, takes the form We use in the following the standard 3+1-dimensional ADM formalism for the spacetime metric as in (2.1), whereÑ ,g andX are lapse, physical metric and shift vector field. The matter quantities appearing in the ADM-Einstein equations in [Re] read Here, 0 refers to the time-function τ.
3.4. Rescaled system. We perform now the rescaling of the macroscopic Einstein equations according to (2.2). All symbols in the following denote the rescaled variables as in (2.2). Then, the Einstein flow in CMCSH gauge reads The relation between the rescaled matter quantities and the original ones is This set of equations is the basis for analyzing the dynamical behaviour of the perturbations of the geometry of the macroscopic space.
Before going on we define our notion of smallness. In what follows we say a solution or data is small or fulfills a smallness condition if for a sufficiently small ε > 0, where Φ b is a fixed constant map. By construction this condition holds for the initial perturbation we consider, then by local stability of the system, the condition holds on a finite time-interval. This justifies to make the smallness assumption to derive the decay estimates in the sense of a standard bootstrap argument.

ESTIMATING THE ENERGY-MOMENTUM TENSOR
In this section we evaluate the matter terms in the Einstein equations in terms of the wave-type-map and the one-forms determining the energy-momentum tensor. The index 0 corresponds to the τ timefunction in this section. We first compute the rescaled energy density . We evaluate T 00 and use the notation tr to compute the trace of the Lie-algebra indices The spatial part is given by We evaluate now the trace-part of η We evaluate the current Here, we require the off-diagonal components of the energy-momentum tensor. Those are given below

L 2 -norms of the Energy-Momentum tensor.
When analyzing the dynamics of the metric variables we require bounds on standard Sobolev norms of the matter variables as listed in (3.18). Those correspond directly to bounds on the Sobolev norms of components of the energy-momentum tensor. We derive those bounds in the following, expressed in terms of the corresponding norms of the matter fields F , A and Φ, respectively. We define some useful norms for this purpose. (4.7) Moreover, F is antisymmetric and the factor |τ| 2 compensates a growth of the 0-components of F relative to the pure spatial components. It needs to be determined/fixed as soon as the decay properties of F are understood. Note here that all objects and derivatives are defined with respect to the rescaled metric g , so that there is no more intrinsic scaling in this energy. . From the comparison of the bases we obtain that the coefficients are related via A τ = dT dτ A T = −τ −1 A T . Analogously, we define for the field Φ a similar Sobolev norm. (4.8) We obtain the following estimates for the components of the energy-momentum tensor as appearing in the matter variables. Recall, X = X /N . Then we find the following LEMMA 4.2. Let ℓ ≥ 3/2. Then the following estimate holds (4.9) and we denote byT the spatial part of the tensor T .
Proof. The estimates follow immediately from the expansions of the energy-momentum tensor components (4.2)-(4.4) and (4.6).
REMARK 4.3. Note that the constant, given that all arguments are uniformly bounded, as we assure by suitable bootstrap assumptions, can be considered a generic constant.
Proof. We evaluate first the term,T 00 . For the first estimate we evaluate (4.10) tr F µ 0 F µ0 = τ 2 (g i j − X i X j )δ mn F i 0,m F 0 j ,n and then apply Sobolov embedding with ℓ > 3/2. We evaluate next the square of F , which is (4.11) The related term containing a square of F , where δ mn is replaced by Φ mn can be decomposed identically by replacing δ by Φ. We consider exemplary the following term Similar decompositions hold for the other terms in the brackets on the right-hand side of the equation forT 00 . The remaining terms in the first line of that equation can be estimated directly. To deduce the full estimate it is sufficient to use the fact that the regularity is high enough to use product estimates for the Sobolev norm and that every time derivative of Φ and every zero-component of F appears with one additional τ factor. We turn now to the estimate forT 0i . We note that the term in the big brackets is identical to the case considered above. Since the last factor is now only a τ −3 and a shift term we obtain the first summand in the estimate. It remains to evaluate the first line of the evaluation ofT 0i . The terms containing derivatives of Φ can immediately be estimated. We evaluate the first term containing F .
These terms yield terms that decay like τ as contained in the estimate. We turn to the estimate forT now noting that the trace ofT can be treated similarly. The term in the large brackets is unchanged and is here multiplied only with a τ −2 factor. This leaves an overall τ 2 factor, which appears in the estimate.

Estimating the matter variables as appearing in the Einstein equations.
The final estimates for the rescaled matter quantities as appearing in the Einstein equations are the following: PROPOSITION 4.4. Let ℓ ≥ 3/2, then we have (4.14) Proof. This is an immediate consequence of the foregoing lemma as well as (4.1) and (4.5) and the definition of η and S.

ELLIPTIC ESTIMATES
We provide in this section the standard elliptic estimates for lapse and shift and their time-derivatives.
PROPOSITION 5.1. Under smallness conditions for the lapse function, a pointwise estimate of the form 0 < N ≤ 3 holds and moreover the following two estimates.
Proof. These estimates are an immediate consequence of elliptic regularity applied to (3.15) and (3.16), respectively and the maximum principle applied to (3.15).

Estimates of the time derivatives.
LEMMA 5.2. Let ℓ ≥ 4. For sufficiently small perturbations, the following estimate holds. (5.2) The constant C depends implicitly on the perturbation via Proof. By differentiation with respect to T the elliptic system implies where 〈D X , Σ, Σ〉 = D i X j Σ i k Σ k j and (5.5) We proceed analogous to [AF17] using the evolution equations for the energy-density and the current, which are independent of the matter model. The divergence identity of the energy momentum tensor in the unrescaled form, ∇ α T αβ (cf. [Re], (2.66), (2.67)) reads with respect to the rescaled variables, ρ =ρ|τ| −3 and  = |τ| −5 , The time derivative of the term containing η, however, requires a detailed evaluation as it depends on the equations of motion for the matter model. We need to estimate the H ℓ−2 -norm of ∂ T (τ 3 η). This term is τ 3 η = 4πτ −2g i jT i j and up to a constant and the factor τ −2 it is evaluated in (4.4). We now take the time derivative of the terms on the right-hand side of (4.4) modulo the τ 2 -factor and replace, if necessary, second time derivatives of the matter fields using the corresponding equations of motion. We do this explicitly for two terms to illustrate the computation and leave the remaining terms to the reader. This computation will provide an estimate for ∂ T τ 3 η H ℓ−2 . The first term we consider explicitly is The norms of the first terms on the right-hand side can directly be estimated. We focus on the evaluation of the last term.
Here we suppress terms that can either directly be estimated or those that can be handled similarly to the one considered explicitly. We proceed with that term.
There are two terms with time derivatives on the right-hand side, which cannot be estimated by the energies. We therefore replace those by the corresponding evolution equations or by suitable quantities estimated in the respective sections on the control of the matter fields. The relation between A, Ψ and ω and the definition of e 0 imply (5.10) We intend to use those equations to replace the left-hand side appearing in the time-differentiated equations by the right-hand sides, which can then be estimated using the corresponding results from section 6.2. We can estimate the Sobolev norms of Ψ and ∂ T Ψ by Lemma 6.8 and 6.9, respectively. The spatial components, i.e. ω and their time-derivatives are estimated using equivalency of energies as stated in (6.21). The term L e 0 L e 0 ω ℓ,n is substituted using equation (6.13). This, in turn, makes terms in Φ appear, for which we use the standard Sobolev norm. Proceeding as described leads to an estimate for the H ℓ−2 -norm of the left-hand side by . The second term from (4.4) that we estimate explicitly is Note that the last two terms are defined in (6.44) and estimated in (6.54). Evaluating the remaining terms of (4.4) after taking the time derivative we conclude an estimate of the following form. (5.13) Here, the suppressed terms are not in factors of the norms of ∂ T N and ∂ T X and therefore contribute directly to the right-hand side of the final elliptic estimate for the norm of ∂ T N . The terms, which are listed explicitly are handled in the following way. We note that the factor of the term ∂ T N H ℓ−2 is small by assumption. Applying elliptic regularity to (5.4) this term can therefore be absorbed in the constant. On the right-hand side however a term with ∂ T X H ℓ−2 remains, which also has a small factor in product. This preliminary estimate is then used in conjunction with the elliptic estimate for (5.5), which contains ∂ T N terms that are replaced by the preliminary estimate. This estimate in turn contains ∂ T X terms on the right-hand side, which can be absorbed using smallness of the factors and we obtain an estimate for ∂ T X H ℓ independent of ∂ T N . This can then in turn be used in the preliminary estimate for ∂ T N H ℓ to obtain the final estimate for ∂ T N H ℓ . 6. ENERGY ESTIMATES 6.1. Energy estimate for the geometry. We define the energy to measure the tracefree part of the second fundamental form and of the difference between the metric and the background metric as in the related work . We recall briefly some necessary notation. The lowest eigenvalue of the Einstein operator corresponding to the specific Einstein metric is denoted by λ 0 . For a relevant lower bound in the present case cf. [Kr15]. The correction constants α = α(λ 0 , δ α ) and c E are given by , c E = 1 λ 0 > 1/9 9(λ 0 − ε ′ ) λ 0 = 1/9 with δ α = 1 − 9(λ 0 − ε ′ ), where 1 >> ε ′ > 0 is a free variable to be chosen below. The energy is defined in the following. For m ≥ 1 let (6.2) The corrected energy is LEMMA 6.1. There exists a δ > 0 and a constant C > 0 such that for δ-small data (g , Σ, A, Φ) the inequality Proof. This is analogous to the previous work [AM-2] taking into account the triviality of the kernel of the Einstein operator.
The relevant energy estimate for the corrected energy is LEMMA 6.2. For sufficiently small E s we have Proof. The proof is analogous to the one in [AF17].
Substituting the norms of the matter quantities by Proposition 4.4 we obtain the energy estimate. PROPOSITION 6.3. Let s ≥ 5/2 and E s be sufficiently small. Then we have 6.2. Energy estimates for the vector potential. With respect to the Lorentzian metric h defined in (2.4), (3.9) can be written as Here and in the rest of the subsection, we omit the index l in equation (3.9) due to convenience. LEMMA 6.4 (Slice-adapted gauge). Let F ∈ Ω 2 ( M ) be exact. Then there exists a unique form A ∈ Ω 1 ( M) with d A = F such that div g ω = 0, ω ⊥ ker(∆ H ),ˆM ΨdV g = 0, (6.8) where ω and Ψ are defined in Section 2.3. In the statement and the proof of the lemma, d is the exterior derivative on M .
Proof. Let B ∈ Ω 1 ( M) such that d B = F . Let f ∈ C ∞ ( M ) with´M f dV g = 0 for each T ∈ I , c ∈ C ∞ (I ) and η ∈ C ∞ (I , Ω 1 (M )) be such that η ∈ ker(∆ H ) for each T ∈ I . Let Demanding the first condition of the lemma yields g i j D i B j +∆ g f = 0 and be-cause´M g i j D i B j dV g = 0, this equation can be uniquely solved at each time. Let ω 1 , . . . , ω L ∈ C ∞ (I , Ω 1 (M )) be for each T an L 2 (g )-orthonormal basis of ker(∆ H ) (Note that the dimension of ker(∆ H ) equals the first Betti number of M . Thus, it does not depend on g ). The second gauge condition is obtained by defining 〈B, ω a 〉 g dV g · ω a . (6.10) The third condition yieldŝ M A(e 0 )dV g =ˆM (B (e 0 ) + d f (e 0 )dV g + ∂ T c ·ˆM N −1 dV g , (6.11) which fixes ∂ T c. Uniqueness of A follows by construction. PROPOSITION 6.5. Let F ∈ Ω 2 ( M ) be exact and assume it solves (6.7). Let A ∈ Ω 1 ( M ) be a potential for F which satisfies the gauge conditions of Lemma 6.4 and let Ψ and ω be as in Lemma 6.4. Then we have the equations (6.13) Proof. Using (2.7), one computes where i e 0 F ∈ C ∞ (I , Ω 1 (M )) is given by where we have used that div g ω = 0. The first formula follows from (6.7). To prove the second formula, we compute, using (2.7) again, Therefore, the second formula again follows from (6.7). REMARK 6.6. Local existence for the system (6.12),(6.13) is argued as follows: One first solves (6.5) by using the Lorentz gauge h λµ ∇ λ A µ = 0. In this gauge, (6.5) becomes for which local existence follows from standard theory. Here, H,h denotes the Hodge wave operator of the metric h. By the construction in the proof of Lemma 6.4, we obtain in a unique way a pair (ω, Ψ) which solves the system (6.12),(6.13). On the other hand, as long as the solution (ω, Ψ) of (6.12),(6.13) is bounded, any corresponding solution A of (6.19) is also bounded: Let B = ω + Ψ · e * 0 and f be a solution of the equation (6.20) Then A = B + d f satisfies the Lorentz gauge and is bounded by construction. The main advantage of the slice-adapted gauge is that it is easier to control the solution (6.12),(6.13) by energy estimates.
The structure of the system (6.12),(6.13) motivates the following energy with k ≥ 1. Note that due to the gauge condition ω ⊥ ker(∆ H ) and elliptic regularity, the L 2 -norm of ω is controlled by E k (ω). LEMMA 6.7. Suppose that F solves (6.7), A ∈ Ω 1 ( M) is a gauged vector potential for F and Ψ and ω are as in Lemma 6.4. Then for k > n/2 + 1 and provided that N is uniformly positive and N H k is bounded by some fixed constant, we have the energy estimate Proof. First recall that the Hodge Laplacian is defined as ∆ H = d d * + d * d . By extending this definition to the exterior algebra Ω * (M ), we may also write ∆ H = (d + d * ) 2 where d + d * is a self-adjoint first-order differential operator acting on the exterior algebra. Fix l ∈ {0, . . . , k − 1}. By integration by parts, (6.23) At first, we compute (6.24) In the following we will make use of the * -notation to denote various contractions between tensors. Therefore, after integration by parts we get Now we estimate the commutator terms. Similarly as in Lemma 2.1, we have [L e 0 , (d + d * )]η = [L e 0 , d * ]η = Π * Dη + D log N * L e 0 η + S * η + Π * D log N * η (6.28) for a general differential form η ∈ C ∞ (I , Ω m (M )). Here, we used the notation S = 2div g Π − Dtr g Π. By induction, we get (6.29) and again by induction, (6.31) Therefore, by the bounds on N , (6.32) Similarly, the second commutator term is estimated aŝ (6.33) Finally, the second last term in (6.25) can be treated by (6.13) and standard estimates. LEMMA 6.8. As long as D log N H k + DΦ H k−1 is small enough, the function Ψ satisfies the estimate Proof. By elliptic regularity and the first equation in Lemma 6.5, (6.36) and by (6.15), Combining these estimates finishes the proof of the lemma. LEMMA 6.9. As long as D log N H k + DΦ H k−1 is small enough and Π H k−2 + Φ −1 H k−2 is uniformly bounded , the function Ψ satisfies the estimate Proof. By differentiating the first equation in Lemma 6.5 in the direction of e 0 and using elliptic regularity, we obtain (6.39) By using Lemma 2.1 and standard estimates, we get (6.40) Using the smallness of D log N H k + DΦ H k−1 , we can absorb the terms containing norms of ∂ e 0 Ψ into the left hand side of the equation. Consequently, by assuming in addition that Π H k−2 + Φ −1 H k−2 is uniformly bounded, we get Using the smallness assumptions again and treating L e 0 (L e 0 ω) H k−2 by the second equation in Lemma 6.5 and standard estimates, we arrive at the estimate of the lemma. PROPOSITION 6.10. We have the energy estimate as long as the norms of the appearing objects are uniformly bounded and N is uniformly positive.
6.3. Energy estimates for the functions. With respect to the metric h, given by (2.4), equation (3.10) is and with respect to the future-directed timelike unit normal e 0 , we can express this equation as where ( * * ) are the linear error terms and ( * ) are the nonlinear error terms. The global existence of a similar system, namely wave maps from a large class of expanding spacetimes, has been studied in [BK].
To write down the right energy, we consider the model ODË Proof. We consider an arbitrary summand of the energy. For convenience, we write u = Φ mn . For the rest of the proof, let l be even, the odd case is similar. By integration by partŝ It remains to consider the commutator terms. At first, we conclude from Lemma 2.1 by induction that for any l ∈ N and any sufficiently regular function f , (6.55) Assuming that D log N H k is uniformly bounded, the four commutator terms can be estimated bŷ (6.56) The statement now follows from combining all the estimates and using (6.44).
LEMMA 6.13. Let u be a function,ū = ffl M udV and u ⊥ = u −ū. Then we have (6.57) In particular, under the assumptions of Proposition 6.12, Proof. Recall that ∂ t g = −2Π − L X g and e 0 = N −1 (∂ T + X ). Then a straightforward computation shows (6.59) The second assertion of the lemma follows from standard estimates and using Π = −Σ+(N −1 −3 −1 )g .
PROPOSITION 6.14. Under the assumptions of Lemma 6.12, we have Proof. This follows from Proposition 6.12, Lemma 6.13, and using the notations S = div g Π − Dtr g Π and Π = −Σ + (N −1 − 3 −1 )g . REMARK 6.15. Note that we have control over Φ −1 H k−1 due to Lemma 6.13 and the exponential decay of the energy that we will obtain. 7. PROOF OF THE MAIN THEOREM 7.1. Preliminaries and local existence. Small perturbations of an initial data set, corresponding to the background solution, are not necessarily CMC. As argued in [FK15], for a related situation, the corresponding maximal globally hyperbolic development contains a CMC surface with data close to the background. A similar argument applies in the present context. Starting from this CMC surface we apply the local existence theory for the reduced system, which is of hyperbolic-elliptic nature. An analysis as in [AM-1] yields a local existence theory for our system and a continuation criterion assuring the existence as long as the Sobolev norms in suitable regularity (H 4 for metric and fields and H 3 for time derivatives) is sufficient. It therefore suffices to establish the energy decay to conclude global existence. 7.2. Global existence. The proof of Theorem 1.1 is an almost immediate consequence of the individual energy estimates for the geometry, the one-forms and the functions. We define a total energy measuring all perturbations simultaneously by (7.1) E tot (g − γ, Σ, ω, Φ) := E 4 (g − γ, Σ) + e −2T E 4 (ω) + E 4 (Φ).
The energy-estimate for the total energy is given by LEMMA 7.1. Under smallness assumption on the perturbation the following estimate holds.
Proof. This estimate is a consequence of Propositions 6.6, 6.10 and 6.14 and the elliptic estimates.
This yields a decay of the total energy of the rate e −2αT , which is used in the following to determine the decay of the individual energies based on the individual energy estimates.
7.3. Decay rates. Appealing to the individual energy estimate for E 4 (ω) this yields the following decay rates.
LEMMA 7.2. For sufficiently small initial perturbations, the following estimates hold.

OTHER RELATED SYSTEMS
In this section we list other well-known models for which our method applies. More precisely, one also obtains nonlinear stability of the 3 + 1-dimensional Milne model as a solution of the following systems.
8.1. The Brans-Dicke system. The Brans-Dicke model is governed by the action where ω denotes the dimensionless coupling constant. The critical points of the Brans-Dicke action are given by This system is obtained from our result by setting F = 0 and assuming T q = S 1 so that Φ mn = φ.
8.2. The Einstein-Wave map system. In order to define the Einstein-Wave map system we also take into account a Riemannian manifold (P, k i j ) and consider a map φ : M → P . This allows us to provide the action for Einstein-wave maps The critical points of the Einstein-wave map system are given by where Γ i j k (φ) are the Christoffel symbols on the Riemannian manifold P . This system is a slight modification of our system in the case of F = 0 since we are now assuming that the map φ takes its values in a Riemannian manifold. If the map is almost constant one can assume that its image is contained in a single coordinate chart such that we can think of it as a set of functions rather than a map between manifolds. The Einstein-wave map system is not directly captured by our main result by setting F = 0 but exactly the same energies can be used.
8.3. The Einstein-Maxwell system. To define the energy for the Einstein-Maxwell system we consider the vector potential A µ and its curvature two-form F µν . The energy functional for the Einstein-Maxwell system is the following The critical points of the Einstein-Maxwell system are given by (8.6) We obtain the Einstein-Maxwell system by setting Φ pq = 0, neglecting the equation for Φ pq and changing some constants on the right hand side of the equation on the metric.