Nonlinear stability of the Milne model with matter

We show that any 3+1-dimensional Milne model is future nonlinearly, asymptotically stable in the set of solutions to the Einstein-Vlasov system. For the analysis of the Einstein equations we use the constant-mean-curvature-spatial-harmonic gauge. For the distribution function the proof makes use of geometric L^2-estimates based on the Sasaki-metric. The resulting estimates on the energy momentum tensor are then upgraded by employing the natural continuity equation for the energy density.

1. Introduction 1.1. Cosmological spacetimes and stability. We consider the following class of cosmological vacuum spacetimes. Let the M be a closed 3-manifold admitting an Einstein metric γ with negative Einstein constant µ = − 2 9 , i.e. (1.1) where the specific value of µ is chosen for convenience. A spacetime of the form ((0, ∞) × M, g) with (1.2) g = −dt 2 + t 2 9 · γ is known as a Milne model and is a solution to the vacuum Einstein equations. Its future nonlinear stability under the vacuum Einstein flow has been shown in  and constitutes the second stability result for the vacuum Einstein equations without symmetry assumptions beside the corresponding one for Minkowski space [CK]. While the stability of the Minkowski spacetime under the vacuum Einstein flow has been generalized to several Einstein-matter systems [BZ,LR,LM,T,LT17] this is not the case for the Milne model. We address this problem for the Einstein-Vlasov system.
1.2. The stability problem for the Einstein-Vlasov system. The Einstein-Vlasov system (EVS) reads where X g denotes the geodesic spray and f a distribution function with domain P x ⊂ T M , the mass-shell of future directed particles for a fixed mass m. It models spacetimes containing ensembles of self-gravitating, collisionless particles and constitutes an almost accurate model for spacetime on large scales, where collisions are negligible and galaxies and galaxy clusters indeed interact solely be their mutual self-gravitation. Its mathematical study in the context of the Cauchy problem dates back to the first works by Rein and Rendall on the evolution of spherically-symmetric perturbations of Minkowski space [RR] and the construction of static nonvacuum solutions . Substantial progress in the study of the EVS happened since then. For a complete overview we refer to the review article by Andréasson [A]. Regarding the nonlinear stability problem, in particular without symmetry assumptions, first results have appeared recently considering different geometric scenarios. Ringström's monumental work, which in particular contains a detailed local-existence theory, addresses the stability problem for exponentially expanding cosmological models [Ri]. These correspond to the presence of a positive cosmological constant in the Einstein equations, which in his case is realized by a scalar field with suitable potential. This has later been extended by Andréasson and Ringström to prove stability of T 3 Gowdy symmetric solutions (in the class of all solutions without symmetry assumptions) [AR]. Furthermore, the stability of Minkowski space for the Einstein-Vlasov system for massless particles has been proven by Taylor [T]. The stability of 2+1-dimensional cosmological spacetimes for the Einstein-Vlasov system has been proven by the second author [F-1, F-2]. We remark that in the physically interesting case of 3+1 dimensions, nonlinear stability results until very recently either required a positive cosmological constant or a restriction to the massless case. A recent series of works then established the stability of Minkowski space for the Einstein-Vlasov system by a vector-field-method approach [FJS15,FJS17, and also independently [LT17].
In the present paper is we establish the first stability result for the Einstein-Vlasov system in 3+1 dimensions in the cosmological case with vanishing cosmological constant. Moreover, to our knowledge, the present work presents the first stability result to an Einstein-matter system with vanishing cosmological constant in the cosmological case. Further stability results for cosmological spacetimes with matter models exist but to our knowledge consider the case of a positive cosmological constant. We refer here to the works of Rodnianski-Speck and Speck on the Einstein-Euler system [RS13,S12], Hadžić-Speck on the Einstein-dust system [HS15], Friedrich on the Einstein-dust system [Fr17] and Olyniyk on the Einstein-fluid system [Ol16].
1.3. Nonvacuum stability of the Milne model -Main theorem. To prove nonlinear stability of any Milne model within the class of solutions to the Einstein-Vlasov system we first extend the rescaling of the geometry by the mean curvature function as done in  to the nonvacuum case by rescaling the momentum variablesp accordingly. The choice of rescaling here is motivated by the behavior of the momentum support for solutions to the transport equation on the background (1.2), which decreases asp ≈ t −2 . The mass-shell relation of massive particles, however, prevents from obtaining a system of autonomous equations, as it occurs for the vacuum system. In the present case, some explicit time functions remain in the rescaled equations, which appear in conjunction with the energy-momentum tensor. We then combine the technique of corrected energies to control the perturbation of the geometry as developed for the vacuum case in  with the technique of L 2 -Sobolev-energies for the distribution function based on the Sasaki metric on the spatial tangent bundle derived in [F-1].
1.3.1. Rescaling. As in [AM-2] we use here a rescaling of the geometric variables (and in addition of the matter quantities) in terms of the mean curvature τ . This rescaling is introduced in (2.8). Moreover, a logarithmic time variable T is then introduced in (2.10). The following discussion and statement of the main theorem is conducted with respect to these variables.
1.3.2. Difficulties in 3+1 dimensions. A fundamental difference to the 2+1-dimensional case considered in [F-1] is the different structure of the matter quantity appearing in the elliptic equation for the lapse function (cf. τ η in equation (2.14)). In dimension 3+1, after the appropriate rescaling, we find that this quantity does not decay faster than e −T . This occurs already on the level of the unperturbed background geometry and implies that Sobolev norms of the gradient of the lapse function only decay as e −T . In view of this slow decay a critical problem arises when the L 2 -Sobolev estimates for the distribution function are considered. In the transport equation the critical term reads (1.4) e T · N ∇ a N p 0 ∂f ∂p a , written in rescaled variables. Roughly analyzed 1 , the decay of the lapse, of the form ∇N ≈ εe −T , where ε denotes the smallness of the initial perturbation, then leads to a small growth of the L 2 -Sobolev energy of the distribution function as e εT . The problem is then apparent if this growth of the matter perturbation couples back into the lapse equation, where it reduces the decay of the gradient of the lapse to εe (−1+ε)T . This cannot be closed in the sense of a suitable bootstrap argument or by an appropriate energy estimate. A correction mechanism for the L 2 -Sobolev energy of the distribution function as used to deal with problematic shift vector terms in the 2+1-dimensional case in [F-1] seems unavailable as the critical terms here do not necessarily appear as an explicit time derivative, which allowed for the correction in [F-1].
1.3.3. A new estimate for the energy density. We resolve the problem of slow decay of the lapse gradient by a different idea. A crucial observation therefore is the fact that the matter term in the lapse equation decomposes as in rescaled variables, where ρ is the rescaled energy density (cf. (2.20)). For Vlasov matter, the remaining term τ 2 η has stronger decay properties due to the explicit τ variable, which can be used to compensate a growth of the L 2 -Sobolev energy of the distribution function. This implies that accepting a small growth for the L 2 -Sobolev energy still yields a decay of τ 2 η ≤ εe (−3+ε)T , which is sufficiently fast. The problematic term is in fact the rescaled energy density ρ. The crucial idea is not to estimate the energy density by the L 2 -Sobolev energy of the distribution function but to use an explicit evolution equation for ρ, which originates from the divergence identity of the energy momentum tensor, ∇ µ T µν = 0. One obtains the evolution equation for the energy density or continuity equation, which in rescaled form (cf. Appendix B) reads where in the setting we consider the last three terms have improved decay from the additional τ factors. This seems to be a particular feature of massive collisionless matter but this structure may also be relevant for other massive matter models. If those terms are estimated by the L 2 -Sobolev energies this additional decay can be used to compensate for the small growth and yields a uniform estimate for the standard Sobolev norm of ρ without the problematic loss. This mechanism allows to close the estimates. It is important to remark that the regularity loss of the evolution equation (1.6) for ρ is compensated by the elliptic regularity of the lapse equation which requires the energy density only at one order of regularity below the top order.
1.3.4. Structure of the proof. The small growth of the L 2 -Sobolev energy of the distribution function, which results from the lapse term, implies that we do not correct this energy as done in [F-1], where we required uniform boundedness. The corresponding energy estimates here are done with respect to the rescaled variables and require higher orders of regularity but are except for these aspects similar to the ones in [F-1]. Also similarly to [F-1] we consider initial data with compact momentum support. We expect that considering non-compact momentum support results in similar decay properties of the system. However, to analyze this issue in detail another additional structural estimate for the transport equation is necessary which is subject to future works on the topic. Regarding the estimates for the perturbation of the geometry we use energy estimates and elliptic estimates according to the vacuum case , where in the present case additional terms due to the matter quanitities appear. For the sake of brevity we derive most estimates under smallness assumptions on the perturbation, which allows us to suppress higher order terms in the perturbation in the estimates and absorb them into uniform constants. Global existence is eventually shown by a bootstrap argument, which implies that for a sufficiently small initial perturbation the smallness assumptions persist throughout the evolution and almost optimal decay holds, if we compare with the vacuum case.
Main theorem. We formulate the main theorem using the terminology of the remainder of the manuscript. The theorem is formulated with respect to the rescaled metric and second fundamental form. After the theorem we clarify the notation used therein.
Theorem 1. Let (M, γ) be a 3-dimensional, compact, Einstein manifold without boundary with Einstein constant µ = − 2 9 and ε decay > 0. Then there exists an ε > 0 such that the future development of the rescaled initial data under the Einstein-Vlasov system is future complete and the rescaled metric and tracefree fundamental form (g, Σ) converge as with decay rates determined by ε decay as in (10.10). In particular, any 3+1-dimensional Milne model is future asymptotically stable for the Einstein-Vlasov system in the class of initial data given above.
The symbols (g, k, Σ, f ) denote the Riemannian metric, the second fundamental form, the tracefree part of k and the distribution function, respectively. τ < 0 is the mean curvature and is related to the time variable in (1.2) via t = −3τ −1 with τ ր 0 being the future direction. B 6,5,5 ε (. , . , .) denotes the ball of radius ε centered at the argument in the set of H 6 (M ) × H 5 (M ) × H Vl,5,3 (T M ) with the canonical Sobolev norms defined further below. Here, H Vl,5,3 denotes the space of distribution functions on T M corresponding to the standard L 2 -Sobolev norms, cf. [F]. H Vl,5,3,c (T M ) is the subset of this space with distribution functions of compact momentum support.
1.4. Remarks. The decay rates (10.10) can be achieved for arbitrarily small ε decay by choosing the perturbation sufficiently small depending on ε decay . This implies that one can get arbitrarily close to the vacuum decay rates which correspond to the case ε decay = 0. The corresponding higher dimensional stability results, which for the vacuum equations have been considered in , are likely resolvable similarly to the case presented herein. In particular, the decay of the matter quantities is expected to be stronger than in the present case. In this sense the 3+1-dimensional case is more difficult.
1.5. Overview on the paper. The remainder of the paper is concerned with the proof of Theorem 1. To simplify the presentation we derive all estimates -hyperbolic and elliptic ones -under smallness assumptions on the solution. These smallness assumptions are compatible with the decay properties of the system and this consistency is then shown in the course of a bootstrap argument. In section 2 we discuss the eigenvalue estimate for the Einstein operator for 3-dimensional negative Einstein metrics, recall the rescaling for the Einstein equations and introduce the rescaling for the matter variables. All relevant equations are collected in section 2 and referred to in the course of the following sections. In section 3 we introduce all relevant norms for the geometric quantities and for the distribution function. In view of these, we introduce the notion of smallness which is a prerequisite for establishing all estimates to follow in their respective concise versions. In the global existence argument this notion of smallness is realized in terms of a suitable bootstrap assumption (cf. (10.3)). In sections 4 and 5 we prove the L 2 -energy estimate and the evolutionary inequality for the bound on the momentum support, respectively. In section 6 we derive the direct energy estimate for the standard Sobolev norm of the energy density ρ of the distribution function. In section 7 we prove elliptic estimates for lapse and shift and their time derivatives. Section 8 contains the energy estimate for the perturbation of the metric and the tracefree part of the second fundamental form. In section 9 we use the elliptic estimates to reduce all evolutionary estimates to a system of estimates solely containing metric, second fundamental form and matter quantities. Basing on these estimates section 10 presents the proof on Theorem 1, which also contains a number of technical remarks on local existence and existence of initial data in the appropriate sense. The appendix contains a collection of formulae used throughout the paper.
Acknowledgements. This project results from early discussions of the authors during the conference Complex Analysis and Dynamical Systems in Acre 2011. D. F. is grateful to Klaus Kröncke for discussions on his results on eigenvalues of the Einstein operator in [Kr15]. D. F. furthermore gratefully acknowledges the support of the Austrian Science Fund (FWF) through the START-Project Y963-N35 of Michael Eichmair as well as through the Project Geometric transport equations and the non-vacuum Einstein flow (P 29900-N27). D.F. acknowledges the hospitality of the Erwin-Schrödinger Institute Vienna during the program Geometric Transport equations in General Relativity.

Preliminaries
We fix for the remainder of the paper a 3-dimensional Einstein manifold (M, γ) with We consider the Einstein operator associated with γ, where • Rh ij = R ikjl h kl for symmetric 2-tensors h (cf. Chapter12D of [B] for more details). The lowest positive eigenvalue of ∆ E plays a crucial role for the construction of suitably decaying energies in the stability problem for the vacuum Einstein flow as demonstrated in [AM-2]. A similar consideration will be relevant for the nonvacuum problem considered below. We denote the lowest positive eigenvalue of ∆ E by λ 0 . The following is an immediate consequence of Kröncke's lower bound on eigenvalues of the Einstein operator (cf. [Kr15]).
Proposition 2. Let (M, γ) be a hyperbolic Einstein 3-manifold with Einstein constant µ = −2/9. Then Proof. From Proposition 3.2 [Kr15] we deduce that the smallest eigenvalue of ∆ E T T , ∆ E restricted to TT-tensors on (M, γ), which we denote by λ 0,T T , obeys This holds, as γ is necessarily of constant scalar curvature and therefore has vanishing Weyl tensor. We show that this can be upgraded to (2.4) as follows. We observe that if an eigenvalue λ of ∆ E obeys (2µ + λ) < 0 or with the present choice λ < 4 9 , then its corresponding eigentensor h λ is TT. This follows as in the proof of Lemma 2.7 in . In particular, the lowest eigenvalue λ 0 either fulfills λ 0 ≥ 4/9 or is in the spectrum of ∆ E T T and in turn fulfills λ 0 ≥ 1/9. A relevant corollary of the above reads Corollary 3. Let (M, g) be a 3-dimensional Einstein manifold with Einstein constant µ = −2/9, then This condition assures that the energy to control the perturbation of the geometry defined below is coercive and allows to avoid introducing a shadow-gauge analog to [AM-2].
2.2. Variables and setup. We use standard index conventions. Roman letters denote spatial indices {1, 2, 3} and greek letters denote spacetime indices {0, 1, 2, 3}. In addition, we use bold roman letters to denote indices on the tangent bundle of T M . This notation is introduced in Section 4.

Standard variables and gauge.
We consider the 3+1-dimensional spacetime in the standard where N , g and X denote the lapse function, the induced Riemannian metric on M and the shift vector field 2 . For the derivation of the Einstein equations in ADM formalism we refer to [Re]. We denote by τ the trace of the second fundamental form k with respect to g and decomposek =Σ + τ 3g . We then impose the CMCSH gauge via where Γ and Γ denote the Christoffel symbols of g and γ, respectively.

Rescaled variables and Einstein's equations.
We rescale the geometry with respect to the mean curvature function τ analogous to the vacuum case . This leaves explicit time-factors as coefficients of the matter variables. We rescale those by rescaling thep-variables (cf. section 2.3 ). The variables with respect to mean curvature time t = τ are denoted by ( g, Σ, N , X), while the rescaled variables are (g, Σ, N, X). We rescale according to (2.8) Then we introduce the logarithmic time We use the notationẊ = ∂ T X,Ṅ = ∂ T N for convenience throughout the manuscript. Also, we denote N = N 3 − 1 and X = X/N . After these modifications the Einstein equations in CMCSH gauge with respect to the rescaled variables take the following form. (2.14) We denote by L X the Lie-derivative with respect to X. Moreover, we recall the decomposition of the curvature term in the spatial harmonic gauge (cf. [AM-2]), The rescaled matter quantities are connected to the unrescaled versions via (2.20) We recall thatρ =Ñ 2 T 00 is the energy density and a = −Ñ T a 0 is the matter current. We also denote T ab = T ab = |τ | −7T ab for later purposes. An important identity, which follows immediately from the definitions above is The decomposition is crucial since the second term on the right-hand side in (2.21) decays fast while the first term is handled differently using the continuity equation as explained in the introduction.
Remark 4. The right-hand sides of the elliptic system for lapse and shift as well as those for the evolution equations decouple into principal terms and terms which can be considered as perturbative and which turn out to decay faster than the leading order terms. To give some orientation about which terms are considered to be principal terms, we have marked those terms by the symbol (⋆) or (⋆⋆). The latter case refers to those terms, which are relevant to establish the decay for the energy measuring the perturbation of the geometry. Terms denoted by (⋆) are for different reasons principal. In the lapse equation, the ρ-term within the η-term has the slowest decay, while in the shift equation, precisely this slow decay is inherited from the lapse equation through the (⋆) term therein. The final principal term to consider is the one in the equation for Σ, where due to regularity conditions we cannot estimate the ρ-term in S by the ρ-energy but we have to use the L 2 -Sobolev energy of the distribution function to estimate this term. This results in a small loss of decay, which is the reason why this term is of worst decay in the respective equation.
2.3. Vlasov matter. We introduce the structures relevant to Vlasov matter and then rescale the energy-momentum tensor and transport equation according to the previous section.
2.3.1. The mass-shell relation. We consider particles of positive mass m = 1 modeled by distribution functions with domain being the mass-shell In particular,p α are canonical coordinates on the tangent bundle. We use the .-notation, since below we introduce rescaled variables. A distribution function f : P → [0, ∞) has the associated energy-momentum tensor where µ Px is the volume form corresponding to the induced metric on P x . We consider the projection of the distribution function under π : and which we refer to as distribution function in the following. We rescale the momentum variables according to Then we express the unrescaled mass-shell relation in (2.22) in coordinates (cf. for instance Section iv in [SZ], equation (37)), which reads and replace all variables by their rescaled counterparts. This yields an expression for p 0 as a function of the p a variables and the metric components in the form An alternative expression is given by is just defined for convenience and does not necessarily have a specific geometric meaning. In addition, p 0 = g 0νp ν = −N p. We derive some useful estimates for p 0 using elementary manipulations. We furthermore use the simplifying notation (2.29) p = N p 0 .
Remark 5. Note, that the rescaled mass-shell relation (2.26) reduces to p 0 = 1 + τ 2 |p| 2 g , when X = 0, N = 1, which corresponds to the background solution. In particular, the constant term under the squareroot, which originates from the mass term, scales like a constant, while the second term decays fast in expanding direction (τ ր 0).
The following lemma contains two useful pointwise estimates on the momentum variable p 0 .
The transport equation. We introduce the transport equation and its rescaling. The transport equation is rescaled via (2.24). To express the transport equation only in terms of the rescaled variables, we require the rescaled Christoffel symbols Γ. The non-rescaled Christoffel symbols read (cf. [Re]) In terms of the rescaled variables the Christoffel symbols are of the form We refer to the latter terms also by the symbols Γ * and Γ * * , respectively, when the indices are suppressed. The fully rescaled transport equation then reads (2.38) which correspond to the natural horizontal and vertical derivatives on T M . The two marked terms are leading order in the sense that term (⋆), among the small terms, has the slowest decay as Γ a contains in particular ∇N , which in combination with τ −1 is of the order of ε. Term (⋆⋆) is the dilution term, driving the downscaling of the momentum support in expanding direction of spacetime and thereby the dilution of the matter variables.
2.4. Energy momentum tensor. The rescaled matter quantities as appearing in the Einstein equations take the following form in terms of the distribution function f .
Remark 7. Recall the definition of the rescaled energy-momentum tensor (2.20) for the identities above. Note, that the expressions for the energy-momentum variables are a consequence of their definition (2.20) and the rescalings (2.8) and (2.24).
2.5. Preview on the decay rates. For the study of the estimates to follow it is important to have an idea about smallness and decay of the rescaled quantities. The quantities (2.45) N − 3, Σ, g − γ and X are small and decay with the rates For the matter terms we have These decay rates are shown to be valid for sufficiently small initial data, where ε is the smallness of the initial perturbation and ε decay > 0 can be chosen arbitrarily small.

Norms and smallness
We introduce all relevant norms for measuring the perturbation of the geometry and the distribution function. Some norms are defined with respect to the fixed Einstein metric γ and others are defined with respect to the rescaled dynamical norm g. As we impose a uniform smallness assumption all these norms are equivalent. We assume for the remainder of the paper that T 0 > 1.
3.1. Constants. We use the symbol C to denote any positive constant, which is uniform in the sense that it does not depend on the solution of the system once a smallness parameter ε for the initial data and an initial time T 0 are chosen. Furthermore, if ε is further decreased or T 0 is increased, C keeps its value.
3.2. Norms -tensor fields. For functions and symmetric tensor fields on M we denote the standard Sobolev norm with respect to the fixed metric γ of order ℓ ≥ 0 by . H ℓ . The corresponding function spaces are denoted by H ℓ = H ℓ (M ).
3.3. Norms -distribution function. We introduce different metrics on T M and related notation necessary for the definition of L 2 -Sobolev energies for the distribution function. This construction is based on the the metric γ on M . In the following section we consider the case when the corresponding construction is based on the rescaled metric g.
The metric γ induces the related Sasaki metric γ on T M via where Dp i = dp i + Γ i jk p j dx k . Recall that Γ denotes the Christoffel symbols of γ. The covariant derivative on the tangent bundle corresponding to γ is denoted by γ ∇. We consider the volume form on T M , We define a weighted version of the Sasaki metric by where we denote p γ = 1 + |p| 2 γ . This metric is necessary to take the norm in the energies to be defined below, which require a weight in the momentum-direction. We define the L 2 -Sobolev energy of the distribution function with respect to Sasaki metric corresponding to the fixed metric γ by The corresponding function spaces are denoted by H Vl,ℓ,µ (T M ). Pointwise estimates are taken with respect to the following L ∞ x L 2 p -norm, which obeys the following lemma.

Smallness.
We define a set of smallness conditions for the dynamical quantities. These are designed to include weights in terms of the time-function to incorporate some decay properties indirectly. These are chosen in a way that in the proof of global existence the smallness conditions serves as a part of the bootstrap assumptions and leaves room to be improved for sufficiently small data and sufficiently large times. We define (3.7) We say a triple (g(τ ), Σ(τ ), f (τ )) is δ-small when (g, Σ, f ) ∈ B 6,5,5 δ,τ (γ, 0, 0). Also, we mostly suppress the dependence on (γ, 0, 0) in the notation. In addition we use the term smallness assumptions if we refer to δ-small data.
3.5. Some immediate estimates. Smallness in the above sense implies smallness of the perturbation for lapse function and shift vector. The following corollary uses the elliptic estimates proven in section 7.
Corollary 9. For any δ > 0 there exists a δ such that Proof. This is an immediate consequence of Proposition 17.

L 2 -estimates for the distribution function
We define the L 2 -Sobolev energy for the distribution function in terms of the Sasaki metric associated with g. Under the present smallness assumptions this energy is equivalent to the norm |||f ||| ℓ,µ . We define the corresponding metrics on T M with respect to g as follows. The metric g induces the related Sasaki metric g on T M via where Dp i = dp i + Γ i jk p j dx k . The covariant derivative corresponding to g is denoted by ∇. We consider the volume form on T M , (4.2) µ g = |g|dx 3 ∧ dp 3 .
We define a weighted version of the Sasaki metric by where we denote p = 1 + |p| 2 g . For explicit computations including the Sasaki metric on the tangent bundle we use indices a, b, . . . ∈ {1, . . . , 6}, where 1, 2, 3 correspond to horizontal directions and 4, 5, 6 to vertical directions. We introduce the frame {θ a } a≤6 = {A 1 , A 2 , A 3 , B 1 , B 2 , B 3 } and we denote the connection coefficients of the Sasaki metric in this frame by Γ (cf. (10.22)). We define the L 2 -Sobolev energy of the distribution function by Remark 10. The choice for the weights, to increase with decreasing level of regularity, is necessary to absorb terms with high weights, which result from the connection coefficients of the Sasaki metric, where momentum weights appear in conjunction with horizontal derivatives. This is discussed in more detail below.
The main energy estimate for the distribution function is given in the following.
Proposition 12. For δ > 0 sufficiently small, (g, Σ, f ) ∈ B 6,5,5 τ,δ and f a solution to the transport equation (2.38) the following estimate holds. (4.5) Proof. The derivation of the energy estimate is a straightforward and technical computation. It follows the lines of the analogous computations in ( [F-1], section 5). We discuss some exemplary steps, which are more general herein.
We take the time derivative of the square of the energy, which yields four leading order terms. The first term results from the time-derivatives hitting the distribution function and reads (4.6) 2 The second leading order term arises from the time-derivative of the volume form, which, when the derivative acts on the time-function in the rescaled momentum variables and reads (4.7) 6 · E 2 µ,ℓ (f ), since µ T M = |g|τ −6 d 3p , wherep is time-independent. The third leading order term occurs when the time derivative hits the time function in the momentum-weight factor and yields The fourth leading order term results from the time-derivative hitting the momentum variable in the inverse g a i b i , when a i , b i ≥ 4 and reads (4.9) for each pair with a i , b i ≥ 4. The non-explicitly listed terms arise when the time-derivative hits the rescaled metric g, which yields terms of the first three types listed in (4.5). We evaluate term (4.6) in the following. (4.10) We first analyze the terms containing commutators of ∂ T with θ a i . Since θ a i is not affected from the rescaling when a i ≤ 3 (all time-factors cancel) the commutator can be estimated by terms arising from ∂ T Γ(g), which yields terms of the form ∇Σ, ∇ 2 X and ∇(N − 3). When a i ≥ 4, i.e. θ a i = B a i −3 we have [∂ T , θ a i ] = −2θ a i . This implies that in the case of a i ≤ 3, term (1) can be estimated by terms included in the first three terms on the right-hand side of (4.5). In the complementary case this results in a term of the form (4.11) − 2θ a 1 ∇ a 2 . . . ∇ a k f, which requires to be canceled for the estimate to hold, as we see below. The terms (3) and (5) again give rise to terms of the form of the three first terms on the right-hand side of (4.5) if a j+1 ≤ 3 and a k ≤ 3. In the complementary case, terms of the form occur, which are canceled by terms arising below. Regarding terms (2) and (4), from (10.22) we observe that these terms yield time derivatives of Γ or of the Riemann tensor Riem, which in combination again yields terms of the form of the first three terms on the right-hand side of (4.5) and terms that arise when the time derivative hits the rescaled momentum variable in the respective cases in (10.22). From these we again obtain leading order terms, which are of the form (4.13) − 2 Γ e a i a 1 ∇ a 2 . . . ∇ e . . . ∇ a k f when e ≤ 3 and 4 < a i + a 1 ≤ 10 or when e ≥ 4 and a i , a 1 ≤ 3; and (4.14) − 2∇ a 1 . . . ∇ a j Γ e a i a j+1 ∇ a j+2 . . . ∇ e . . . ∇ a k f, when e ≤ 3 and 4 < a i + a j+1 ≤ 10 or when e ≥ 4 and a i , a j+1 ≤ 3. Both types of terms are cancelled by terms arising below.
It remains to consider term (6), where ∂ T f is replaced by the transport equation yielding the following term. (4.15) We begin with the most important term to evaluate, which is here marked by (⋆). This term is relevant for the cancellation of all non-perturbative terms above. Before we start the computation we derive a few simple commutators. The following identities hold. (4.16) We evaluate now the term from above.
According to the commutators above, the second term in the previous line vanishes if a k ≤ 3 or cancels the second term in (4.12). We proceed with the first term. (4.18) According to the previous step, the second term cancels the corresponding term from (4.12) and the third term on the right-hand side cancels the corresponding term from (4.14).
Continuing with the first term on the right-hand side of the previous equation and further commuting p i B i to the left, we obtain terms cancelling all terms in (4.12) and (4.14). Then we are left with the term Integration by parts yields three types of terms. The first term arises when B i acts on p i and cancels (4.7). The second term results from B i acting on p and cancels (4.8). Finally, the term arising from B i acting on g a i b i when a i , b i ≥ 4 cancels (4.9).
It remains to consider the remaining terms in (4.15). When estimating the term corresponding to the first term in (4.15) we use the estimate (4.20) |p| g p ≤ G .
The corresponding term in the estimate (4.5) is |τ |G . Note that compact support is necessary for this. Otherwise we would obtain an additional factor |τ | −1 , which would leave no decay for this term. To outline the estimates in more detail we consider one particular term from (4.15) and claim the other terms can be handled in a similar way. We sketch (4.21) where we suppress all mixed terms. Commuting the operator p u Γ e u B e to the front we obtain a term of the form . The corresponding integral, after an integration by parts, yields the term Γ * * H ℓ in (4.5). The remaining terms, after commuting p u Γ e u B e to the front, are schematically of the form Note that the momentum variables in front of the Riemann tensor, which arises as part of the Γ terms, can increase while appearing as coefficients of ∇ i with i ≤ 3. In this case, the weights in the energy, appearing for lower numbers of derivatives, allow for these terms to be estimated by the energy. All the remaining terms arising from (4.15) can be estimated similarly.

Estimating the energy momentum tensor.
Lemma 13. Under smallness assumptions and ℓ ≥ 4 the following estimates hold.
Proof. We begin by estimating an integral of the form T M F · G(|p| g )µ T M for functions G, F on T M to explain the number of momentum weights. Let µ ≥ 2, then (4.25) Depending on the additional momentum factors in G, which are of order one for ρ and  and two for the other quantities, this explains the order of weights, necessary in the energies. In the above computation F represents the term where derivatives have acted on the distribution function and other quantities in the matter variables. We discuss how to estimate these terms in the following. Covariant derivatives of matter quantities correspond to horizontal derivatives under the momentumintegral by the following identity, Similar identities hold, for f replaced by f p 0 etc. and for higher derivatives. For higher derivatives, we obtain not the full covariant derivative of the Sasaki metric. The additional terms arising from the Riemann tensor in (10.22) can however be added and substracted where the additional terms are lower order and due to the smallness condition, can be absorbed into the constants. Finally, if the horizontal derivative hits the momentum variables such as p or p 0 we use the formulae and estimate the arising shift vector terms using the smallness condition by the constants.

Control of the momentum support
Using the characteristic system associated with the rescaled transport equation we derive an estimate on the supremum of the outer radius of the support of the distribution in momentum space.
The characteristic system corresponding to the rescaled transport equation (2.38) reads We define the auxiliary quantity (5.2) G(T, x, p) ≡ |p| 2 g . Using the characteristic system we compute the derivative of G along a given characteristic. This yields where it is important to recall that the rescaled momentum variables are time-dependent. Invoking the corresponding estimates for p 0 and |p| g (p 0 ) −1 in (2.30), (2.31) we deduce the following estimate.
Lemma 14. Under smallness assumptions, the following estimate holds for any characteristic.
We define the supremum of the values of G in the support of f at a fixed time T by From the estimate for individual characteristics above, we derive an estimate for G , which serves as a bound for momenta in the support of the distribution function.
Proposition 15. Under smallness assumptions we obtain Proof. For any characteristic in the support of f we obtain an inequality of the form

Energy estimates from the divergence identity
The key quantity, which provides improved estimates for the energy density is the standard L 2 -Sobolev energy for the rescaled energy density ρ with respect to the dynamical metric g on (M, g). This energy reads We derive the energy estimate for ̺ ℓ in the following. We denote (6.2) ∂ 0 ≡ ∂ T + L X as this combination of derivatives naturally appears when taking the time-derivative of norms, which are taken w.r.t. the volume form µ g (cf. below). A part of the divergence-identity for the energy momentum tensor with 0-component reads in its rescaled form (cf. (10.16)). Two identities relevant for the energy estimate in rescaled form are (6.4) ∂ 0 g ab = −2N Σ ab + 2(1 − N/3)g ab . and (6.5) for a function u on M (cf. [CM]). Moreover, for a function u. This identity arises from the corresponding unrescaled one by multiplication with −τ . Next, we derive the standard energy estimate for this energy.
Proposition 16. Let ℓ ≥ 4 (6.7) Proof. We take the time derivative of one of the summands of the square of the energy, which takes the following form. (6.8) The first three terms on the right-hand side contribute to the first line of the estimate. We proceed with the evaluation of the term I.
Using the commutator formula above we obtain (6.9) where we use the notation (6.10) (N k cb )] . The second term on the right-hand side can be estimated by terms of the form (6.11) which yield the third term in the estimate. We continue with estimating the final term II using (6.3). (6.12)

Elliptic estimates
We derive in this section elliptic estimates on the lapse function, the shift vector and their respective time derivatives.
Proposition 17. Under smallness conditions, for the lapse function, a pointwise estimate of the form 0 < N ≤ 3 holds and moreover the following two estimates.
Proof. The pointwise estimate for the lapse follows from the lapse equation and the maximum principle. The two following estimates are a straightforward consequence from elliptic regularity applied to the elliptic system for lapse and shift.
Furthermore we require estimates for the time derivatives of of the lapse function and shift vector. These are given in the following lemma.
Lemma 18. The following estimates hold under smallness conditions, for T sufficiently large and ℓ ≥ 4. (7.2) Proof. Both estimates follow from standard elliptic regularity estimates and the elliptic system for (∂ T N, ∂ T X), which is deduced from the elliptic system for (N, X) by taking the derivative with respect to ∂ T . This system reads Here we use ., ., . to denote any suitable contraction of a number of tensor fields, where the specific structure of indices does not matter. Due to the time derivative of η and the terms containing ∂ T N explicitly in the equation for ∂ T X we do the estimates in two steps. Note furthermore, that we do not aim at the sharpest possible estimates and allow rather rough but brief expressions where we absorb many terms into the constants.
From elliptic regularity and equation (7.3) we obtain (7.5) Using the smallness we can absorb the last line of the previous equation into the left-hand side and obtain a formally identical estimate where the last line is not present. The term including the time derivative of ρ is treated using the evolution equation (10.16). This yields Now, we estimate the remaining term using the corresponding formula (10.23). Invoking the smallness assumption and the fact that when taking derivatives of the explicit function of the momentum only yields terms with an additional smallness factor, then reduces the number of relevant terms to an estimate of the following form. (7.7) The term in the second line can be absorbed in the constant in estimate (7.5) by the largeness of T . Before concluding the estimate for ∂ T N we require the estimate for ∂ T X to replace the corresponding terms in the previous estimate. We therefore turn to the equation for ∂ T X and apply elliptic regularity which yields the following first estimate, where we again absorb several terms in the constant due to the smallness criterion. (7.8) The last term on the right-hand side can be absorbed into the constant by the smallness assumption. The term containing the time derivative of  can be estimated using (10.16) by At this point the estimate for ∂ T X is not complete, since there are still ∂ T N terms on the righthand side. We return to the estimate for ∂ T N and absorb the corresponding terms in the estimate and then finish the estimate for ∂ T X. Plugging (7.8) without the last term on the right-hand side into (7.7) and the resulting estimate into (7.5), without the last line on the right-hand side, we observe that every ∂ T N term on the right-hand side comes with a |τ | 3 and consequently can be absorbed into the constant. This proves the estimate for ∂ T N , which in particular is independent of ∂ T X. Then, in turn, plugging the final estimate for ∂ T N into the estimate for ∂ T X and simplifying the estimates with respect to the smallness criteria finishes the proof.
8. Energy estimate -geometry 8.1. Decomposing the evolution equations. We decompose the evolution equations into their principle parts and higher order terms, which are eventually treated as bulk terms.
The evolution equations can be rewritten to the following system.
where the bulk terms obey estimates of the form H s < ε for ε sufficiently small. Using the elliptic estimates for lapse and shift we obtain the following estimates for the bulk terms.
Energy. We define the energy for the tracefree part of the second fundamental form and the metric perturbation below. The choice is identical to the vacuum case considered in [AM-2] and we briefly recall the relevant aspects and point out the improvements in 3+1 dimensions compared to the higher dimensional case. The definition of the energies, which include a correction factor to obtain a suitable decay estimate, depends on the lowest eigenvalue of the Einstein operator corresponding to the specific Einstein metric, λ 0 . Due to the lower bound (2.4) we only distinguish between two cases here. We define the correction constant α = α(λ 0 , δ α ) by (8.4) α = 1 λ 0 > 1/9 1 − δ α λ 0 = 1/9, where δ α = 1 − 9(λ 0 − ε ′ ) with 1 >> ε ′ > 0 remains a variable to be determined in the course of the argument to follow. By fixing ε ′ once and for all, δ α can be made suitable small when necessary. The corresponding correction constant, relevant for defining the corrected energies is defined by We are now ready to define the energy for the geometric perturbation. For m ≥ 1 let Then, the corrected energy for the geometric perturbation is defined by Under the imposed conditions, the energy is coercive.
Lemma 19. There exists a δ > 0 and a constant C > 0 such that for δ-small data (g, Σ, f ) the inequality Proof. The proof is analogous to the corresponding Lemma 7.2 in . The difference consists in the fact that in the 3+1 dimensional setting here, the kernel of the Einstein operator consists only of the zero-tensor (cf. Corollary 3). This implies that the projection operator necessary in Lemma 7.2 [AM-2] is not necessary in the present case.
The energy estimate for the corrected energy is given in the following.
Lemma 20. Under a smallness assumption on E s we have Proof. The proof is analogous to the proof of Lemma 7.6 in [AM-2]. The only difference results from the additional matter term in the evolution equation for Σ. As a direct consequence of the equation, this yields terms of the types (8.10) which can straightforwardly be estimated by yielding the claim.

Total energy estimate
With the individual energy estimates for geometry and matter variables at hand these require to be synchronized in view of their different decay inducing terms. For this purpose we define a total energy with explicit weight functions in time and bound all elliptic variables in terms of this energy. We then derive energy estimates under the smallness assumption on ρ 4 (f ), G and the total energy which are the key estimates to establish the global existence result further below.
9.1. Total energy. We define the total energy including the matter energy and the energy for the metric perturbation.
We choose now all auxiliary constants in the following way. For a given ε decay < 1 we choose positive constants (δ α , δ E , δ E , ε tot ) such that hold. For small ε decay this is achieved, when δ E is almost one and δ E is sufficiently small relative to δ E such that δ E + δ E < 1 holds. We define a uniform constant C, that bounds all constants C in previous estimates from above by (9.3) 10 · C 3 ≤ C.

Preparations.
We gather now a number of simplifying lemmas to reduce the length of the final energy estimate. We express in the following all relevant norms in terms of the energies E s , E 5,4 (f ), ̺ 4 (f ) and G . For the norms appearing in the energy estimate for the L 2 -energies we have Lemma 22. Under suitable smallness assumptions the following estimates hold.
In total, 9.3. Estimates for ̺ 4 (f ). We begin with an estimate for the auxiliary energy of the energy density.
Lemma 23. For δ-small data with δ sufficiently small, the folllwing estimate holds.
Proof. From Proposition 16, using Lemma 13 and Proposition 17, we obtain (9.7) Estimating by the total energy and integrating yields Then, Gronwall's lemma yields the claim.
9.4. Estimate on G . For the bound on the support of the momentum variables we obtain the following estimate.
Lemma 24. For T 0 > 1 and under the δ-smallness assumption for δ sufficiently small, the following estimate holds. (9.9) Proof. The estimate follows directly from Proposition 15 in combination with Lemma 22.
9.5. Estimate -total energy. We proceed with an estimate on the total energy under a smallness assumption on the auxiliary energy.
Proof. Taking the time derivative of the total energy, using the estimate for the energy for the perturbation of the geometry, Lemma 20, and the estimate for the L 2 -Sobolev energy of the distribution function, Proposition 12, we obtain (9.14) .
The terms resulting from the energy estimate for E 6 are denoted by numbers (1.i). The term (1.1) results from the decay inducing term in the estimate (8.9) and the time derivative of the time-weight function. The term (1.2) results from any matter term in the estimate (8.9), where we have to estimate by the L 2 -norm since the regularity is up to the order s − 1 = 5. Note that the time-weight function is distributed to re-obtain the properly weighted energies as they appear in the total energy. Finally, term (1.3) results from the higher order term. The terms resulting from the energy estimate for E 5,4 (f ) are denoted by numbers (2.i). Term (2.1) results from the time derivative of the time-weight function. Term (2.2) bounds all terms from estimate (4.5), which result from the term τ −1 N −1 Γ * , which is estimated using (9.5) where only the term with ρ 4 (f ) is considered, all other terms are of higher order in energy and are absorbed into the term (2.4) except for the term τ G , which is estimated by (2.3).
Using the smallness conditions appropriately, the previous estimate reduces to tot . Here, terms (1.1) and (2.1) appear as before and provide decay inducing terms. Terms (1.3) and (2.4) are absorbed in the higher order term. Invoking smallness conditions (9.10), (9.11) and (9.12) allows us to bound the sum of terms (1.2), (2.2) and (2.3) by ε tot E tot . This yields tot , which under the conditions (9.2) on the auxiliary constants yields the claim.

Global existence and completeness
In this final section we present the proof of Theorem 1 based on the estimates in the previous sections.
10.1. Preliminaries. We consider initial data at time T 0 , which is close to the induced data of the Milne model at T = T 0 . The data is not necessarily CMC initial data. We argue below why it is sufficient to consider only CMC initial data and consider this case for now. The existence of a local-in-time solution for CMC initial data close to the Milne geometry has been developed in [F] and we adapt the local-existence theory therein to our present notation and variables. The local existence theorem (Theorem 4.2, [F]) assures existence of a unique local solutions for initial data (g 0 , k 0 , f 0 ) ∈ H 6 × H 5 × H 5 Vl,3 , which is the regularity assumed in the present case.
Moreover, this solution is depending on the initial data in a continuous sense, which allows to increase T 0 suitably and assume smallness at the increased T 0 without loss of generality. We denote the smallness parameter according to which we express smallness of the initial data in the sense of B 6,5,5 ε 0 by ε 0 . To establish global existence we require a continuation criterion analogous to Theorem 8.1 in [F]. It is important to specify this to our present situation where we consider the rescaled system in 3+1-dimensions. If we replace the non-rescaled system in [F] by the rescaled equations (2.14) -(2.17), the smallness, which has to be assured to continue the solution translates to for a fixed ε loc > 0. This means, either the maximal interval of existence is infinite or the bound above is attained as this time is approached. In particular, starting with sufficiently small initial data, if this smallness persists throughout the evolution, global existence is automatically assured. This persistence is shown for the initial data we consider, which according to the previous discussion guarantees existence of the solution.
10.2. Existence of a CMC surface. Considering sufficiently small initial data which is not necessarily CMC, the maximal globally hyperbolic development under the Einstein-Vlasov system is, locally in time, as close to the background geometry as desired in a suitable regularity [Ri]. The existence of a CMC surface in such a spacetime can be shown along the lines of the corresponding argument in the vacuum case presented for instance in [FK15].
10.3. Guaranteeing the smallness condition on an open interval. From the local Cauchy stability by choosing the initial data sufficiently small we can assure existence of the solution up to T 0 and smallness at T 0 such that condition (9.12) holds at T 0 . We choose the new initial data at T 0 small such that Since all estimates are uniform in the sense that they do not depend on the smallness of the initial data once ε 0 is chosen sufficiently small, we can further decrease ε 0 in the course of the argument. The same holds for increasing T 0 .
By local existence T + > T 0 exists. Note that the condition (9.12) holds automatically at later times.
10.4. Improving the bootstrap conditions -Global existence. We show in the following that if ε 0 > 0 is sufficiently small then T + = ∞.
Using the previous decay result in combination with Lemma 23 yields Choosing now ε 0 sufficiently small, this implies (9.10) with a strict inequality.
Finally, we invoke Lemma 24, which, using the previous estimates, takes the form (10.7) This implies < ε tot /3, by choosing ε 0 sufficiently small. In total, we have shown that for sufficiently small ε 0 the estimates (9.10) and (9.11) hold with strict inequalities on (T 0 , T + ) and In particular, all relevant norms remain sufficiently small as T → T + and by the continuation criterion the solution can be extended to (T 0 , T + + ε) for a small ε, where (9.10) and (9.11) hold on this extended interval. A standard continuity argument then implies T + = ∞.
10.5. Decay and asymptotic stability. From the decay of the total energy the decay rates of the individual quantities can be inferred to read ∂ T g ab = −g ac g bd ∂ T g cd (10.11) Also relevant for time differentiation of energies is the following formula for the time derivative of the Christoffel symbols.
10.6.1. Curvature of the tangent bundle. The connection coefficients Γ of the Sasaki metric g with respect to the connection basis A a = D x a , B a = ∂ p a take the following form.