The Global Nonlinear Stability of Minkowski Space for the Massless Einstein–Vlasov System

Minkowski space is shown to be globally stable as a solution to the Einstein–Vlasov system in the case when all particles have zero mass. The proof proceeds by showing that the matter must be supported in the “wave zone”, and then proving a small data semi-global existence result for the characteristic initial value problem for the massless Einstein–Vlasov system in this region. This relies on weighted estimates for the solution which, for the Vlasov part, are obtained by introducing the Sasaki metric on the mass shell and estimating Jacobi fields with respect to this metric by geometric quantities on the spacetime. The stability of Minkowski space result for the vacuum Einstein equations is then appealed to for the remaining regions.


Introduction
It is of wide interest to understand the global dynamics of isolated self-gravitating systems in general relativity. Without symmetry assumptions, problems of this form present a great challenge even for systems arising from small data. In the vacuum, where no matter is present, the global properties of small data solutions were first understood in the monumental work of Christodoulou-Klainerman [10]. They show that Minkowski space is globally stable to small perturbations of initial data, i.e. the maximal development of an asymptotically flat initial data set for the vacuum Einstein equations which is sufficiently close to that of Minkowski space is geodesically complete, possesses a complete future null infinity and asymptotically approaches Minkowski space in every direction (see also Lindblad-Rodnianski [26], Bieri [2], and also Section 1.15 where these results, along with other related works, are discussed in more detail).
In the presence of matter, progress has been confined to models described by wave equations. 1 Here collisionless matter, described by the Einstein-Vlasov system, is considered. This is a model which has been widely studied in both the physics and mathematics communities; see the review paper of Andréasson [1] for a summary of mathematical work on the system. New mathematical difficulties are present since the governing equations for the matter are now transport equations, though in the case considered here, where the particles have zero mass and hence travel through spacetime along null curves, the decay properties of the function describing the matter are compatible in a nice way with those of the spacetime metric.
The Einstein-Vlasov system takes the form The unknown is a Lorentzian manifold (M, g) together with a particle density function f : P → [0, ∞), defined on a subset P ⊂ T M of the tangent bundle of M called the mass shell. The function f (x, p) describes the density of the matter at x ∈ M with velocity p ∈ P x ⊂ T x M. Here (x μ , p μ ) denote coordinates on the tangent bundle T M with p μ conjugate to x μ , so that (x, p) denotes the point p μ ∂ x μ | x in T M. The Ricci curvature and scalar curvature of (M, g) are denoted Ric, R respectively. The integral in (1) is taken with respect to a natural volume form, defined later in Section 2.2. The vector field X ∈ (T T M) is the geodesic spray, i.e. the generator of the geodesic flow, of (M, g). The Vlasov equation (2) therefore says that, given (x, p) ∈ T M, if γ x, p denotes the unique geodesic in M such that γ x, p (0) = x,γ x, p (0) = p, then f is constant along (γ x, p (s),γ x, p (s)), i.e. f is preserved by the geodesic flow of (M, g). Equation (2) is therefore equivalent to for all s ∈ R such that the above expression is defined, where exp s : T M → T M is the exponential map defined by exp s (x, p) = (γ x, p (s),γ x, p (s)).
In the case considered here, where the collisionless matter has zero mass, f is supported on the mass shell P := {(x, p) ∈ T M | p is null and future directed}, a hypersurface in T M. In this case one sees, by taking the trace of (1), that the scalar curvature R must vanish for any solution of (1)- (2) and the Einstein equations reduce to The main result is the following. Theorem 1.1 Given a smooth asymptotically flat initial data set for the massless Einstein-Vlasov system suitably close to that of Minkowski Space such that the initial particle density function is compactly supported on the mass shell, the resulting maximal development is geodesically complete and possesses a complete future null infinity. Moreover the support of the matter is confined to the region between two outgoing null hypersurfaces, and each of the Ricci coefficients, curvature components and components of the energy momentum tensor with respect to a double null frame decay towards null infinity with quantitative rates.
The proof of Theorem 1.1, after appealing to the corresponding result for the vacuum Einstein equations, quickly reduces to a semi-global problem. This reduction is outlined below and the semi-global problem treated here is stated in Theorem 1.2. Theorem 1.1 extends a result of Dafermos [12] which establishes the above under the additional restricted assumption of spherical symmetry. Note also the result of Rein-Rendall [29] which treats the massive case in spherical symmetry, where all of the particles have mass m > 0 (i.e. f is supported on the set of future pointing timelike vectors p in T M such that g( p, p) = −m 2 ). The main idea in [12] was to show, using a bootstrap argument, that, for sufficiently late times, the matter is supported away from the centre of spherical symmetry. By Birkhoff's Theorem the centre is therefore locally isometric to Minkowski space at these late times and the extension principle of Dafermos-Rendall [14] (see also [15]) then guarantees that the spacetime will be geodesically complete.
In these broad terms, a similar strategy is adopted here. The absence of good quantities satisfying monotonicity properties which are available in spherical symmetry, however, makes the process of controlling the support of the matter, and proving the semi-global existence result for the region where it is supported, considerably more involved. The use of Birkhoff's Theorem and the Dafermos-Rendall extension principle also have to be replaced by the much deeper result of the stability of Minkowski space for the vacuum Einstein equations. The use of the vacuum stability result, which is in fact appealed to in two separate places, is outlined below.

The Uncoupled Problem
It is useful to first recall what happens in the uncoupled problem of the Vlasov equation on a fixed Minkowski background. Let v = 1 2 (t +r ), u = 1 2 (t −r ) denote standard null Fig. 1 The projection of the support of f in the uncoupled problem coordinates on Minkowski space R 3+1 (these form a well defined coordinate system on the quotient manifold R 3+1 /SO (3) away from the centre of spherical symmetry {r = 0}) and suppose f is a solution of the Vlasov equation (2) with respect to this fixed background arising from initial data with compact support in space. From the geometry of null geodesics in Minkowski space it is clear that the projection of the support of f to the spacetime is related to the projection of the initial support of f as depicted in the Penrose diagram in Figure 1.
In particular, for sufficiently late advanced time v 0 the matter will be supported away from the centre {r = 0}, and there exists a point q ∈ R 3+1 /SO (3), lifting to a (round) 2-sphere S ⊂ R 3+1 , with r (q) > 0 such that where J − (S) denotes the causal past of S and π : P → M denotes the natural projection.

Initial Data and First Appeal to the Vacuum Result
Recall that initial data for the Einstein-Vlasov system (1)-(2) consists of a 3-manifold with a Riemannian metric g 0 , a symmetric (0, 2) tensor K and an initial particle density function f 0 satisfying the constraint equations, for j = 1, 2, 3, where div 0 , tr 0 , R 0 denote the divergence, trace and scalar curvature of g 0 respectively, and T 00 , T 0 j denote (what will become) the 00 and 0 j components of the energy momentum tensor. See [30] for a discussion of initial data for the Einstein-Vlasov system. The topology of will here be assumed to be that of R 3 . The issue of constructing solutions to the constraint equations (5) will not be treated here. A theorem of Choquet-Bruhat [4] guarantees that, given such an initial data set, a solution to (1)-(2) will exist locally in time.
The initial density function f 0 is assumed to have compact support. It will moreover be assumed that f 0 and a finite number of its derivatives will be small initially. The precise condition is given in Section 5. Note the assumption of compact support for f 0 is in both the spatial variable x, and in the momentum variable p. As will become evident, the compact support in space is used in a crucial way. The assumption of compact support in momentum is made for simplicity and can likely be weakened. 2 Let B ⊂ be a simply connected compact set such that π(supp( f | P )) ⊂ B, where P denotes the mass shell over . By the domain of dependence property of the Einstein-Vlasov system the development of the complement of B in , D + ( B), will solve the vacuum Einstein equations, The stability of Minkowski space theorem for the vacuum Einstein equations then guarantees the stability of this region. See Klainerman-Nicolò [21] where exactly this situation is treated. In particular, provided g 0 , K satisfy a smallness condition 3 in B (i.e. they are suitably close to the g 0 , K of Minkowski space), there exists a future complete, outgoing null hypersurface N in this region which can be foliated by a family of 2-spheres, {S u 0 ,v } parameterised by v, approaching the round 2-sphere as v → ∞. Moreover the Ricci coefficients and curvature components of the spacetime will decay to their corresponding Minkowski values and, by taking g 0 , K suitably small, certain weighted quantities involving them can be made arbitrarily small on N . It will be assumed that g 0 , K are sufficiently small so that the precise conditions stated in Theorem 5.1 are satisfied on N . A second appeal to a form of the stability of Minkowski space result in the vacuum (which can be shown to also follow from the Christodoulou-Klainerman Theorem [10] using upcoming work) will be made in Section 1.4 below.

Cauchy Stability
By Cauchy stability for the Einstein-Vlasov system (see Choquet-Bruhat [4] or Ringström [30]), Cauchy stability for the geodesic equations and the considerations of Section 1.1, provided the initial data on are taken sufficiently small, there exists a 2sphere S ⊂ M and an incoming null hypersurface N such that S ⊂ N , Area(S) > 0, π(supp( f )) ∩ S = ∅, and In other words, the existence of the point q in the Penrose diagram of Figure 1 is stable. It can moreover be assumed that the N above and N intersect in one of the 2-spheres of the foliation of N , where v 0 can be chosen arbitrarily large. The induced data on N can be taken to be arbitrarily small, provided they are sufficiently small on .

A Second Version of the Main Theorem and Second Appeal to the Vacuum Result
A more precise version of the main result can now be stated. A final version, Theorem 5.1, is stated in Section 5. This is depicted in Figure 2. Theorem 1.1 follows from Theorem 1.2 by the considerations of Section 1.2, Section 1.3, and by another application of the vacuum stability of Minkowski space result with the induced data on a hyperboloid contained between the null hypersurfaces {u = u f } and {u = u f − 1}. The problem of stability of Minkowski space for the vacuum Einstein equations (6) with hyperboloidal initial data was treated by Friedrich [17], though his result requires the initial data to be asymptotically simple. This is,  in general, inconsistent with the induced data arising from Theorem 1.2. 6 Whilst a proof of the hyperboloidal stability of Minkowski space problem with initial data compatible with Theorem 1.2 can most likely be distilled from the work [10], there is currently no precise statement to appeal to. In future work it will be shown how one can alternatively appeal directly to [10] by extending the induced scattering data at null infinity and solving backwards, in the style of [13].

Theorem 1.2 Given characteristic initial data for the massless Einstein-Vlasov sys
A precise formulation of Theorem 1.1, including an explicit statement of the norms used in the first appeal to the vacuum result in Section 1.2 and the Cauchy stability argument of Section 1.3, will not be made here. The assumptions made in Theorem 5.1, the final version of Theorem 1.2, will be given some justification at various places in the introduction however. The remainder of the paper will concern Theorem 1.2, and in the remainder of the introduction its proof will be outlined. The greatest new difficulty is in obtaining a priori control over derivatives of f . The approach taken involves introducing the induced Sasaki metric on the mass shell P and estimating certain Jacobi fields on P in terms of geometric quantities on the spacetime (M, g). This approach is outlined in Section 1.14 below.
Note that the analogue of Theorem 1.2 for the vacuum Einstein equations (6) follows from a recent result of Li-Zhu [25].

The Bootstrap Argument
The main step in the proof of Theorem 1.2 is in obtaining global a priori estimates for all of the relevant quantities. Once they have been established there is a standard procedure for obtaining global existence, which is outlined in Section 12. The remainder of the discussion is therefore focused on obtaining the estimates.
Moreover, using a bootstrap argument, it suffices to show that if the estimates already hold in a given bootstrap region of the form {u 0 ≤ u ≤ u } ∩ {v 0 ≤ v ≤ v }, depicted in Figure 3, then they can be recovered in this region with better constants independently of u , v . This is extremely useful given the strongly coupled nature of the equations.
The better constants in the bootstrap argument arise from either estimating the quantities by the initial data on {v = v 0 } and {u = u 0 } or by 1 v 0 , and using the smallness of the initial data and the largeness of v 0 . Recall that, in the setting of Theorem 1.1, both the largeness of v 0 and the smallness of the induced data on N = {u = u 0 }, N = {v = v 0 } arise by taking the asymptotically flat Cauchy data on to be suitably small.

The Double Null Gauge
The content of the Einstein equations is captured here through the structure equations and the null Bianchi equations associated to the double null foliation (u, v). The constant u and constant v hypersurfaces are outgoing and incoming null hypersurfaces respectively, and intersect in spacelike 2-spheres which are denoted S u,v . This choice of gauge is made due to its success in problems which require some form of the null condition 7 to be satisfied. 8 See, for example, [7,9,13,21,27].
The foliation defines a double null frame (see Section 2.1) in which one can decompose the Ricci coefficients, which satisfy so called null structure equations, the Weyl (or conformal) curvature tensor, whose null decomposed components satisfy the null Bianchi equations, and the energy momentum tensor (which, by the Einstein equations (4), is equal to the Ricci curvature tensor).
It is the null structure and Bianchi equations which will be used, together with the Vlasov equation (2), to estimate the solution. Following the notation of [13,27], the null decomposed Ricci coefficients will be schematically denoted . Two examples are the outgoing shearχ , which is a (0, 2) tensor on the spheres S u,v , and the renormalised outgoing expansion trχ − 2 r , which is a function on the spacetime, renormalised using the function r so that the corresponding quantity in Minkowski space will vanish.
The null decomposed components of the Weyl curvature tensor will be schematically denoted ψ and the null decomposed components of the energy momentum tensor will be schematically denoted T . This schematic notation, together with the p-index notation described in Section 1.8 below, will be used to convey structural properties of the equations which are heavily exploited later.

The Schematic Form of the Equations
The null structure equations for the Ricci coefficients , which are stated in Section 2.5, take the following schematic form, 7 See Section 1.8. 8 It is well known that the Einstein equations in the harmonic gauge do not satisfy the classical null condition of [20]. Despite this fact, it has been shown by Lindblad-Rodnianski [26] that one can still prove stability of Minkowski space for the vacuum Einstein equations in this gauge. One could therefore imagine adopting a similar strategy to approach the current problem in the harmonic gauge. Note the recent work of Fajman-Joudioux-Smulevici [16] on the development of a vector field method for relativistic transport equations, which could play an important role in such an approach.
Here / ∇ 3 and / ∇ 4 denote the projections of the covariant derivatives in the incoming and outgoing null directions respectively to the spheres S u,v . The 1 r terms appear in the equations for the outgoing and incoming expansions trχ − 2 r , trχ + 2 r , which are renormalised using the function r . Each satisfies exactly one of the two form of equations (7) and hence are further decomposed as (3) or (4) depending on whether they satisfy an equation in the / ∇ 3 or / ∇ 4 direction respectively. It should be noted that there are further null structure equations satisfied by the Ricci coefficients which take different forms to (7), some of which will make an appearance later.
The Weyl curvature components ψ can be further decomposed into Bianchi pairs, defined in Section 3.1, which are denoted (ψ, ψ ) (examples are (ψ, ψ ) = (α, β) or (β, (ρ, σ ))). This notation is used to emphasise a special structure in the Bianchi equations, which take the form, Here / D denote certain angular derivative operators on the spheres of intersection of the double null foliation, and ∇T schematically denote projected covariant derivatives of T in either the 3, 4 or angular directions.
The Ricci coefficients can be estimated using transport estimates for the null structure equations (7) since derivatives of do not appear explicitly on the right hand sides of the equations. The transport estimates are outlined below in Section 1.11 and carried out in detail in Section 10. Note that using such estimates does, however, come with a loss, namely the expected fact that angular derivatives of live at the same level of differentiability as curvature is not recovered. This fact can be recovered through a well known elliptic procedure, which is outlined below in Section 1.12 and treated in detail in Section 11. One cannot do the same for the curvature components and the Bianchi equations (8) due to the presence of the / Dψ terms on the right hand sides. In order to obtain "good" estimates for the Bianchi equations one must exploit the special structure which, if S denotes one of the spheres of intersection of the null foliation, takes the following form, i.e. the adjoint of the operator / D ψ is − / D ψ . Using this structure, if one contracts the / ∇ 3 ψ equation with ψ and adds the / ∇ 4 ψ equation contracted with ψ , the terms involving the angular derivatives will cancel upon integration and an integration by parts yields energy estimates for the Weyl curvature components. It is through this procedure that the hyperbolicity of the Einstein equations manifests itself in the double null gauge. These energy estimates form the content of Section 9 and, again, are outlined below in Section 1.10.
We are therefore forced (at least at the highest order) to estimate the curvature components in L 2 . All of the estimates for the Ricci coefficients here will also be L 2 based. In order to deal with the nonlinearities in the error terms of the equations, the same L 2 estimates are obtained for higher order derivatives of the quantities and Sobolev inequalities are used to obtain pointwise control over lower order terms. 9 To do this, a set of differential operators D is introduced which satisfy the commutation principle of [13]. This says that the "null condition" satisfied by the equations (which is outlined below and crucial for the estimates) and the structure discussed above are preserved when the equations are commuted by D, i.e. D and Dψ satisfy similar equations to and ψ. The set of operators D is introduced in Section 3.3.
As they appear on the right hand side of the equations for ψ, , the energy momentum tensor components T are also, at the highest order, estimated in L 2 . These estimates are obtained by first estimating f using the Vlasov equation. It is important that the components of the energy momentum tensor, and hence also f , are estimated at one degree of differentiability greater than the Weyl curvature components ψ. The main difficulty in this work is in obtaining such estimates for the derivatives of f . See Section 1.14 for an outline of the argument and Section 8 for the details.

The p-Index Notation and the Null Condition
The discussion in the previous section outlines how one can hope to close the estimates for and ψ from the point of view of regularity. Since global estimates are required, it is also crucial that all of the error terms in the equations decay sufficiently fast in v (or equivalently, since everything takes place in the "wave zone" where r := v − u + r 0 is comparable to v, sufficiently fast in r ) so that, when they appear in the estimates, they are globally integrable. For quasilinear wave equations there is an algebraic condition on the nonlinearity, known as the null condition, which guarantees this [20]. By analogy, we say the null structure and Bianchi equations "satisfy the null condition" to mean that, on the right hand sides of the equations, certain "bad" combinations of the terms do not appear. There is an excellent discussion of this in the introduction of [13]. As they are highly relevant, the main points are recalled here.
Following [13], the correct hierarchy of asymptotics in r for , ψ and T is first guessed. This guess is encoded in the p-index notation. Each , ψ, T is labelled with a subscript p to reflect the fact that r p | p |, r p |ψ p |, r p |T p | are expected to be uniformly bounded. 10 Here | · | denotes the norm with respect to the induced metric on the 2spheres / g. The weighted L 2 quantities which will be shown to be uniformly bounded will imply, via Sobolev inequalities, that this will be the case at lower orders.
In Theorem 5.1, the precise formulation of Theorem 1.2, it is asymptotics consistent with the p-index notation which will be assumed to hold on the initial outgoing hypersurface N = {u = u 0 }. In the context of Theorem 1.1, recall the use of the Klainerman-Nicolò [21] result in Section 1.2. The result of Klainerman-Nicolò guarantees that, provided the asymptotically flat Cauchy data on has sufficient decay, there indeed exists an outgoing null hypersurface in the development of the data on which asymptotics consistent with the p-index notation hold.

Geometry of Null Geodesics and the Support of f
If the Ricci coefficients are assumed to have the asymptotics described in the previous section then it is straightforward to show that u f can be chosen to have the desired property that f = 0 on the mass shell over any point x ∈ M with u(x) ≥ u f − 1. In fact, it can also be seen that the size of the support of f in P x , the mass shell over the point x ∈ M, will decay as v(x) → ∞. This decay is important as it is used to obtain the decay of the components of the energy momentum tensor. The argument for obtaining the decay properties of supp( f ) is outlined here and presented in detail in Section 7.
The decay of the size of the support of f in P x can be seen by considering the decay of components of certain null geodesics. Suppose first that γ is a future directed null geodesic in Minkowski space emanating from a compact set in the hypersurface {t = 0} such that the initial tangent vectorγ (0) is contained in a compact set in the mass shell over {t = 0}. One can show that, iḟ where e 1 = ∂ θ 1 , e 2 = ∂ θ 2 , e 3 = ∂ u , e 4 = ∂ v is the standard double null frame in Minkowski space, then the bounds, hold uniformly along γ for some constant C. 11 The bounds (9) will be assumed to hold in supp( f ) in the mass shell over the initial hypersurface {v = v 0 } in Theorem 5.1, the precise formulation of Theorem 1.2. In the setting of Theorem 1.1, the bounds (9) can be taken to hold on the hypersurface N = {v = v 0 } in view of the Cauchy stability argument of Section 1.3 and the fact 11 One sees these are the correct asymptotics for p 1 (s), p 2 (s) by using the three angular momentum Killing vector fields of Minkowski space 1 , 2 , 3 and the fact that, if K is a Killing vector, g(γ , K ) is constant along γ , to see that the angular momentum of the geodesic, is conserved along γ . This fact together with the mass shell relation 4 p 3 p 4 = / g AB p A p B and the fact that p 4 does not decay along a null geodesic (which can be seen in Minkowski space by looking at the geodesic equations and noting thatṗ 4 (s) ≥ 0) gives the required asymptotics for p 3 . that they hold globally in supp( f ) for the uncoupled problem of the Vlasov equation on a fixed Minkowski background.
The idea is now to propagate the bounds (9) from the initial hypersurface {v = v 0 } into the rest of the spacetime. If e 1 , . . . , e 4 now denotes the double null frame of (M, g) (defined in Section 2.1), one then uses the geodesic equations, for a null geodesic γ withγ (s) = p μ (s)e μ | γ (s) , a bootstrap argument and the pointwise bounds r p | p | ≤ C to see thaṫ The estimates (9) follow by integrating along γ since dr ds ∼ p 4 (0). Finally, to show the retarded time u f can be chosen as desired, let u(s) denote the u coordinate of the geodesic γ at time s. Then and hence |u(s) − u(0)| ≤ C for all s ∈ [0, ∞), for some constant C.

Global Energy Estimates for the Curvature Components
The global energy estimates for the Weyl curvature components can now be outlined. They are carried out in detail in Section 9. The Bianchi equations take the schematic form, where c is a constant (which is different for the different ψ p ) and E p is an error which will decay, according to the p notation, like 1 r p . Similarly, E p + 3 2 is an error which will decay like 1 Recall from equation (8) that the errors E p and E p + 3 2 contain linear terms involving , nonlinear terms of the form · ψ and · T , and projected covariant derivatives of components of the energy momentum tensor ∇T . Using (10) to compute Div r w |ψ p | 2 e 3 , Div r w |ψ p | 2 e 4 , after summing a cancellation will occur in the terms involving angular derivatives, as discussed in Section 1.7, and they can be rewritten as a spherical divergence. If the weight w is chosen correctly, a cancellation 12 also occurs in the ctrχ ψ p term (which, since trχ looks like 2 r to leading order, cannot be included in the error E p + 3 9 Page 14 of 177 where Div denotes the spacetime divergence and B denotes a spacetime "bulk" region bounded to the past by the initial characteristic hypersurfaces, and to the future by constant v and constant u hypersurfaces. See Figure 3. Note that this procedure will generate additional error terms but they can be treated similarly to those arising from the errors in (10) and hence are omitted here. See Section 9 for the details.
If the curvature fluxes are defined as, then by the divergence theorem, when the above identity (11) is summed over all Bianchi pairs (ψ p , ψ p ), the left hand side becomes Due to the relation between the weights w(ψ p , ψ p ) and p, p , and the bounds assumed for and T through the bootstrap argument, the right hand side of (11) can be controlled by, for some constant C (which, of course, arises from inserting the bootstrap assumptions). It is this step where one sees the manifestation of the null condition in the Bianchi equations. Dropping the F 2 u 0 ,u (v) term on the left yields, and hence, by the Grönwall inequality, F 1 v 0 ,v (u) can be controlled by initial data and the term C v 0 . Returning to the inequality, and inserting the above bounds for F 1 v 0 ,v (u), F 2 u 0 ,u (v) can now also be similarly controlled.

Global Transport Estimates for the Ricci Coefficients
Turning now to the global estimates for the Ricci coefficients, which are treated in detail in Section 10, in the p-index notation the null structure equations take the form, where again E p is an error which decays, according to the p-index notation, like 1 r p and E p+2 decays like 1 r p+2 . Recall from equation (7) that E p and E p+2 contain linear terms involving , ψ, T , and quadratic terms of the form · . The / ∇ 4 p equations can be rewritten as To estimate the (4) one then uses the identity, for a function h on M, where the trχ h term comes from the derivative of the volume form on S u,v , with The r −2 factor serves to cancel the trχ term (which, recall, behaves like 2 r and so is not globally integrable in v). Hence, since the volume form is of order r 2 . Integrating in v from the initial hypersurface {v = v 0 } then gives, Note that the error E p+2 is integrated over a u = constant hypersurface. These are exactly the regions on which the integrals of the Weyl curvature components were controlled in Section 1.10, and it is for this reason the curvature terms in the error E p+2 can be controlled in (12).
Since the volume form is of order r 2 , the bound (12) is consistent with (4) p decaying like 1 r p and, after repeating the above with appropriate derivatives of (4) p , this pointwise decay can be obtained using Sobolev inequalities on the spheres.
It is not a coincidence that the p 2 coefficient of trχ (4) p in the / ∇ 4

(4)
p equation is exactly that which is required to obtain 1 r p decay for (4) p . In fact some of the (4) p will decay faster than this but the other null structure equations are required, along with elliptic estimates, to obtain this. It is therefore the p 2 coefficient which determines the p index given here to the (4) as it restricts the decay which can be shown to hold using only the / ∇ 4 (4) equations. Note the difference with [13] where the authors are not constrained by this coefficient as they there integrate "backwards" from future null infinity.
Turning now to the equations in the 3 direction, the (3) p quantities are estimated using the identity, It does not now matter that trχ only decays like 1 r since the integration in u will only be up to the finite value u f .

Suppose first that
where E p+1 decays like 1 r p+1 and E 0 p decays like 1 r p but only contains Weyl curvature, energy momentum tensor and (4) terms which have already been estimated (the energy momentum tensor estimates are outlined below as they present the greatest difficulty but in the logic of the proof are estimated first). Then, Integrating from u 0 to u and inserting the bootstrap assumptions and the previously obtained bounds for E 0 p , the Grönwall inequality then gives, where ε 0 controls the size of the initial data. Note that it was important that the only error terms which have not already been estimated are of the form E p+1 , and not E p , in order to gain the 1 v 0 smallness factor. It turns out that there is a reductive structure in the null structure equations so that, provided they are estimated in the correct order, each (3) satisfies an equation of the form (13) where E 0 p now also contains (3) terms which have been estimated previously. Hence all of the (3) can be estimated with smallness factors.

Elliptic Estimates and Ricci Coefficients at the Top Order
The procedure in Section 1.11 is used to estimate the Ricci coefficients, along with their derivatives at all but the top order, in L 2 of the spheres of intersection of constant u and constant v hypersurfaces. The derivatives of Ricci coefficients at the top order are estimated only in L 2 on null hypersurfaces. These estimates are obtained using elliptic equations on the spheres for some of the Ricci coefficients, coupled to transport equations for certain auxilliary quantities. This procedure is familiar from many other works (e.g. [9,10]) and forms the content of Section 11. It should be noted that these estimates are only required here for estimating the components of the energy momentum tensor. If one were to restrict the semi-global problem of Theorem 1.2 to the case of the vacuum Einstein equations (6) then the estimates for the Ricci coefficients and curvature components could be closed with a loss (i.e. without knowing that angular derivatives of Ricci coefficients lie at the same degree of differentiability as the Weyl curvature components) as only the null structure equations of the form (7) would be used, and these elliptic estimates would not be required. See Section 1.7.

Global Estimates for the Energy Momentum Tensor Components
At the zeroth order the estimates for the energy momentum tensor components follow directly from the bounds (9), which show that the size of the region supp( f | P x ) ⊂ P x on which the integral in (1) is taken is decaying as r (x) → ∞, and the fact that f is conserved along trajectories of the geodesic flow. For example, using the volume form for P x defined in Section 2.2, if sup {v=v 0 } | f | ≤ ε 0 , In fact, provided the derivatives of f can be estimated, the estimates for the derivatives of T are obtained in exactly the same way.

Global Estimates for Derivatives of f
A fundamental new aspect of this work arises in obtaining estimates for the derivatives of f . Recall from Section 1.7 that, in order to close the bootstrap argument, it is crucial that the energy momentum tensor components T , and hence f , are estimated at one degree of differentiability greater than the Weyl curvature components, i.e. k derivatives of f must be estimated using only k − 1 derivatives of ψ. Written in components with respect to the frame 13 e 1 , e 2 , e 3 , e 4 , ∂ p 1 , ∂ p 2 , ∂ p 4 for P, the Vlasov equation (2) takes the form, where μ νλ denote the Ricci coefficients of M. See (24)-(28) below. One way to estimate derivatives of f is to commute this equation with suitable vector fields and integrate along trajectories of the geodesic flow. If V denotes such a vector field, commuting will give, where E is an error involving terms of the form V ( μ νλ ). At first glance this seems promising as derivatives of the Ricci coefficients should live at the same level of differentiability as the Weyl curvature components ψ. This is not the case for all of the and hence, commuting twice with angular derivatives and using the elliptic estimates described in Section 1.12 will only give estimates for two angular derivatives of b by first order derivatives of ψ and T . The angular derivatives of the spherical Christoffel symbols / , see (22) below, which also appear when commuting the Vlasov equation give rise to similar issues.
Whilst it may still be the case that E as a whole (rather than each of its individual terms) can be estimated just at the level of ψ, a different approach is taken here in order to see more directly that derivatives of f can be estimated at the level of ψ. This approach, which is treated in detail in Section 8, is outlined now.
Consider again a vector V ∈ T (x, p) P. Recall the form of the Vlasov equation (3). Using this expression for f and the chain rule, for any s, and hence, if J (s) 13 Recall e 1 , e 2 , e 3 , e 4 is the double null frame for M, defined in Section 2.1 using the (u, v, θ 1 , θ 2 ) coordinate system on M. Note the slight abuse of notation here as e 1 , e 2 , e 3 , e 4 act only on functions on M, whilst f is a function on P. As e 1 = ∂ θ 1 in the (u, v, θ 1 , θ 2 ) coordinate system, e 1 f is used to denote ∂ θ 1 f in the (u, v, θ 1 , θ 2 , p 1 , p 2 , p 4 ) coordinate system for P. Similarly for e 2 f, e 3 f, e 4 f .
If s < 0 is taken so that π(exp s (x, p)) ∈ {v = v 0 } then the expression (14) relates a derivative of f at (x, p) to a derivative of f on the initial hypersurface. It therefore remains to estimate the components of J (s), with respect to a suitable frame for P, uniformly in s and independently of the point (x, p). The metric g on the spacetime M can be used to define a metric on the tangent bundle T M, known as the Sasaki metric [33], which by restriction defines a metricĝ on the mass shell P. See Section 4 where this metric is introduced. With respect to this metric trajectories of the geodesic flow s → exp s (x, p) are geodesics in P and, for any vector V ∈ T (x, p) P, J (s) := d exp s | (x, p) V is a Jacobi field along this geodesic (see Section 4). Therefore J (s) satisfies the Jacobi equation, where∇ denotes the induced connection on P, andR denotes the curvature tensor of (P,ĝ). Equation (15) is used, as a transport equation along the trajectories of the geodesic flow, to estimate the components of J . The curvature tensorR can be expressed in terms of (vertical and horizontal lifts of) the curvature tensor R of (M, g) along with its first order covariant derivatives ∇ R. See equation (90). At first glance the presence of ∇ R again appears to be bad. On closer inspection, however, the terms involving covariant derivatives of R are always derivatives in the "correct" direction so that they can be recovered by the transport estimates, and the components of J , and hence V f , can be estimated at the level of ψ.
The above observations of course only explain how one can hope to close the estimates for T from the point of view of regularity. In order to obtain global estimates for the components of J one has to use the crucial fact that, according to the p-index notation, the right hand side of the Jacobi equationR(X, J )X , when written in terms of ψ, T , p 1 , p 2 , p 3 , p 4 , decays sufficiently fast as to be twice globally integrable along s → exp s (x, p). This can be viewed as a null condition for the Jacobi equation and is brought to light through further schematic notation introduced in Section 8.2.
The fact that the right hand side of (15) has sufficient decay in r is perhaps not surprising. Consider for example the term inR(X, J )X . Here γ is a geodesic in M such that exp s (x, p) = (γ (s),γ (s)) and J h is a vector field along γ on M such that, together with another vector field J v along γ , with Hor (γ ,γ ) and Ver (γ ,γ ) denoting horizontal and vertical lifts at (γ ,γ ) (defined in Section 4). The slowest decaying ψ and T are those which contain the most e 3 vectors. Whenever such ψ and T arise in (16) however, they will typically be accompanied by p 3 (s), the e 3 component ofγ (s), which (recall from Section 1.9) has fast 1 r (s) 2 decay. Similarly the non-decaying p 4 (s), the e 4 component ofγ (s), can only appear in (16) accompanied by the ψ and T which contain e 4 vectors and hence have fast decay in r . In particular, potentially slowly decaying terms involving p 4 (s) multiplying the ψ and T which contain no e 4 vectors do not arise in (16).
Finally, since J f now itself is also conserved along s → exp s (x, p), second order derivatives of f can be obtained by repeating the above. If J 1 , J 2 denote Jacobi fields corresponding to vectors V 1 , V 2 at (x, p) respectively, then, In order to control V 2 V 1 f (x, p) it is therefore necessary to estimate the J 2 derivatives of the components of J 1 along s → exp s (x, p). This is done by commuting the Jacobi equation (15) and showing that the important structure described above is preserved. The Jacobi fields which are used, and hence the vectors V used to take derivatives of f , have to be carefully chosen. They are defined in Section 8.3.
Note that this procedure can be repeated to obtain higher order derivatives of f . Whilst the pointwise bounds on ψ at lower orders mean that lower order derivatives of f can be estimated pointwise, at higher orders this procedure will generate terms involving higher order derivatives of ψ and hence higher order derivatives of T must be estimated in L 2 on null hypersurfaces. In fact, at the very top order, T is estimated in the spacetime L 2 norm.

Related Previous Stability Results in General Relativity
There are several related previous works on the stability of Minkowski space for the Einstein equations coupled to various matter models. Without simplifying symmetry assumptions, the first such work was that of Christodoulou-Klainerman [10]. They show that, given an initial data set for the vacuum Einstein equations, satisfying an appropriate asymptotic flatness condition, which is sufficiently close to that of Minkowski space, the resulting maximal development is geodesically complete, possesses a complete future null infinity, and asymptotically approaches Minkowski space with quantitative rates. The result, more fundamentally, provided the first examples of smooth, geodesically complete, asymptotically flat solutions to the vacuum Einstein equations, other than Minkowski space itself. The existence of such spacetimes if far from trivial. The proof relies on foliating the spacetimes they construct by the level sets of a so called maximal time function, along with another, null, foliation by the level sets an optical function. Detailed behaviour of the solutions are obtained, along with various applications including a rigorous derivation of the Bondi mass loss formula.
The proof of Christodoulou-Klainerman was generalised by Zipser [36], who showed that the analogue of their theorem holds for electromagnetic matter described by the Maxwell equations. The proof of this generalisation again relies on foliating the spacetimes by the level hypersurfaces of a maximal time and optical function.
The Christodoulou-Klainerman proof was later revisited by Klainerman-Nicolò [21] who showed the stability of the domain of dependence of the complement of a ball in a standard spacelike hypersurface in Minkowski space. Their smallness condition on initial data is similar to that of Christodoulou-Klainerman, however the Klainerman-Nicolò proof is based on a double null foliation, defined as the level hypersurfaces of two, outgoing and incoming, optical functions. The proof of Theorem 1.2 in this work is based on a similar approach and moreover, since the Klainerman-Nicolò result is appealed to in its proof, the smallness condition required in Theorem 1.1 is similar to that of [21]. See Section 1.2.
A new proof of the stability of Minkowski space for the vacuum Einstein equations using the harmonic gauge, the gauge originally used by Choquet-Bruhat [3] to prove local existence for the vacuum Einstein equations, was developed by Lindblad-Rodnianski [26]. Their proof essentially reduces to a small data global existence proof for a system of quasilinear wave equations and, despite the equations failing to satisfy the classical null condition of Klainerman [20], is relatively technically simple. The proof moreover requires a weaker asymptotic flatness condition on the data, compared to [10], and also allows for coupling to matter described by a massless scalar field. The asymptotic behaviour obtained for the solutions is less precise, however, than in [10].
The Christodoulou-Klainerman proof was returned to again by Bieri [2], who imposes a weaker asymptotic flatness condition on the initial data, in terms of decay, and is able to close the proof using fewer derivatives of the solution than [10]. The proof again, as in [10], is based on a maximal-null foliation of the spacetimes.
The stability of Minkowski space problem for the Einstein-Maxwell system, as studied by Zipser, was returned to by Loizelet [24], this time using the harmonic gauge approach of Lindblad-Rodnianski. The harmonic gauge approach was also used by Speck [34], who considers the Einstein equations coupled to a large class of electromagnetic equations, which are derivable from a Lagrangian and reduce to the Maxwell equations in an appropriate limit. A recent result of LeFloch-Ma [23] on the problem for the Einstein-Klein-Gordon system also uses the harmonic gauge approach (see also [35]).
Finally, there are more global stability results for the Einstein equations with a positive cosmological constant, for example the works of Friedrich [17], Ringström [30] and Rodnianski-Speck [31]. A more comprehensive list can be found in the introduction to the work of Hadzić-Speck [19].

Outline of the Paper
In the next section coordinates are defined on the spacetime to be constructed, and on the mass shell P. The Ricci coefficients and curvature components are introduced along with their governing equations. In Section 3 the schematic form of the quantities and equations are given. Three derivative operators are then introduced which are shown to preserve the schematic form of the equations under commutation. Some facts about the Sasaki metric are recalled in Section 4 and are used to describe certain geometric properties of the mass shell. A precise statement of Theorem 1.2 is given in Section 5, along with the statement of a bootstrap theorem. The proof of the bootstrap theorem is given in the following sections. The main estimates are obtained for the energy momentum tensor components, Weyl curvature components and lower order derivatives of Ricci coefficients in Sections 8, 9 and 10 respectively. The estimates for the Ricci coefficients at the top order are obtained in Section 11. The results of these sections rely on the Sobolev inequalities of Section 6, and the decay estimates for the size of supp( f | P x ) ⊂ P x as x approaches null infinity from Section 7. The fact that the retarded time u f can be chosen to have the desired property, stated in Theorem 1.2, is also established in Section 7. Finally, the completion of the proof of Theorem 1.2, through a last slice argument, is outlined in Section 12.

Basic Setup
Throughout this section consider a smooth spacetime is a manifold with corners and g is a smooth Lorentzian metric on M such that (M, g), together with a continuous function f : P → [0, ∞), smooth on P Z , where Z denotes the zero section, satisfy the Einstein-Vlasov system (1)-(2).

Coordinates and Frames
A point in M will be denoted (u, v, θ 1 , θ 2 ). It is implicitly understood that two coordinate charts are required on S 2 . The charts will be defined below using two coordinate charts on S u 0 ,v 0 = {u = u 0 } ∩ {v = v 0 }. Assume u and v satisfy the Eikonal equation Following [9,21], define null vector fields Then extend to u > u 0 by solving This defines coordinates (u, v, θ 1 , θ 2 ) on the region D(U 1 ) defined to be the image of U 1 under the diffeomorphisms generated by L on {u = u 0 }, then by the diffeomorphisms generated by L. Coordinates can be defined on another open subset of the spacetime by considering coordinates in another region U 2 ⊂ S u 0 ,v 0 and repeating the procedure. These two coordinate charts will cover the entire region of the spacetime in question provided the charts U 1 , U 2 cover S u 0 ,v 0 . The choice of coordinates on U 1 , U 2 is otherwise arbitrary. The spheres of constant u and v will be denoted S u,v and the restriction of g to these spheres will be denoted / g. A vector field V on M will be called an Similarly for (r, 0) tensors. A one form ξ is called an S u,v one form if ξ(L) = ξ(L) = 0. Similarly for (0, s), and for general (r, s) tensors.
In these coordinates the metric takes the form where b is a vector field tangent to the spheres S u,v , which vanishes on the initial hypersurface {u = u 0 }. Note that, due to the remaining gauge freedom, can be specified on {u = u 0 } and {v = v 0 }. Since, in Theorem 5.1, it is assumed that 1 2 − 1 and r 3 | / ∇ 4 log | are small on {u = u 0 }, it is convenient to set = 1 on {u = u 0 } so that they both vanish.
Integration of a function φ on S u,v is defined as where τ 1 , τ 2 is a partition of unity subordinate to D U 1 , D U 2 at u, v.
Define the double null frame and let ( p μ ; μ = 1, 2, 3, 4), denote coordinates on each tangent space to M conjugate to this frame, so that the coordinates (x μ , p μ ) denote the point . This then gives a frame, {e μ , ∂ p μ | μ = 1, 2, 3, 4}, on T M. The Vlasov equation (2) written with respect to this frame takes the form where μ νλ are the Ricci coefficients of g with respect to the null frame (18). For f as a function on the mass shell P, this reduces to, whereμ now runs over 1, 2, 4, and p 1 , p 2 , p 4 denote the restriction of the coordinates p 1 , p 2 , p 4 to P, and ∂ pμ denote the partial derivatives with respect to this restricted coordinate system. Using the mass shell relation (21) below one can easily check, Note that Greek indices, μ, ν, λ, etc. will always be used to sum over the values 1, 2, 3, 4, whilst capital Latin indices, A, B, C, etc. will be used to denote sums over only the spherical directions 1, 2. In Section 8 lower case latin indices i, j, k, etc. will be used to denote summations over the values 1, . . . , 7.
Remark 2.1 A seemingly more natural null frame to use on M would be Dafermos-Holzegel-Rodnianski [13] use the same "unnatural" frame for regularity issues on the event horizon. The reason for the choice here is slightly different and is related to the fact that ω, defined below, is zero in this frame.

Null Geodesics and the Mass Shell
Recall that the mass shell P ⊂ T M is defined to be the set of future pointing null vectors. Using the definition of the coordinates p μ and the form of the metric given in the previous section one sees that, since all of the particles have zero mass, i.e. since f is supported on P, the relation is true in the support of f . The identity (21) is known as the mass shell relation. The mass shell P is a 7 dimensional hypersurface in T M and can be parameterised by coordinates (u, v, θ 1 , θ 2 , p 1 , p 2 , p 4 ), with p 3 defined by (21).
To make sense of the integral in the definition of the energy momentum tensor (1) one needs to define a suitable volume form on the mass shell, P x , over each point Since P x is a null hypersurface it is not immediately clear how to do this. Given such an x, the metric on M defines a metric on T x M, −4dp 3 dp 4 + / g AB (x)dp A dp B , which in turn defines a volume form on T x M, 2 det / gdp 3 ∧ dp 4 ∧ dp 1 ∧ dp 2 .
A canonical one-form normal to P x can be defined as the differential of the function x (X ) := g(X, X ).
See Section 5.6 of [32]. The energy momentum tensor at x ∈ M therefore takes the form

Ricci Coefficients and Curvature Components
Following the notation of [9] (see also [10,21]), define the Ricci coefficients The null second fundamental forms χ, χ are decomposed into their trace and trace free parts Note that due to the choice of frame, since e 3 is an affine geodesic vector field, ω := 1 2 g(∇ e 3 e 4 , e 3 ) = 0. Also note that in this frame ζ A := 1 2 g(∇ e A e 4 , e 3 ) = −η A . The Christoffel symbols of (S u,v , / g) with respect to the frame e 1 , e 2 are denoted / C AB , Define also the null Weyl curvature components Here 14 is the Weyl, or conformal, curvature tensor of (M, g) and * W denotes the hodge dual of W , where is the spacetime volume form of (M, g).
Define the S u,v (0,2)-tensor / T to be the restriction of the energy momentum tensor defined in equation (1) to vector fields tangent to the spheres S u,v : Similarly let / T 3 , / T 4 denote the S u,v 1-forms defined by restricting the 1-forms T (e 3 , ·), T (e 4 , ·) to vector fields tangent to the spheres S u,v :

The Minkowski Values
For the purpose of renormalising the null structure and Bianchi equations, define the following Minkowski values of the metric quantities using the function r := v −u +r 0 , with r 0 > 0 a constant chosen to make sure r ≥ inf u Area(S u,v 0 ), 14 Recall that the scalar curvature R vanishes for solutions of the massless Einstein-Vlasov system.
where γ is the round metric on the unit sphere. Similarly, define and let / •C AB denote the spherical Christoffel symbols of the metric / g • = r 2 γ , so that, where / ∇ • is the Levi-Civita connection of / g • . These are the only non-identically vanishing Ricci coefficients in Minkowski space. All curvature components vanish, as do all components of the energy momentum tensor.
Note that the function r in general does not have the geometric interpretation as the area radius of the spheres S u,v . Note also that

The Renormalised Null Structure and Bianchi Equations
The Bianchi equations, written out in full using the table of Ricci coefficients, take the form 15 15 See [9] for a detailed derivation in the vacuum case. Recall that ζ = −η and ω = 0 in the frame used here, and that the scalar curvature R vanishes for solutions of the massless Einstein-Vlasov system.
Here for an S u,v 1-form ξ ,/ ∇ξ denotes the transpose of the derivative of ξ , The left Hodge-dual * is defined on S u,v one forms and (0, 2) S u,v tensors by respectively. Here / denotes the volume form associated with the metric / g and, for a (0, 2) S u,v tensor ξ , The symmetric traceless product of two S u,v one forms is defined by and the anti-symmetric products are defined by for two S u,v one forms and S u,v (0, 2) tensors respectively. Also, and the norm of a (0, n) S u,v tensor The notation | · | will also later be used when applied to components of S u,v tensors to denote the standard absolute value on R. See Section 6. It will always be clear from the context which is meant, for example if ξ is an S u,v 1-form then |ξ | denotes the / g norm as above, whilst |ξ A | denotes the absolute value of ξ(e A ).
The null structure equations for the Ricci coefficients and the metric quantities in the 3 direction, suitably renormalised using the Minkowski values, take the form and in the 4 direction Through most of the text, when referring to the null structure equations it is the above equations which are meant. The following null structure equations on the spheres will also be used in Section 11, where K denotes the Gauss curvature of the spheres (S u,v , / g). The additional propagation equations forχ,χ , will also be used in Section 11 to derive propagation equations for the mass aspect function μ, μ defined later. Here/ T = / T − / T 34 / g is the trace free part of / T . The following first variational formulas for the induced metric on the spheres will also be used, where L denotes the Lie derivative. There are additional null structure equations but, since they will not be used here, are omitted.

The Schematic Form of the Equations and Commutation
In this section schematic notation is introduced for the Ricci coefficients, curvature components and components of the energy momentum tensor, which is used to isolate the structure in the equations that is important for the proof of Theorem 1.2. A collection of differential operators is introduced and it is shown that this structure remains present after commuting the equations by any of the operators in the collection. This section closely follows Section 3 of [13] where this notation was introduced.

Schematic Notation
Consider the collection of Ricci coefficients 16 which are schematically denoted , 16 The quantities 1 2 − 1, b and / g − / g • are, of course, metric quantities. They are however, in Section 10, estimated systematically along with the Ricci coefficients. Any discussion of the "Ricci coefficients" from now on will hence implicitly refer also to these metric quantities.
Note that the are normalised so that each of the corresponding quantities in Minkowski space is equal to zero. In the proof of the main result it will be shown that each converges to zero as r → ∞ in the spacetimes considered. Each will converge with a different rate in r and so, to describe these rates, each is given an index, p, to encode the fact that, as will be shown in the proof of the main result, r p | p | will be uniformly bounded. The p-indices are given as follows, so that 1 schematically denotes any of the quantitiesχ , / g − / g • , b, η, etc. It may be the case, for a particular p , that lim r →∞ r p | p | is always zero in each of the spacetimes which are constructed here. This means that some of the Ricci coefficients will decay with faster rates than those propagated in the proof of Theorem 1.2. Some of these faster rates can be recovered a posteriori.
The notation (3) will be used to schematically denote any for which the corresponding null structure equation of (39)-(48) it satisfies is in the / ∇ 3 direction, Similarly, (4) will schematically denote any for which the corresponding null structure equation of (39)-(48) is in the / ∇ 4 direction, p will schematically denote any p which has also been denoted (3) . So, for example,χ may be schematically denoted p will schematically denote any p which has also been denoted (4) . Consider now the collection of Weyl curvature components, which are schematically denoted ψ, Each ψ is similarly given a p-index, to encode the fact that, as again will be shown, r p |ψ p | is uniformly bounded in each of the spacetimes which are constructed. When deriving energy estimates for the Bianchi equations in Section 9, a special divergence structure present in the terms involving angular derivatives is exploited. For example, the / ∇ 3 α equation is contracted with α (multiplied by a suitable weight) and integrated by parts over spacetime. The / ∇ 4 β equation is similarly contracted with β and integrated by parts. When the two resulting identities are summed, a cancellation occurs in the terms involving angular derivatives leaving only a spherical divergence which vanishes due to the integration on the spheres. The / ∇ 3 α equation is thus paired with the / ∇ 4 β equation. To highlight this structure, consider the ordered pairs, Each of these ordered pairs will be schematically denoted (ψ p , ψ p ), with the subscripts p and p as in (56), and referred to as a Bianchi pair.
The components of the energy momentum tensor are schematically denoted T , and each T is similarly given a p-index, to encode the fact that r p |T p | will be shown to be uniformly bounded. Finally, for a given p ∈ R, let h p denote any smooth function h p : M → R, depending only on r , which behaves like 1 r p to infinite order, i.e. any function such that, for any k ∈ N 0 , there is a constant C k such that r k+ p |(∂ v ) k h p | ≤ C k , where the derivative is taken in the (u, v, θ 1 , θ 2 ) coordinate system. In addition, the tensor field

The Schematic Form of the Equations
Using the notation of the previous section, the null structure and Bianchi equations can be rewritten in schematic form. For example the null structure equation (40) can be rewritten, Here and in the following, p 1 · p 2 denotes (a constant multiple of) an arbitrary contraction between a p 1 and a p 2 . In the estimates later, the Cauchy-Schwarz inequality | p 1 · p 2 | ≤ C| p 1 || p 2 | will always be used and so the precise form of the contraction will be irrelevant. Similarly for h p 1 p 2 .
Rewriting the equations in this way allows one to immediately read off the rate of decay in r of the right hand side. In the above example one sees that / ∇ 3 2 is equal to a combination of terms whose overall decay is, according to the p-index notation, like 1 r 2 , consistent with the fact that applying / ∇ 3 to a Ricci coefficient does not alter its r decay (see Section 3.3). Each of the null structure equations can be expressed in this way.

Proposition 3.1 (cf. Proposition 3.1 of [13]). The null structure equations (39)-(48)
can be written in the following schematic form where 17 This proposition allows us to see that the right hand sides of the / ∇ 3 p equations behave like 1 r p , whilst the right hand sides of the / ∇ 4

(4)
p equations behave like 1 r p+2 . This structure will be heavily exploited and should be seen as a manifestation of the null condition present in the Einstein equations.  (58) is not contained in the error since trχ behaves like 1 r and so this term only behaves like 1 r p+1 . This would thus destroy the structure of the error. It is not a problem that this term appears however since, when doing the estimates, the following renormalised form of the equation will always be used: This can be derived by differentiating the left hand side using the product rule, substituting equation (58) and using the fact that (trχ • − trχ) (4) p can be absorbed into the error.
It is not a coincidence that the coefficient of this term is always p 2 , it is the value of this coefficient which decides the rate of decay to be propagated for each (4) p . This will be elaborated on further in Section 10. This was not the case in [13]; they have more 17 Stricly speaking, the terms listed in the error can be "worse" than the terms which actually appear in the p ], ψ p may actually refer to some ψ q where q > p. In fact, in some of the null structure equations no curvature term actually appears. freedom since they are integrating backwards from null infinity, rather than towards null infinity, and so can propagate stronger decay rates for some of the (4) p . These stronger rates could be recovered here using the ideas in Section 11, however it is perhaps interesting to note that the estimates can be closed with these weaker rates.
The Bianchi equations can also be rewritten in this way.
where / D denotes the angular operator appearing in equation for the particular curvature component under consideration 18

The error terms take the form
where D is used to denote certain derivative operators which are introduced in Section 3.3.
When applied to T p , the operators D should not alter the rate of decay so again this schematic form allows one to easily read off the r decay rates of the errors. This structure of the errors will again be heavily exploited. The first summation in E 4 [ψ p ] can in fact actually always begin at p + 2 except for in E 4 [β] where the term η # · α appears. Also the terms, can be upgraded to, . These points are important and will be returned to in Section 9. 18 So, for example, / Dψ p = / ∇⊗β in the / ∇ 3 α equation, and / Dψ p = − / ∇ρ + * / ∇σ in the / ∇ 4 β equation, etc.

The Commuted Equations
As discussed in the introduction, the Ricci coefficients and curvature components will be estimated in L 2 using the null structure and Bianchi equations respectively 19 . In order to deal with the nonlinearities some of the error terms are estimated in L ∞ on the spheres. These L ∞ bounds are obtained from L 2 estimates for higher order derivatives via Sobolev inequalities. These higher order L 2 estimates are obtained through commuting the null structure and Bianchi equations with suitable differential operators, showing that the structure of the equations are preserved, and then proceeding as for the zero-th order case. It is shown in this section that the structure of the equations are preserved under commutation.
It is also necessary to obtain higher order estimates for components of the energy momentum tensor in order to close the estimates for the Bianchi and Null structure equations. Rather than commuting the Vlasov equation, which leads to certain difficulties, these estimates are obtained by estimating components of certain Jacobi fields on the mass shell. See Section 8.
Define the set of differential operators { / ∇ 3 , r / ∇ 4 , r / ∇} acting on covariant S u,v tensors of any order 20 , and let D denote an arbitrary element of this set. These operators are introduced because of the Commutation Principle of [13]: Commutation Principle: Applying any of the operators D to any of the , ψ, T should not alter its rate of decay.
This will be shown to hold in L 2 , though until then it serves as a useful guide to interpret the structure of the commuted equations.
If ξ is an S u,v tensor field, D k ξ will be schematically used to denote any fixed In order to derive expressions for the commuted Bianchi equations in this schematic notation, the following commutation lemma will be used. Recall first the following lemma which relates projected covariant derivatives of a covariant S u,v tensor to derivatives of its components. 19 The former in L 2 on the spheres, the latter in L 2 on null hypersurfaces. 20 Note that / ∇ 3 and r / ∇ 4 preserve the rank of a tensor, whilst r / ∇ takes tensors of rank (0, n) to tensors of rank (0, n + 1).
The commutation lemma then takes the following form.
Proof The proof of the first identity follows by writing and using where the last line follows by using equation (23) to write, Similarly for the second one uses and for the third, If / R denotes the curvature tensor of (S u,v , / g), the last follows from writing, and the fact that, The above Lemma implies that the terms arising from commutation take the following schematic form, The commuted Bianchi equations can then be written as follows. where and, for p = 1, 2, p denotes some fixed sum of contractions of h, , ψ and T such that p decays, according to the p-index notation, like 1 r p . Explicitly Note the presence of the first order derivative of in 2 , whilst 1 and 2 contain only zeroth order terms. Remark 3.7 By the commutation principle and induction, it is clear that the first two terms in the error (63) preserve the structure highlighted in Proposition 3.3. In the remaining terms, it is essential that D k ψ p and D k−1 ψ p appear contracted with 2 and 2 , rather than 1 . It will be clear in the proof below that it is the special form of the operators that cause this to occur. Since, for each Bianchi pair (ψ p , ψ p ), it is the case that p ≥ p + 1 2 , the 1 · D k ψ p and 1 · D k−1 ψ p terms in E 4 [D k ψ p ] still preserve the form of the error.
Similarly looking at the error (62), it is clear that the expected r decay will be preserved from Proposition 3.3.
It will also be important later that 1 and 2 do not contain any derivatives of ψ or , whilst 2 only contains first order derivatives.
Proof of Proposition 3. 6 The proof proceeds exactly as in Proposition 3.4 of [13], though one does need to be careful since some of the quantities decay slightly weaker here. We consider only the k = 1 case. A simple induction argument completes the proof for k > 1.
Consider first the / ∇ 4 ψ p equations. Using the schematic form of the commutation formulae (59), Now the Raychaudhuri equation 22 , and the Bianchi equation for / ∇ 4 ψ p imply that, Note the cancellation. Hence where Similarly, using again the schematic expressions (59), and hence, Finally, where the Gauss equation (49) has been used for the third equality. Note also the cancellation which occurs in the first equality. Hence, The schematic expressions for the / ∇ 3 ψ p equations follow similarly.
Similarly, the commuted null structure equations can be schematically written as follows.
Proposition 3.8 (cf. Proposition 3.5 of [13]). For any integer k ≥ 1 the commuted null structure equations take the form, where, and again, Proof The proof is similar to that of Proposition 3.6, though slightly simpler as there are no terms involving / D.
p only appear multiplying terms which decay like 1 r 2 . Note also that again 1 , 2 contain no derivative terms, whilst 2 contains only first order derivatives.

The Sasaki Metric
The Lorentizian metric g on M induces a metric, g, on T M, known as the Sasaki metric, which in turn induces a metric on P by restriction. The metric on T M was first introduced in the context of Riemannian geometry by Sasaki [33]. Certain properties of this metric will be used when estimating derivatives of f later. The goal of this section is to define the metric and compute certain components of its curvature tensor in terms of the curvature of (M, g). It is then shown that trajectories of the null geodesic flow of (M, g) are geodesics in P (or more generally that trajectories of the full geodesic flow are geodesics in T M) with respect to this metric and derivatives of the exponential map are Jacobi fields along these geodesics. This fact will be used in Section 8 to estimate derivatives of f . Most of this section is standard and is recalled here for convenience.

Vertical and Horizontal Lifts
Given (x, p) ∈ T M, g (x, p) is defined by splitting T (x, p) T M into its so-called vertical and horizontal parts. This is done using the connection of g on M. Given To define the horizontal lift of v at (x, p), first let c : Extend p to a vector field along c by parallel transport using the Levi-Civita connection of g on M, It is straightforward to check this is independent of the particular curve c, as long as Given the coordinates p 1 , . . . , p 4 on T x M conjugate to e 1 , . . . , e 4 , the double null frame on M, one has a frame for T M given by e 1 , . . . , and where λ μν are the Ricci coefficients of the frame e 1 , . . . , e 4 .

Example 4.1
The generator of the geodesic flow, X , at (x, p) ∈ T M is given by The vertical and horizontal subspaces of T (x, p) T M are defined as, One clearly has the following.

Proposition 4.2
The tangent space to T M at (x, p) can be written as the direct sum Since each vector in T (x, p) T M can be uniquely decomposed into its horizontal and vertical components, the following defines g on all pairs of vectors in T (x, p) T M.

The Connection and Curvature of the Sasaki Metric
Since the Sasaki metric g is defined in terms of the metric g on M, the connection and curvature of g can be computed in terms of the connection and curvature of g. The computations are exactly the same as in Riemannian geometry. See [22].
where ∇ is the connection and R is the curvature tensor of (M, g).
where R denotes the curvature tensor of g, and R the curvature tensor of g.
The proofs of Proposition 4.4 and Proposition 4.5 follow by direct computation. See [22] and also [18]. The remaining components of R can be computed similarly but are not used here.
One important property of the Sasaki metric is the following. Proof The tangent vector to a trajectory of the geodesic flow is given by the generator X . As noted above, this is given at (x, p) ∈ T M by A trajectory of the geodesic flow takes the form (γ (s),γ (s)) where γ is a geodesic in M. Hence, by Proposition 4.4,

Curvature of the Mass Shell
Proposition 4.7 IfR denotes the curvature of the mass shell P then, if (x, p) ∈ T M and X, Y, Z ∈ T x M, the following formula for certain components ofR are true.
where V = ∂ p 3 is transverse to the mass shell P.
Proof Throughout N = ∂ p 4 + p 3 p 4 ∂ p 3 + p A p 4 ∂ p A will denote the normal to the mass shell, P, such that g(N , V ) = −2.
Each identity can be shown by first writing the curvature of P in terms of the curvature of T M, N and V . If A, B, C ∈ (T T M) denote vector fields on T M then, since,∇ where∇ is the induced connection on P, one easily deduces, To obtain the first identity note that, by Proposition 4.4, and so Finally, by Proposition 4.5, and the formula follows from Proposition 4.5.
For the second identity note that, as above,

Derivatives of the Exponential Map
Recall the definition of the exponential map (or geodesic flow) for (x, p) ∈ T M, where γ x, p is the unique geodesic in M such that γ x, p (0) = x,γ x, p (0) = p.
Derivatives of the particle density function f are estimated using the fact that derivatives of the exponential map are Jacobi fields as follows. Consider (x, p) ∈ T M and V ∈ T (x, p) T M. Using the Vlasov equation, s (x, p)), and the chain rule one obtains, exp s (x, p)).
By Proposition 4.6, s → exp s (x, p) is a geodesic in P (or in T M). Below it will be shown that J := d exp s | (x, p) (V ) is a Jacobi field along this geodesic, and moreover J (0) and (∇ X J )(0) are computed. By taking s < 0 so that exp s (x, p) lies on the initial hypersurface {u = u 0 }, this then gives an expression for V ( f ) in terms of initial data which can be estimated using the Jacobi equation. In practice it is convenient to split V into its horizontal and vertical parts.
Here X , the generator of the geodesic flow, is tangent to the curve s → exp s (x, p).
, is a Jacobi field along exp s (x, p).
Since exp 0 (x, p) = (x, p) is the identity map, is again a Jacobi field. Clearly, as above, The first derivative can again be computed, using Proposition 4.4, as follows.

Characteristic Initial Data
In Theorem 5.1 below, characteristic initial data, prescribed on the hypersurfaces {v = v 0 }, {u = u 0 }, satisfying a certain smallness condition is considered. Of course, in the setting of Theorem 1.1, such data arises as induced data on two transversely intersecting null hypersurfaces, whose existence is guaranteed by a Cauchy stability argument and an application of a result of Klainerman-Nicolò [21] on the vacuum equations (6). See Section 1.2 and Section 1.3 where this argument is discussed. Characteristic initial data for Theorem 5.1 can, however, be prescribed independently of the setting of Theorem 1.1. Suppose "free data", consisting of a "seed" S u,v -tensor density of weight −1,/ g v 0 , on [u 0 , u f ]× S 2 , a "seed" S u,v -tensor density of weight −1,/ g u 0 , on [v 0 , ∞)× S 2 , and a compactly supported function f 0 : P| {v=v 0 } → [0, ∞), along with certain quantities on the sphere of intersection S u 0 ,v 0 , are given. Here P| {v=v 0 } denotes the mass shell over the initial hypersurfaces {v = v 0 }. The characteristic constraint equations for the system (1)-(2) take the form of ordinary differential equations and can be integrated to give all of the geometric quantities , ψ, T , along with their derivatives, on {v = v 0 } and {u = u 0 } once the above "free data" is prescribed. These geometric quantities, on {v = v 0 } and {u = u 0 }, are what is referred to as the "characteristic initial data" in the statement of Theorem 5.1. Appropriate smallness conditions can be made for the "free data" (and their derivatives), along with appropriate decay conditions for the seed/ g u 0 , in order to ensure the conditions of Theorem 5.1 are met.
The prescription of such characteristic "free data", and the determining of the geometric quantities from them, will not be discussed further here. The interested reader is directed to [9], where this is discussed in great detail in a related setting. See also [13].

The Main Existence Theorem
Define the norms where the second summation is taken over The main theorem can now be stated more precisely as follows.
Theorem 5.1 There exists a v 0 large and an ε 0 small such that the following holds. Given smooth characteristic initial data for the massless Einstein-Vlasov system and the data on {u = u 0 } satisfy, Here κ and the variables are defined as certain combinations of Ricci coefficients and Weyl curvature components in Section 11, andẼ 1 , . . . ,Ẽ 7 is a frame for P defined by,Ẽ where E 1 , . . . , E 7 is a frame for P defined in Section 8.5. Suppose also that, in supp f | P| {v=v 0 } , for some fixed constants C p 1 , . . . , C p 4 independent of v 0 , and that, in each of the two spherical coordinate charts, the components of the metric satisfy, for some constant C uniformly on the initial hypersurfaces {u = u 0 }, {v = v 0 }.

Then there exists a unique solution of the Einstein-Vlasov system
where C is a constant which can be made arbitrarily small provided ε 0 and 1 v 0 are taken sufficiently small. Moreover one also has explicit decay rates for the size of supp f | P x ⊂ P x as v(x) → ∞ and explicit bounds on weighted L 2 norms of the variables. See Section 7 and Section 11 respectively. Finally, if u f was chosen sufficiently large, f = 0 on the mass shell over any point The L 4 norms of the Weyl curvature components are required for the Sobolev inequalities on the null hypersurfaces. See Section 6.

Bootstrap Assumptions
The proof of Theorem 5.1 is obtained through a bootstrap argument, so consider the following bootstrap assumptions for Ricci coefficients for k = 0, 1, 2, the spherical Christoffel symbols, for k = 0, 1, 2, for Weyl curvature components and for the energy momentum tensor components, where C is some small constant 23 . Moreover, since a derivative of b appears in the expression for / ∇ 4 e A , consider also the bootstrap assumption for an additional deriva- and also for / − / • at the top order, Recall that / − / • is a geometric object, an S u,v (1, 2) tensor, and so its covariant derivatives are well defined.
Note that, since the volume form of S u,v grows like r 2 , (66) is consistent with the expectation that p behaves like 1 r p . Moreover, (67), (71) is consistent with the expectation that | / − / • | decays like 1 r , or equivalently (by Proposition 6.3 below) that the components / C AB − / •C AB behave like 1 with respect to r . Since the / •C AB behave like 1, this implies the components / C AB also behave like 1 and hence that r | / | ≤ C where, These pointwise bounds for lower order derivatives are derived from the bootstrap assumptions (67), (71) via Sobolev inequalities in Section 6. The covariant derivatives of / are defined in each coordinate system as, Finally, note also that (70) is consistent with b = 1 . Since b is only estimated on an outgoing null hypersurface at the top order though, the Sobolev Inequalities of the next section only allow us to conclude unlike at lower orders where the Sobolev inequalities will give, r |b|, r |Db| ≤ C.
Here and throughout the remainder of the paper C will denote a numerical constant which can change from line to line.

The Bootstrap Theorem
Theorem 5.1 will follow from the following bootstrap theorem, Theorem 5.2, via a last slice argument.

Theorem 5.2
There exist ε 0 , C small and v 0 large such that the following is true. Given initial data satisfying the restrictions of Theorem 5.1, let A denote a characteristic rectangle of the form

Sobolev Inequalities
The Sobolev inequalities shown in this section will allow one to obtain L ∞ estimates on the spheres for quantities through the L 2 bootstrap estimates (66)-(71). They are shown to hold in the setting of Theorem 5.2, i. e. for x ∈ A, and are derived from the isoperimetric inequality for each sphere S u,v : if f is a function which is integrable on S u,v with integrable derivative, then f is square integrable and where denotes the average of f on S u,v , and I (S u,v ) denotes the isoperimetric constant of S u,v : where the supremum is taken over all domains U ⊂ S u,v with C 1 boundary ∂U in S u,v .
The following Sobolev inequalities are standard, see e.g. Chapter 5.2 of [9]. Lemma 6.1 Given a compact Riemannian manifold (S, / g), let ξ be a tensor field on S which is square integrable with square integrable first covariant derivative. Then ξ ∈ L 4 (S) and where I (S) := max{I (S), 1}, and C is a numerical constant independent of (S, / g).

Lemma 6.2 If ξ is a tensor field on S such that
, where again C is independent of (S, / g).
Under the assumption that the components of / g satisfy, for some constant C > 0. It follows that and hence there exist constants c, C > 0 such that Using this fact, Lemma 6.1 and Lemma 6.2 can be combined to give Thus, in order to gain global pointwise control over the Ricci coefficients, curvature components and energy momentum tensor components, it remains to gain control over the isoperimetric constants I (S u,v ). We first show the above bounds on the components of / g hold under the following bootstrap assumptions. Let A ⊂ A be the set of points x ∈ A such that, trχ for some constant C 0 , for all points y ∈ A such that u(y) ≤ u(x), v(y) ≤ v(x).

Proposition 6.
3 If x ∈ A then, in each of the two spherical charts defined in Section 2.1, Proof Recall the first variation formula (54) which implies that, and hence, This gives, and hence, using the assumption that for some constants C, c > 0, where γ is the round metric, and the bootstrap assumptions (75)-(77), Let λ and denote the eigenvalues of / g AB r 2 such that 0 < λ ≤ . There exists v = (v 1 , v 2 ) such that max{|v 1 |, |v 2 |} and, Hence, and This implies that, Returning now to the first variational formula (54), summing over A, B and using the above bounds for the components of χ this gives, Using again the bootstrap assumptions (75)-(77) and the fact that, the Grönwall inequality implies, The first result then follows from the fact that, The result for / g AB follows from this and the bounds on det / g above.
If ξ is a (0, k) S u,v tensor such that |ξ | ≤ C r p , then this proposition implies that the components of ξ satisfy This fact will be used in Section 8, together with the bootstrap assumptions and Sobolev inequalities, to give bounds on the components of the Ricci coefficients, Weyl curvature components and energy momentum tensor components.

Lemma 6.4 If x ∈ A then, for u
so that the constant I (S u,v ) in Lemma 6.1 and Lemma 6.2 is equal to 1.
Proof The proof proceeds as in Chapter 5 of [9].
Combining the equation (74) for k = 0, 1. In particular, by taking C to be sufficiently small, the bootstrap assumptions (75)-(77) can be recovered with better constants. Hence A ⊂ A is open, closed, connected and non-empty, and therefore A = A. Note that, provided C is taken to be sufficiently small, this implies that This fact will be used throughout. Finally, to obtain pointwise bounds for curvature on the spheres from bounds on , an additional Sobolev inequality is required.

Lemma 6.5 If ξ is an S u,v tensor then, for any weight q,
. Lemma 6.5 together with Lemma 6.2, equation (73) and the bound on the isoperimetric constant combine to give the inequalities, +|r / ∇ξ | 2 + |ξ | 2 dμ S u,v dv 1 2 , and sup u 0 ≤u ≤u for any weight w. The bootstrap assumptions (68) then give the following pointwise bounds on curvature, sup for k = 0, 1. Finally, note also that, whilst (78) give the pointwise bounds r |D k b| ≤ C, for k = 0, 1, the bootstrap assumption (70) together with Lemma 6.5 give the additional pointwise bounds, and (67), (71) give, for k = 0, 1, as discussed at the end of Section 5.

Geometry of Null Geodesics and the Support of f
The decay of the components of the energy momentum tensor come from the decay of the size of the support of f in P x as r (x) → ∞. The estimates on the size of the support of f form the content of this section. It will also be shown that, provided u f is chosen suitably large, the matter is supported to the past of the hypersurface {u = u f − 1}.
Recall that the results of this section are shown in the setting of Theorem 5.2, so that they hold for points x ∈ A. Throughout this section γ will denote a null geodesic emanating from {v = v 0 } (so that v(γ (0)) = v 0 ) such that (γ (0),γ (0)) ∈ supp( f ). The tangent vector to γ at time s can be written with respect to the double null frame aṡ γ (s) = p A (s)e A + p 3 (s)e 3 + p 4 (s)e 4 .
Note that in the next section notation will change slightly (γ (0) there will be a point in {v > v 0 } ∩ π(supp( f )) and the parameter s will always be negative).

Proposition 7.1 For such a geodesic γ ,
for all s ≥ 0 such that γ (s) ∈ A, provided v 0 is taken suitably large.
The proof of Proposition 7.1 is obtained through a bootstrap argument, so suppose s 1 ∈ (0, ∞) is such that, for all 0 ≤ s ≤ s 1 . Clearly the set of all such s 1 is a non-empty, closed, connected subset of (0, ∞). The goal is to show it is also open, and hence equal to (0, ∞), by improving the constants.
The following fact, proved assuming the above bootstrap assumptions hold, is used for integrating the geodesic equations in the proof of Proposition 7.1.

Proof of Proposition 7.1 The geodesic equation for p 4 ,
written using the notation for the Ricci coefficients introduced in Section 2.3 takes the formṗ Using the fact that the pointwise bounds for give, the bootstrap assumptions (80)-(82) then imply that, Hence, for any s ∈ [0, s 1 ], Taking v 0 , and hence r (0), sufficiently large then gives, improving the bootstrap assumption (80). Consider now the geodesic equation for p 3 (s), Recalling that,ṙ this then gives, where the mass shell relation (21) has been used to obtain the cancellation. Inserting the pointwise bounds for the components of and the bootstrap assumptions (80)-(82), this gives, Again, integrating from s = 0 gives, and using the assumption on p 3 (0), If v 0 , and hence r (0) is sufficiently large this gives, hence improving the bootstrap assumption (81). Finally, consider the geodesic equation for p A (s), for A = 1, 2, which similarly gives, and using the bootstrap assumptions (80)-(82), Integrating and again taking v 0 large similarly gives, improving the bootstrap assumption (82). The set of all s 1 such that (80)-(82) hold for all 0 ≤ s ≤ s 1 is therefore a non-empty, open, closed, connected subset of (0, ∞), and hence equal to (0, ∞).
Finally we can show that π(supp( f )) is contained in {u ≤ u f − 1} for some u f large.

for all s ≥ 0 provided u f is chosen sufficiently large and γ (s) ∈ A.
Proof Recall thatu(s) = p 3 (s) 2 . Since 2 ≥ 1 2 , Proposition 7.1 implies that and so

Estimates for the Energy Momentum Tensor
Recall the notation from Section 3, and the set A from Theorem 5.2. The main result of this section is the following.

Proposition 8.1 If x ∈ A then, for u
and and for k ≤ 4, for some constant C.
Recall from Section 5 that ε 0 describes the size of the data. The main difficulty in the proof of Proposition 8.1, and in fact the main new difficulty in this work, is estimating derivatives of f . In Section 8.1, Proposition 8.1 is reduced to Proposition 8.2, a statement about derivatives of f . In particular, in Section 8.1 it is seen how the zeroth order bounds, r p |T p | ≤ C, are obtained using the results of Section 4. A collection of operatorsD, which act on functions h : P → R, akin to the collection D introduced in Section 3.3, is defined and used in the formulation of Proposition 8.2. In Section 8.2 additional schematic notation is introduced. This notation is used throughout the remainder of Section 8. In Section 8.3 seven more operators, V (0) , . . . , V (6) , are introduced and Proposition 8.2 is further reduced to Proposition 8.6, which involves bounds on combinations of the six operators V (1) , . . . , V (6) applied to f . The main observation is that the Vlasov equation can be used to replace V (0) f := r Hor (x, p) (e 4 ) f with a combination of operators from V (1) , . . . , V (6) (such that the coefficients have desirable weights) applied to f . The operators V (1) , . . . , V (6) , in Section 8.4, are then used to define corresponding Jacobi fields (with respect to the Sasaki metric, defined in Section 4) J (1) , . . . , J (6) . In Section 8.5 two frames for the mass shell, {E i } and {F i }, are defined. In Sections 8.6 and 8.7 bounds for the components, with respect to the frame {E i }, of the Jacobi fields, along with their derivatives, are obtained. It is transport estimates for the Jacobi equation which are used to obtain these bounds. For such estimates it is convenient to use the parallel frame {F i }, and it is therefore also important to control the change of frame matrix , also defined in Section 8.5. In Section 8.8 it is shown how Proposition 8.6 follows from the bounds on derivatives of the components of the Jacobi fields obtained in Sections 8. 6  In particular This notation will be used throughout this section. Finally, it is assumed throughout this section that x ∈ A.

Estimates for T Assuming Estimates for f
Consider the set of operators {e 3 , re 4 , re 1 , re 2 }. The notationD will be used to schematically denote an arbitrary operator in this set. These operators act on functions h : P → R on the mass shell where, for example, in the coordinate system (u, v, θ 1 , θ 2 , p 1 , p 2 , p 4 ) on P (as usual it is assumed we are working in one of the two fixed spherical coordinate charts). Given a collection of derivative operators from the set {r / ∇, / ∇ 3 , r / ∇ 4 }, say D k , this will act on (0, m) S u,v tensors and give a (0, l + m) S u,v tensor, where l ≤ k is the number of times r / ∇ appears in D k . LetD k C 1 ,...,C l denote the corresponding collection of derivative operators inD where, in D k (e C 1 , . . . , e C l ), each r / ∇ 4 is replaced by re 4 , each / ∇ 3 is replaced by e 3 , and each r / ∇ C i is replaced by re C i . So for example if k = 4 and Using the results of Section 7, the proof of the k = 0 case of Proposition 8.1 can immediately be given. First, however, Proposition 8.2, a result aboutD derivatives of f , is stated. The full proof of Proposition 8.1, assuming Proposition 8.2, is then given after Proposition 8.3, Proposition 8.4 and Lemma 8.5, which relate D derivatives of T toD derivatives of f .

Proposition 8.2 If x ∈ A then, for u
In the above, The proof of Proposition 8.2 is given in Sections 8.3-8.8.

Proposition 8.3 Given h :
where the last line is true for C = 1, 2.
Proof This follows by directly computing the derivatives of each T p . For example, hp 3 p 4 det / g p 4 dp 1 dp 2 dp 4 = 4e 4 P x hp 3 det / gdp 1 dp 2 dp 4 = 4 P x (e 4 (h) + htrχ ) p 3 det / gdp 1 dp 2 dp 4 hp 3 p A det / g p 4 dp 1 dp 2 dp 4 The other derivatives are similar.

Proposition 8.4
For any k ≥ 1, and E T p = 0.
Here l ≤ k is the number of times r / ∇ appears in D k .
Proof The proof follows by induction by writing and using the previous proposition.
Proof Recall that Using the first variation formula (54), and hence, by the Cauchy-Schwarz inequality. Similarly, using (55), and hence, Finally, hence similarly, The higher order derivatives follow similarly.
Proof of Proposition 8.1 For k = 0 the result follows using the bounds on p 1 , p 2 , p 3 , p 4 in supp( f ), and hence, since f is preserved by the geodesic flow, This fact and the above bounds imply f p 4 det / g(x) p 4 dp 1 dp 2 dp 4 ≤ Cε 0 r (x) 2 , as det / g(x) ≤ Cr(x) 2 . One then easily sees, Moreover, For first order derivatives of T p , Proposition 8.4 and the pointwise bounds on , / , r / ∇b imply that, and hence Proposition 8.2 and Lemma 8.5 imply, Similarly for the second order derivatives, the pointwise estimates on , / , r / ∇b, D , D / , Dr / ∇b imply that, and hence, by Proposition 8.3, Proposition 8.2, Lemma 8.5 and the above bounds for |DT p | therefore give, For the third and fourth order derivatives, Proposition 8.2, Lemma 8.5, Proposition 8.3, the pointwise estimates for , / , r / ∇b, D , D / , Dr / ∇b and the pointwise bounds on T , DT , D 2 T obtained above similarly give and The estimates (83) and (85) now follow using Proposition 8.2 and the bootstrap assumptions for derivatives of , r / ∇b and / .
To obtain (84) first compute, The result then follows by integrating from u 0 to u and using (85).

Schematic Notation
To deal with some of the expressions which arise in the remainder of this section it is convenient to introduce further schematic notation. Like the previous schematic notation introduced in Section 3.1, this notation will make it easy to read off the overall r decay of complicated expressions. Throughout most of this section we will consider a point (x, p) ∈ P ∩ supp( f ) and the trajectory of the geodesic flow s → exp s (x, p) through this point. The trajectory will be followed backwards to the initial hypersurface, so s will be negative. Note that where γ is the unique geodesic in M such that γ (0) = x,γ (0) = p. The expressions exp s (x, p) and (γ (s),γ (s)) will be used interchangeably. Alsoγ μ (s) and p μ (s) will both be used to denote the μ component ofγ (s) with respect to the frame e 1 , e 2 , e 3 , e 4 . Note the slight change in notation from Section 7 where γ (0) lay on the initial hypersurface {v = v 0 } and s was positive.
Recall from Section 7 that, for such a geodesic,γ 4 (s) will remain bounded in s (and in fact will be comparable toγ 4 (0)), whilstγ 1 (s),γ 2 (s),γ 3 (s) will all decay like 1 r (s) 2 . The notationγ 0 will be used to schematically denoteγ 4 (s) andγ 2 to schematically denote any ofγ 1 (s),γ 2 (s),γ 3 (s), Certain vector fields, K , along γ will be considered later. If K 1 , . . . , K 4 denote the components of K with respect to the frame 24 1 r e 1 , 1 r e 2 , e 3 , e 4 , so that, then it will be shown that, for all such K considered, K 3 will always be bounded along γ , and K 1 , K 2 , K 4 can grow at most like r (s). Therefore K 0 will be used to schematically denote K 3 and K −1 will schematically denote K 1 , K 2 , K 4 , so it is always true that, Finally, let / e, / e −1 schematically denote the following quantities, where e A denotes the S u,v one-form / g(e A , ·), for A = 1, 2. This notation will be used for schematic expressions involving the components of Weyl curvature components, Ricci coefficients and energy momentum tensor components. If ξ is a (0, k) S u,v tensor, write, to schematically denote any of the components, Note that, if r p |ξ | / g ≤ C, then the bounds, | / g AB | ≤ C r 2 imply that, where | · | here denotes the usual absolute value on R. For example, the sum of the components α AB + / T AB decays like 1 r 2 and is schematically written, where in the summation p 1 = 4, p 2 = −2 for each of the two terms. Looking at the summation on the right hand side, it is straightforward to read off that each term should decay like 1 r 2 . Similarly, if ξ is a (k, 0) S u,v tensor, write, to denote any of the components For example, allowing ψ to also denote ψ , T to also denote T etc., where in the summation now p 1 = 4, p 2 = 2 for each of the two terms. Again, the subscript of the summation on the right hand side allows us to immediately read off that the components α AB , / T AB decay like 1 to schematically denote any of the components, For example, and it can immediately be read off from the subscript of the summation that the components α A B , / T A B decay like 1 r 4 . Note that, in this notation, it is clearly not necessarily the case that no e 1 , e 2 , e 1 , e 2 appear in the expression ξ / e 0 .

Vector Fields on the Mass Shell
Consider the vectors V (1) , . . . , V (6) ∈ T (x, p) P defined by, , for A = 1, 2. The proof of Proposition 8.2 reduces to the following.

Proposition 8.6 At any point x ∈
The vectors V (1) , . . . , V (6) , together with V (0) given by, form a basis for T (x, p) P. They are preferred to the operatorsD introduced in Section 8.1 as, in view of Proposition 4.8, it is much more natural to work with vectors divided into their horizontal and vertical parts. It will be shown below that |V (i) f | is uniformly bounded for i = 0, . . . , 6.

Remark 8.7 It is not the case that |Hor
It is for this reason the term p 4 r ∂ p A also appears in V (A) . A cancellation will later be seen to occur in these two terms, so that |V (A) f | is uniformly bounded.
In Section 8.4, the vector fields V (1) , . . . , V (6) are used to define corresponding Jacobi fields, J (1) , . . . , J (6) , along exp s (x, p). The boundedness of |V (i) f | will follow from bounds on the components of J (i) for i = 1, . . . , 6. Whilst it is true that |V (0) f | is uniformly bounded, the appropriate bounds for the components of the Jacobi field corresponding to V (0) do not hold. This is the reason V (0) derivatives of f are treated separately in the proof of Proposition 8.2 below, and do not appear in Proposition 8.6. See also the discussion in Remark 8.8. Similarly, the components of the Jacobi field corresponding to p 4 ∂ p 4 do not satisfy the appropriate bounds. The r Hor (x, p) (e 4 ) term in V (4) appears for this reason. The bound on | p 4 ∂ p 4 f | can easily be recovered from the bound on |V (4) f | using the below observation that the Vlasov equation can be used to re-express r Hor (x, p) (e 4 ) f in terms of other derivatives of f .
Below is a sketch of how Proposition 8.2 follows from Proposition 8.6. (x, p) ∈ supp( f ) is fixed. The goal is to estimateD k C 1 ,...,C l f whereD ∈ {re 1 , re 2 , e 3 , re 4 }, k ≤ 4, and l ≤ k is the number of times re 1 or re 2 appears inD k . Since V (0) , . . . , V (6) span T (x, p) P, clearlyD k C 1 ,...,C l f can be written as a combination of terms of the form V (i 1 ) . . . V (i k ) f , where k ≤ k and i 1 , . . . , i k = 0, . . . , 6. It remains to check that the V (0) can be eliminated and then that the coefficients in the resulting expressions behave well. This is done in the following steps.

Proof of Proposition 8.2 Recall the point
(1) First rewriteD k C 1 ,...,C l f as whereD ∈ {e 1 , e 2 , e 3 , re 4 }, the C k are constants and p k are powers such that p k ≤ l. The terms in the sum are lower order and so, by induction, can be viewed as having already been estimated "at the previous step", so they are ignored from now on. The r l factor in the first term will vanish when the norm is taken with the metric / g and so is also ignored. (2) Rewrite eachD in terms of the vectors V (0) , . . . , V (6) defined above, (4) .
(3) In the resulting expression bring all of the coefficients of the vectors V (0) , . . . , V (6) out to the front to get where the d i 1 ...i k are combinations of h p terms and derivatives of components of , / and b. Clearly d i 1 ...i k involves at most k − k derivatives of the components of and / , and k − k + 1 derivatives of components of b. Moreover, using the bootstrap assumptions from Section 5 they are bounded with respect to r .
commute to bring the innermost to the inside. Relabelling if necessary, this gives for some (new) d j 1 ... j k as above.
Rewrite this right hand side in terms of V (1) , . . . , V (6) and repeat Step (3) to bring the coefficients to the outside of the expression. (6) Repeat steps (4) and (5) to eliminate all of the V (0) terms to leavê where the d i 1 ...i k have the correct form as above.
The result now clearly follows from Proposition 8.6.

Remark 8.8
The vector V (0) = r Hor (x, p) (e 4 ) should be compared to the vector In the context of term is not the dominant one in S. Recall also that v here is comparable to r . Given any vector V ∈ T (x, p) P, in the sections to follow it is shown that there exists a vector field J on P such that J f satisfies the Vlasov equation and J coincides with V at the point (x, p) ∈ P. For V (1) , . . . , V (6) it will be shown that the corresponding J are all of size 1 (independent of the point x) at the initial hypersurface {v = v 0 }. Whilst this is not the case for the vector field J corresponding to S, a form of the following observation was used in the proof of Proposition 8.2. The vector field J corresponding to S has a large component, but this component points in the X direction and hence vanishes when applied to f . A manifestation of this fact is brought to light through the fact that [X M , S] = X M , where X M denotes the generator of the null Minkowski geodesic flow and, with a slight abuse of notation, S now denotes the vector field S = vHor (x, p) (e 4 ) + uHor (x, p) (e 3 ) on the mass shell over Minkowski space.
Note that this observation is not specific to the massless case, i.e. the identity [X M , S] = X M is still true if X M now denotes the Minkowski geodesic flow restricted to the hypersurface P m = {(x, p) ∈ T M | p future directed, g Mink ( p, p) = −m 2 } for m > 0. A form of this observation is used in the work [16].

The Jacobi Fields
Define vector fields J (1) , . . . , J (6) along the trajectory of the geodesic flow s → exp s (x, p) by (exp s (x, p)), by the chain rule, Taking s < 0 so that exp s (x, p) lies on the mass shell over the initial hypersurface {v = v 0 }, this relates V (i) f (x, p) to intial data. By Proposition 4.8, J (i) is in fact a Jacobi field and hence J (i) (s) can be controlled using the Jacobi equation.
Note that so far the Jacobi fields are only defined along the trajectory s → exp s (x, p). Since higher order derivatives of f will be taken it is necessary to define them in a neighbourhood of the geodesic s → exp s (x, p) in P. They are in general defined differently depending on what the higher order derivatives to be taken are. When considering the quantity for 2 ≤ k ≤ 4 the Jacobi fields are extended so that as follows.

Two Frames for P and Components of Jacobi Fields
Let s * ≥ 0 be the time such that π(exp −s * (x, p)) ∈ {v = v 0 }, where π : P → M is the natural projection. The definition of the Jacobi fields of Section 8.4 imply that, for 1 ≤ k ≤ 4, and so Proposition 8.6 will follow from appropriate estimates for J (i k ) . . .
Recall from Section 2 that p 1 , p 2 , p 4 denote the restrictions of p 1 , p 2 , p 3 to P and ∂ p 1 , ∂ p 2 , ∂ p 4 denote the corresponding partial derivatives with respect to the (u, v, θ 1 , θ 2 , p 1 , p 2 , p 4 ) coordinate system for P. For every (x, p) ∈ P define the frame E 1 , . . . , E 7 of horizontal and vertical vectors for P by Recall the expressions (19), which imply, Ver (x, p) e 1 + / g 1A p A 2 p 4 e 3 , Ver (x, p) e 2 + / g 2 A p A 2 p 4 e 3 , The vectors 1 r e A , for A = 1, 2, are used rather than the vectors e A which grow like r . For J ∈ {J (1) , . . . , J (6) } let J j denote the components of J with respect to this frame. So Define also the frame F 1 , . . . , F 7 for P along s → exp s (x, p) by π(exp s (x, p)) denotes the geodesic in M, so that exp s (x, p) = (γ (s),γ (s)), and Par (γ ,γ ) denotes parallel transport along (γ (s),γ (s)).
Let denote the change of basis matrix from {F i } to {E i }, so that Note that, at s = 0, Remark 8.9 In the following, when tensor fields on P are written in components, these will always be components with respect to the frame {E i }. So if J is a Jacobi field then J i denote the components such that When writing components with respect to the parallelly transported frame {F i }, the matrix will always be used. So, Latin indices i, j will always run from 1, . . . , 7.
It will be necessary in the following sections to estimate the components of and −1 , and certain derivatives, along (γ (s),γ (s)). It is therefore necessary to derive equations satisfied by the components of and −1 .

Proposition 8.10
The components of the matrices and −1 satisfy the equations, and This can be written as the system of equations (87).
Similarly, writing This yields the system (88).

Now,
s * are assumed to be bounded pointwise by assumption. See Theorem 5.1. Higher order derivatives can similarly be expressed in terms of derivatives of components of Jacobi fields. This is discussed further in Section 8.8. It is hence sufficient to just estimate the derivatives of components, for k = 1, . . . , 4, j = 1, . . . , 7. These estimates are obtained in the next two subsections in Propositions 8.17,8.22,8.26,8.27. In Section 8.8 they are then used to prove Proposition 8.6.

Pointwise Estimates for Components of Jacobi Fields at Lower Orders
For the fixed point (x, p) ∈ P ∩ supp( f ), recall that s * = s * (x, p) denotes the parameter time s such that π(exp −s * (x, p)) ∈ {v = v 0 }. The goal of this section is to show that the components of the Jacobi fields J (1) , . . . , J (6) , with respect to the frame E 1 , . . . , E 7 , are bounded, independently of (x, p), at the parameter time s = −s * , and then similarly for the first order derivatives J (i 2 ) (J (i 1 ) j ), for i 1 , i 2 = 1, . . . , 6, j = 1, . . . , 7. Second and third order derivatives of the Jacobi fields are estimated in Section 8.7.
The estimates for the components of J = J (1) , . . . , J (6) are obtained using the Jacobi equation which in components takes the form, There are two important structural properties of the right hand side of equation (89), essential for obtaining good global estimates for the Jacobi fields. The first involves the issue of regularity. Given that derivatives of components of the energy momentum tensor appear in the Bianchi equations as error terms, it is important to estimate derivatives of T , and hence the components of J (1) , . . . , J (6) , at one degree of differentiability greater than the Weyl curvature components ψ. It is therefore important that the right hand side of equation (89) has the correct structure to allow the components of the Jacobi fields to be estimated at this level of regularity. 25 The second important property concerns the issue of decay. Since equation (89) is used to estimate the Jacobi fields globally, it is important that the right hand side is twice globally integrable. Recall thatR denotes the curvature tensor of the induced Sasaki metic on P and, by Proposition 4.7,R(X, J )X can be expressed in terms of the curvature tensor R of (M, g). In order to check that the right hand side of (89) indeed has the two above properties,R is expressed in terms of R, which is then expanded in terms of ψ, T . It is then written using the schematic notation of Section 8.2, which allows one to easily read off the decay and regularity properties. First, by Proposition 4.7, Here J h and J v are defined by | (γ ,γ ) ).
For a vector Y ∈ T γ M, T Ver (γ ,γ ) (Y ) denotes the projection of Ver (γ ,γ ) (Y ) to P. So, for each (y, q) ∈ P on the trajectory of the geodesic flow through (x, p), It is tempting to view J h , J v as vector fields on M, though this is not strictly the case as the value of J h | (y,q) , J v | (y,q) depends not only on y but also on q. Some care therefore needs to be taken here.
In view of the above discussion regarding the regularity of the right hand side of (89), the presence of the derivatives of R in the expression (90) seem, at first glance, to be bad. On closer inspection however, one sees that such terms are always horizontal or vertical lifts of derivatives of R in theγ direction. Since the components of J (1) , . . . , J (6) are estimated by integrating the Jacobi equation (89) twice in s, the fact that the potentially problematic terms, when integrated in s, in fact lie at the same level of differentiability as R can be taken advantage of. In other words, the derivatives of R appearing on the right hand side of the Jacobi equation (89) always point in exactly the correct direction so that transport estimates for the Jacobi equation (89) recover this loss. In order to exploit this fact, it is convenient to rewrite the expression (90), using Proposition 4.4, as, The above observations explain how (89) can be used to give good estimates for J (1) , . . . , J (6) from the point of view of regularity. In order to obtain global estimates however, it is also important to see that (91) has the correct behaviour in r so as to be twice globally integrable. This can be seen by rewriting (91) in terms of ψ, T and using the bootstrap assumptions (68), (69), along with the the asymptotics for p 1 , p 2 , p 3 , p 4 obtained in Section 7, being sure to allow certain components of J to grow like r . Consider, for example, just the first term Hor (γ ,γ ) R(γ , J h )γ in (91). Recall the identity, For a vector field K along γ in M, let K μ denote the components of K with respect to 1 r e 1 , 1 r e 2 , e 3 , e 4 , Using the form of the metric in the double null frame, Hence, This can schematically be written as which, as the schematic notation now makes clear, decays like r − 5 2 and hence is twice globally integrable. Note that the summation on the right hand side of (93) can actually always begin at 3, except for terms involving β. This fact is important when estimating higher order derivatives of Jacobi fields in Section 8.7 and will be returned to then.
Recall, from Section 8.2, that K −1 is used to schematically denote K 1 , K 2 , or K 4 . It is important to denote K 1 , K 2 , K 4 as such since the E 1 , E 2 and E 4 components of some of J (1) , . . . , J (6) will be allowed to grow at rate r .
Clearly, in order to use the Jacobi equation (89) to estimate the components of J (1) , . . . , J (6) , several additional points need to first be addressed. First, it is obviously important to understand how the matrix , along with its inverse −1 , behaves along (γ ,γ ). Moreover, since the terms in R (X, J )X i involving derivatives of R have to be integrated by parts, it is also important to understand how the derivative of , d ds , behaves along (γ ,γ ). An understanding of the behaviour of these matrices is obtained in Proposition 8.14 below using the equations (87), (88). Since the compo-nents of∇ X E i appear in equations (87), (88), they are written in schematic notation in Proposition 8.12.
Secondly, it is necessary to understand the initial conditions, for the Jacobi equation (89). These initial conditions are computed in Proposition 8.15. The reader is encouraged on first reading to first setR(X, J )X equal to zero (i.e. to consider the Jacobi fields on a fixed Minkowski background) in order to first understand the argument in this simpler setting. The Jacobi equation (89) can, in this case, be explicitly integrated and explicit expressions for the F 1 , . . . , F 7 components of the Jacobi fields, (J i i j )(s), can be obtained. It is clear that, even in this simplified setting, in order to obtain the appropriate boundedness statements, certain cancellations must occur in certain terms arising from J j (0) and certain terms arising from (∇ X J ) j (0) for some of the Jacobi fields. Lemma 8.11 below is used to exploit these cancellations in the general setting.
Finally, it is convenient to write some remaining quantities appearing in the expression (91) in schematic notation. This is done in Proposition 8.13.
The zeroth order estimates for the components of J (1) , . . . , J (6) are then obtained in Proposition 8.17, with Lemma 8.16 being used to make the presentation more systematic.
To obtain estimates for first order derivatives of the components of J (i 1 ) , for i 1 = 1, . . . , 6, the Jacobi equation (89) is commuted with J (i 2 ) , for i 2 = 1, . . . , 6. The fact that J (i 2 ) is a Jacobi field along s → exp s (x, p) = (γ (s),γ (s)), a curve in P whose tangent vector is X , guarantees that [X, J (i 2 ) ] = 0, i.e. J (i 2 ) commutes with d ds . It is now crucial to ensure that the schematic form of the error terms is preserved on applying J (i 2 ) , e.g. J (i 2 ) (R(X, J )X ) j must have the same, globally twice integrable, behaviour in r as (R(X, J )X ) j , for j = 1, . . . , 7. Moreover, in obtaining the zeroth order estimates, it was important that the bound was true in order that the right hand side of the Jacobi equation (89) could be twice integrated in s. It is therefore also important to also ensure that, Note that this property is completely independent of the behaviour in r . Proposition 8.18 is motivated by showing these properties of the error terms are preserved. In order for this to be so, it quickly becomes apparent that, at the zeroth order, it is necessary to show that |J (i 1 ) j | ≤ C p 4 for j = 5, 6, 7. Also, on inspection of Proposition 8.18, one sees that J (γ A ) contains terms of the form, The presence of such terms means that, in order to see that J (i 2 ) (γ A ) has the correct 1 r 2 behaviour, it is necessary to ensure that, for each J = J (1) , . . . , J (6) , either J 4+A is not merely bounded at s = −s * , but behaves like 1 r along (γ ,γ ), or that an appropriate cancellation occurs between the J 4+A and J A terms. It is hence necessary, at the zeroth order, to not just show boundedness of the components at s = −s * , but to understand their behaviour along (γ ,γ ) in more detail. To gain this understanding it also becomes necessary to understand properties of the change of frame matrices, and −1 , in more detail. See Proposition 8.14 and Proposition 8.17. In order to further commute the Jacobi equation (89), to estimate second and third order derivatives of the components of the Jacobi fields, it is also necessary to understand more detailed properties of first order derivatives of the components of the Jacobi fields.
As preliminaries to the estimates for the first order derivatives of the components of the Jacobi fields, which are treated in Proposition 8.17, relevant properties of first order derivatives of and −1 are understood in Proposition 8.20, along with relevant properties of first order derivatives of the initial conditions for the Jacobi equation in Proposition 8.21.
The following Lemma, recall, will be used to exploit certain cancellations in terms arising from the initial conditions for the Jacobi equation (89).
Proof Note thatṙ and so Recall that | 1 2 p 3 (s)| ≤ C r 2 p 4 (s), and the geodesic equation for p 4 , which implies that, Hence, Note that Lemma 8.11 in particular implies that and also, In the following two propositions, terms arising in the equations (87), (88) for and −1 , and in the expression (91) forR(X, J )X are respectively written in schematic form.
Note that the second summation in each line guarantees that terms involving Weyl curvature components and energy momentum tensor components decay slightly better than the others. This fact is important and will be returned to in Proposition 8. 23. Note also that, if r / ∇b is just regarded as D 1 , the terms involving r / ∇b above also decay slightly better. This extra decay is important at higher orders because of the weaker bounds we have for D 3 r / ∇b.
Proof Using Proposition 4.4 and the table of Ricci coefficients (24)-(28), one derives, , where, for a vector Y on M, T Ver (γ ,γ ) (Y ) denotes the projection of Ver (γ ,γ ) (Y ) to P. The terms involving the curvature R of M can be found explicitly in terms of ψ, T by setting K = e A , e 3 or e 4 in the expression (92). For example, in the above expression for∇ X E 4+A , The result follows by inspecting each of the terms and writing in schematic notation.
Using the bootstrap assumptions for the pointwise bounds on / , / ∇b, , ψ, T and the fact that r p |γ p | ≤ p 4 , Proposition 8.12 in particular gives, for i, j = 1, . . . , 7, for all i = 1, . . . , 4, j = 5, 6, 7 or vice versa, and for A, B = 1, 2. These facts are crucial for showing the schematic form of the error term in the Jacobi equation (see Proposition 8.18 below) is preserved after taking derivatives. Recall, by Proposition 7.1, for all s ∈ [−s * (x, p), 0], for some constants c, C which are independent of the point (x, p) ∈ P ∩ supp( f ). So p 4 in the above expressions can either be taken to be evaluated at time s or time 0.
Proposition 8. 13 In schematic notation, Proof Using the table of Ricci coefficients (24)-(28) one can compute, In the next proposition, estimates for the components of the matrices , −1 are obtained.
Proof The proof proceeds by a bootstrap argument. Assume, for some s ∈ [−s * , 0], that for i, j = 1, . . . , 7 and that Choose C 1 > 4 large so that C 1 > 4C, and v 0 large so that 1 The In the next proposition the initial conditions for the Jacobi equation (89) are computed.

Proposition 8.15
The Jacobi fields J (1) , . . . , J (6) , along with their first order derivatives in the X direction, take the following initial values. Hor for A = 1, 2, and J (7) Hor The expressions involving the curvature tensor R of (M, g) can be written explicitly in terms of ψ and T using the expression (92).
Proof The proof follows directly from Proposition 4.8.
The components of the Jacobi fields J (1) , . . . , J (6) can now be estimated along (γ (s),γ (s)). Recall that it is important, in order to show the schematic form of the Jacobi equation is preserved after commuting with Jacobi fields, to identify the leading order terms of some of the components. The leading order term of J (3) 3 is also identified in order to carry out a change of variables in the proof of Proposition 8.6. See Section 8.8.
The following lemma will be used.
for A = 1, 2, and Proof The proof follows by using the fundamental theorem of calculus to write, and using the estimates for the components of −1 from Proposition 8.14.

s) .
If i = 1, 2 then, using the fact that (see the proof of Lemma 8.11), Lemma 8.11 and Proposition 8.15 imply that for A = 1, 2, and that, for j = 1, . . . , 4, j = A, and for j = 5, 6, 7, j = 4 + A. Lemma 8.16 then gives, for B = 1, 2, and Similarly, using the fact that, etc., the bounds for J (4) , and, for A = 1, 2, and can be obtained. It remains to recover the bootstrap assumptions (96)-(98) with better constants. It will be shown that, for each J = J (1) , . . . , J (6) , for A = 1, 2, and where C(C 1 ) is a constant depending on C 1 . By integrating the Jacobi equation (89) and taking C small depending on C(C 1 ), the bootstrap assumptions (96) The components of this term with respect to E 1 , . . . , E 4 are exactly the components of R(γ , J h )γ with respect to the frame 1 r e 1 , 1 r e 2 , e 3 , e 4 for M. From the schematic expression (93), the pointwise bounds on the components 26 of ψ, T and the fact that r p |γ p (s)| ≤ C p 4 (s), one immediately sees that, r (s) 5 2 , and hence, by Proposition 8.14, r (s) 3 2 , for j = 1, 2, 3, 4. Other than those in the bottom line, the remaining horizontal components in the expression (91) can be treated similarly using also Proposition 8.13 and the pointwise bounds on the components of , / , / ∇b. For the term in the bottom line of (91), writê where μ runs from 1 to 4 in the summations. The horizontal components of the second term of (102) can be estimated exactly as the others using Proposition 8.12. For the components of the first term, write, for j = 1, . . . , 7. Then using again the schematic expression (93), the pointwise bounds on ψ, T and Proposition 8.14, the terms in the Jacobi equation (89) which the first term of (102) give rise to can be estimated, The vertical terms in (91) are similarly estimated as follows. Notice now that, ignoring the term∇ in the bottom line of (91) for now, each term contains at least threeγ factors and moreover that, since Proposition 8.13 guarantees that the terms involving ∇γ e α gain an extra power of decay. Similarly, since ∇γγ = 0, one can check, and hence the terms involving X (γ α ) also gain an extra power of r decay. Similarly for the vertical terms arising from the second term in (102), by Proposition 8.12. Since is quadratic in R this term also decays better. Hence, using also Proposition 8.14, one sees all the vertical terms in (91), still ignoring the final term in the bottom line, can be controlled by 27 For the final term, writê where λ runs over 1,2,4 andλ(1) = 5,λ(2) = 6,λ(4) = 7. The components of the second term can be estimated as before (with the additional r decay) by Proposition 8.12. The components of the first term can again be estimated after integrating, 27 They will actually decay like 1 r 7 2 , but 1 r 3 is sufficient.
Suppose now i 1 , i 2 = 1, . . . , 6. Since J (i 2 ) is a Jacobi field along (γ ,γ ), a curve with tangent vector field X , it is true that [X, J (i 2 ) ] = 0 and the Jacobi equation for the components of J (i 1 ) can be commuted with J (i 2 ) to give, The goal now is to repeat the proof of Proposition 8.17 to get pointwise estimates for J (i 2 ) (J (i 1 ) j ) along (γ ,γ ). It is first necessary to show that the schematic form of R(X, J (i 1 ) )X is preserved after differentiating the components with respect to J (i 2 ) . As with the K notation introduced in Section 8.2, for J = J (1) , . . . , J (6) , let the components be schematically denoted as follows, By Proposition 8.17, it is always true that for some constant C. 28
The last point is immediate from Lemma 3.4.
Note that it is the terms in the last line of (103) which make it necessary to keep track of the leading order terms in some of the Jacobi fields.
Remark 8. 19 One easily sees that the last point from Proposition 8.18 is true at higher orders, i.e., for k ≥ 1. This fact will be used later when estimating higher order derivatives of the Jacobi fields.
Using the bounds on the components of J (1) , . . . , J (6) from Proposition 8.17, Proposition 8.18 in particular guarantees that for J = J (1) , . . . , J (6) . In order to mimic the strategy used to obtain the zeroth order estimates of the components of the Jacobi fields, estimates for J (i 2 ) ( k j ) along (γ ,γ ) are first obtained, then the initial conditions for J (i 2 ) (J (i 1 ) j ) are computed in Proposition 8.20 and Proposition 8.21 respectively.
The next proposition gives pointwise estimates for the initial conditions for the commuted Jacobi equation. As was the case for the uncommuted equation, the leading order terms of some of the components have to be subtracted first.
for i = 4, and that, by Proposition 8.18, one has the same pointwise bounds for V (i 2 ) applied to the terms involving curvature as one does for the terms involving curvature alone.
for all s ∈ [s, 0], for all i 1 , i 2 = 1, . . . , 6, where C 1 is a large constant which will be for A = 1, 2, i 1 , i 2 = 1, . . . , 6 otherwise, and for all i 1 , i 2 = 1, . . . , 6, for all s ∈ [s, 0]. The remainder of the proof proceeds exactly as that of Proposition 8.17, recalling that [X, J (i 2 ) ] = 0. By Proposition 8.18 and Proposition 8.20 one has the same bounds for, the right hand side of the commuted Jacobi equation, as for the uncommitted equation since the bootstrap assumptions of Section 5 and the Sobolev inequalities give pointwise bounds for Dψ, DT , D , D / , Dr / ∇b.

L 2 Estimates for Components of Jacobi Fields at Higher Orders
To estimate , the Jacobi equation needs to commuted three and four times respectively. This will generate terms involving two and three derivatives of Ricci coefficients, Weyl curvature components and energy momentum tensor components. The higher order derivatives of the components of the Jacobi fields must therefore be estimated in L 2 . They will additionally only be estimated after integrating in momentum space, i.e. after integrating over P x .
Given p) to be the parameter time such that π(exp s v (x, p)) ∈ {v = v }, where π : P → M is the natural projection map. In this notation, where s * (x, p) is defined in Section 8.5.
The goal of this section is to show that, for all i 1 , i 2 , i 3 , i 4 = 1, . . . , 6, j = 1, . . . , 7, the quantities , by up to two and three derivatives of Ricci coefficients, curvature components and energy momentum tensor components respectively. It will then be possible to estimate the quantities (104) after taking appropriate weighted square integrals.
The case where two derivatives of the components of the Jacobi fields are taken will first be considered. Mimicking again the proof of the zeroth order estimates, second order derivatives of the matrices and −1 are first estimated, followed by estimates for second order derivatives of the initial conditions for the Jacobi equation in Proposition 8.23 and Proposition 8.25 respectively. The following Proposition should therefore be compared to Proposition 8.14 and Proposition 8.20.
By Proposition 8.18, Proposition 8.17 and Proposition 8.12, Hence, integrating equation (106) from v to v(x) and using the fact that, it follows that Now choose C 1 large so that C 1 > 4C, where C is the constant appearing in the above inequality, and v 0 large so that C 1 v ≤ 1. Then, The proof of the second part follows by making the bootstrap assumption, for all v ≤ṽ ≤ v(x), for j = 1, . . . , 4, k = 5, 6, 7 or vice versa. The proof proceeds as before, now using the fact that and for l = 1, . . . , 4, j = 5, 6, 7 or vice versa. The proof for −1 k j is identical.
The next proposition gives estimates for the initial conditions for the commuted Jacobi equation. Again, the leading order terms of some of the components have to be subtracted first.
Proof Again follows from considering expressions for V (i 1 ) ,∇ X V (i 1 ) , differentiating the components and using the fact that, for i = 4, as in Proposition 8.21.
where the terms arising from T p ( For the second term on the right hand side of (108), For the final terms, Hence, and hence, dṽ.
Integrating each term on the right hand side of (109) then gives, Taking v 0 large so that C(1+C 1 ) v 0 ≤ C 1 2 , this then recovers the first bootstrap assumption with a better constant. The other bootstrap assumptions can be recovered similarly.
where they hold is non-empty, open and closed and hence equal to Finally, at the very top order, we have the following.
The proof then follows from the considerations of Section 8.3.

Estimates for Weyl Curvature Components
The Weyl curvature components ψ are estimated in L 2 on null hypersurfaces through weighted energy estimates for the Bianchi equations. The main proposition of this section, Proposition 9.3, will show that, at any point x ∈ A (see Theorem 5.2), the bootstrap assumptions for curvature (68) can be retrieved with better constants. Each Bianchi pair is assigned a weight q, q(α, β) = 5, q(β, (ρ, σ )) = 4, q((ρ, σ ), β) = 2, q(β, α) = 0. (111) The energy estimates will be derived by integrating the following identities over a spacetime region.
Terms of the form r q h 1 |β| 2 , which would appear if any weight other than q = 2 had been chosen, would not have the correct form to be absorbed by the error term in the expression for Div r 2 |β| 2 e 4 . The proof follows by computing using the bootstrap assumptions for η, η and the upper bound for . In Lemma 9.4 below it will be shown that Repeating this for each of the identities (113),(114),(115) for k = 0, 1, . . . , s and summing then gives Proof For the sake of brevity, unless specified otherwise will denote the integral Consider first the errors in the / ∇ 3 Bianchi equations. Recall from Proposition 3.3 and Proposition 3.6 that The first term in E 3 [ψ p ] will contribute terms of the form h 1 D k ψ p to the error where 0 ≤ k ≤ k (recall that Dh 1 = h 1 ) and these can be dealt with easily The second term in E 3 [ψ p ] will contribute terms of the form D k 1 p 1 · D k 2 ψ p 2 where p 1 + p 2 ≥ p and 0 ≤ k 1 , k 2 ≤ k. Note also that, since k 1 + k 2 = k, at most one of k 1 or k 2 can be greater than 1. Assume first that k 1 ≤ 1.
Suppose ψ p = α, β, then q(ψ p , ψ p ) = 2 p − 4 and where the first line follows from the Sobolev inequality (78) (and the fact that k 1 ≤ 1) and uses p 1 + p 2 ≥ p. The third line uses the fact that q = 2 p − 4 (and recall that |α| r 2 ≤ C v 0 ). If ψ p = α or β then q(ψ p , ψ p ) = 2 p − 3 and the second term in the second line above would be r q+1 |D k ψ p | 2 which can't be controlled by the last line. The sum in the error (118), however, begins at p + 1 2 for α and β, and so in the first line in the above would have p 1 + p 2 ≥ p + 1 2 . Using this fact these terms can be controlled. If k 1 > 1 then it must be the case that k 2 ≤ 1. The above steps can then be repeated but using the Sobolev inequality (79) for D k 2 ψ p 2 . For ψ p = α, β then get where the last line now uses the fact that by the bootstrap assumption (66) and the fact that v ∼ r in the "wave zone". Similarly for ψ p = α, β. The third term in E 3 [ψ p ] will contribute terms of the form h p 1 D k T p 2 to E 3 [D k ψ p ] where 1 ≤ k ≤ k + 1 and p 1 + p 2 ≥ p. Recall that, if ψ p = α or β then actually p 1 + p 2 ≥ p + 1 2 . If ψ p = α, β then q = 2 p − 4 and, by Proposition 8.1. Similarly, if ψ p = α or β, then q = 2 p − 3 and p 1 ≥ p − p 2 + 1 2 , so r q ≤ r p+ p 2 − 7 2 and, r q h p 1 D k T p 2 · D k ψ p ≤ C r 2 p 2 −4 |D k T p 2 | 2 + r 2 p−3 |D k ψ p | 2 The final term in E 3 [ψ p ] contributes terms of the form D k 1 p 1 · D k 2 T p 2 with 0 ≤ k 1 , k 2 ≤ k, k 1 + k 2 = k and p 1 + p 2 ≥ p, or p 1 + p 2 ≥ p + 1 2 if ψ p = α or β. These terms can be dealt with as before using the fact that either k 1 ≤ 1 or k 2 ≤ 1, and the pointwise bounds for T p , DT p , D 2 T p from Proposition 8.1.
The final terms in E 3 [D k ψ p ], i.e. the terms of the form D k 1 1 · D k 2 ψ p etc. can be dealt with similarly (since all of the terms in 1 are zeroth order and the numerology for the Sobolev inequalities still work out).
The errors in the / ∇ 4 Bianchi equations can be dealt with in a similar manner. First recall that ≤ C r 2 p 2 −4 |D k 2 ψ p 2 | 2 + r 2q−2 p +1 |D k ψ p | 2 The second summation in E 4 [ψ p ], (120), will contribute terms of the form , where 1 ≤ k ≤ k + 1 and p 1 + p 2 ≥ p + 2. Since 2q(ψ p , ψ p ) − 2 p ≤ q (ψ p ), The final summation in (120) contributes terms of the form D k 1 p 1 · D k 2 T p 2 to E 4 [D k ψ p ], where 0 ≤ k 1 , k 2 ≤ k, k 1 + k 2 = k and p 1 + p 2 ≥ p + 2. These terms can be treated similarly using the fact that either k 1 ≤ 1 or k 2 ≤ 1, and the pointwise bounds for T p 2 , DT p 2 from Proposition 8.1. The remaining terms in E 4 [D k ψ p ] can again be dealt with similarly. It is important to note that 1 and 2 both contain zero-th order derivatives only of and ψ. Whilst

Transport Estimates for Ricci Coefficients
In this section the Ricci coefficients are estimated in L 2 on each of the spheres S u,v through transport estimates for the null structure equations. This is done by using the identities, which hold for any scalar function h, and with h = r 2 p−2 |D k p | 2 .
can appear on the right hand side of the estimates and be dealt with by the Grönwall The equation in the e 3 direction follows from the fact that, e 3 / •C AB = 0, and, See Lemma 4.1 of [9]. The equation in the e 4 direction can similarly be derived using the fact that, Proof Equation (126) takes the schematic form, h p 1 · D p 2 · 1 + / g The estimates for D k / − / • with k ≤ 2 then follow exactly as in Proposition 10.2 (in fact these are even easier since there are no borderline terms). The estimates for D 2 / ∇ 3 / − / • follow from applying D 2 to equation (126), and the estimates for D 2 r / ∇ 4 / − / • follow from multiplying equation (125) by r and applying D 2 .
This recovers the bootstrap assumptions (67) and the (71) for when D 3 = (r / ∇) 3 . This remaining case will be recovered in the next section.

Ricci Coefficients at the Top Order
The goal of this section is to estimate D 3 r / ∇b and (r / ∇) 3 / − / • . This will recover all of the bootstrap assumptions of Section 5. In order to do this, D 3 r / ∇ p must be estimated for most of the other Ricci coefficients p . Recall the set A from Theorem 5.2.

Propagation Equations for Auxiliary Variables
and the S u,v 1-form, where, of Section 6 give the pointwise bounds 37 + C ε 0 + Again see Chapter 7 of [9]. The identity (143) then gives, , again using |K | ≤ C r 2 . Multiplying by r 4 and using the k = 1 estimate then gives the estimate for k = 2.
For k = 3 can similarly compute ( / ∇( / ∇ξ) s ) s to get By the Sobolev inequality, Similarly, Inserting this into the above inequality and multiplying by r 6 gives the result for k = 3.
Finally, for k = 4, one similarly gets, Using the Sobolev inequality as above get, by Proposition 11.6. Similarly, and, For the remaining terms, Inserting into the above estimate and multiplying by r 8 gives the result for k = 4.
If ξ is a symmetric trace free (0, 2) S u,v tensor then it suffices to know only its divergence.
the characteristic data on {v = v 0 } if t * < u f + v 0 ), Theorem 12.1 and Theorem 12.2 imply that a smooth solution to the mixed Cauchy, characteristic initial value problem exists in the region M t * +ε for some small ε > 0. This is depicted in Figure 4.
Since the bootstrap assumptions (66)-(71) hold in u 0 ≤ u ≤ u , v 0 ≤ v ≤ v for all u , v such that u + v ≤ t * , Theorem 5.2 implies they in fact hold with the better constant C 2 . Then, taking ε smaller if necessary, by compactness of t * and continuity they will hold for all u , v with u + v ≤ t * + ε (with constant C). This contradicts the maximality of t * and hence the solution exists in the entire region.