Asymptotic Stability of Minkowski Space-Time with Non-compactly Supported Massless Vlasov Matter

We prove the global asymptotic stability of the Minkowski space for the massless Einstein–Vlasov system in wave coordinates. In contrast with previous work on the subject, no compact support assumptions on the initial data of the Vlasov field in space or the momentum variables are required. In fact, the initial decay in v is optimal. The present proof is based on vector field and weighted vector field techniques for Vlasov fields, as developed in previous work of Fajman, Joudioux, and Smulevici, and heavily relies on several structural properties of the massless Vlasov equation, similar to the null and weak null conditions. To deal with the weak decay rate of the metric, we propagate well-chosen hierarchized weighted energy norms which reflect the strong decay properties satisfied by the particle density far from the light cone. A particular analytical difficulty arises at the top order, when we do not have access to improved pointwise decay estimates for certain metric components. This difficulty is resolved using a novel hierarchy in the massless Einstein–Vlasov system, which exploits the propagation of different growth rates for the energy norms of different metric components.

, is one of the most important results in mathematical relativity. There are by now several well-established strategies to address this problem, such as the original approach of [12] or the one by Lindblad and Rodnianski [26] based on the formulation of the Einstein equations in wave coordinates. These pioneering works were generalized in different ways to more general sets of initial perturbations as well as to various Einstein-matter models [5,15,19,27,23,37,40,20,21]. On the other hand, not all Einstein-matter systems have Minkowski space as an attractor. The Einstein-dust system leads to the well known Oppenheimer-Snyder collapse for initial data arbitrarily close to Minkowski space, while the Euler equations will generally lead to the formation of shocks 1 even in the absence of coupling with gravity. A realistic matter model which is widely used in general relativity and avoids shock formation on any fixed background spacetime is that of collisionless matter considered in Kinetic theory, which, when coupled to gravity, constitutes the Einstein-Vlasov system (EVS). In the case when the individual particles in the ensemble are massive this system models distributions of stars, galaxies or galaxy clusters and constitutes an accurate model for the large scale structure of spacetime. It admits a large variety of nontrivial static solutions [29,30,4,3,22] which are potential attractors other than Minkowski space.
The study of the nonlinear stability problem for Minkowski space for the EVS was initiated by Rein and Rendall in the spherically symmetric setting [28] and recently established without symmetry restrictions for certain complementary regimes of initial perturbations [15,27]. Other stability results for the massive EVS were established in the cosmological setting [1,31].
1.2. The massless Einstein-Vlasov system. The EVS is also used to model ensembles of self-gravitating photons or other massless particles, when the corresponding mass parameter m is set to zero. The system then takes the following form, T g (f )(x, v) = 0, ∀(x, v) ∈ P, for (M, g) a Lorentzian manifold and f a massless Vlasov field. Here, T g denotes the Louville vector field and P ⊂ T ⋆ M is the fiber bundle consisting of all the future light cones of the spacetime. We refer to P also as the co-mass shell 2 . The fibre of P over x ∈ M we denote by π −1 (x) and dµ π −1 (x) is the natural volume form on π −1 (x) arising from the metric g. For a comprehensive geometric introduction to relativistic Vlasov fields, see for example [33]. While the massless system formally differs from the massive system only by changing the support of f from timelike to null vectors, the behaviour of its solutions differs substantially in several key points. The first stability result of Minkowski space for the massless EVS in spherical symmetry was established by Dafermos [13] and later generalised to the case without any symmetry assumptions by Taylor [39]. In both cases, initial data are restricted to distributions of particles with compact support in momentum variables and space. This implies in particular that the particles stay in the wave zone, while the spacetime remains vacuum in interior and exterior regions. For a global existence result in spherical symmetry without necessarily small (but strongly outgoing) initial data cf. [17]. Note that, for initial data 1 On the other hand, shock formation can be avoided in the presence of accelerated expansion [36,32,35,18]. 2 This is a small abuse of language, since the particles have no mass here.
with generic momenta, a smallness assumption is nevertheless necessarily required since the massless system does possess steady states for sufficiently large data [2].
In the present paper we consider the nonlinear stability problem of Minkowski spacetime for the Einstein-Vlasov system with massless particles without any compact support assumptions, neither for the distribution function nor for the metric perturbation. This removes any restrictions related to the semi-global features observed in [13,39] and allows for arbitrary initial particle distributions including standard Maxwellians, which are excluded by compact momentum support assumptions. Moreover, metric perturbations and matter field are coupled initially in all regions and the propagation of these general initial conditions is captured by the solutions we consider. For the metric, the spatial decay rates of the initial perturbations we consider coincide with those of [26].
1.3. The main result. Our main theorem can be summarized as follows.
Theorem 1.1. (Main theorem, rough version) Consider smooth and asymptotically flat initial data (Σ 0 ,g,k,f ), where Σ 0 ≈ R 3 , to the massless Einstein-Vlasov system which are sufficiently close to the ones of Minkowski spacetime (R 3 , δ, 0, 0). Then, the unique maximal Cauchy development (M, g, f ) arising from such data is geodesically complete and asymptotically approaches Minkowski spacetime.
For a more precise statement, we refer to Subsection 2.3. In the massive case, metric perturbations and particles travel at different speeds, in particular in a uniform sense when velocities are bounded away strictly from the speed of light. In contrast, for the massless system this decoupling does not occur, which creates substantial new difficulties 3 in comparison with the massive system. We resolve these problems by a number of new techniques in the realm of the vector-field-method for relativistic transport equations [16] discussed in the following. ?
1.4. The vector field method for transport equations and technical aspects. The vector field method for relativistic transport equations was developed recently to provide a robust technique which yields sharp estimates on velocity averages of kinetic matter in spacetimes with geometries close to Minkowski spacetime [16]. It is based on the commutation properties of complete lifts of Killing fields of Minkowski spacetime with the transport operator. The method has the additional feature to be compatible with the corresponding method for the wave equation introduced by Klainerman, which constitutes the foundation of most stability results of Minkowski spacetime. For a classical version cf. [37]. The vector field method for transport equations has in the meantime been applied successfully to the Vlasov-Nordström system [14] and the massive Einstein-Vlasov system in [15]. In a serie of papers, [6,7,8,9], the method has also been extended to the Vlasov-Maxwell system in various contexts, in particular, without the need of any compact support assumptions.
In the present paper, we apply the method to the massless Einstein-Vlasov system. In particular, we introduce fundamental improvements, which are tailored to the structure of the system in the massless case, which we will lay out in the following.
1.4.1. Null structures. The vector field method is based on the commutation properties of the transport operator T g with the complete lifts of Killing fields of Minkowski spacetime. The pertubation of the transport operator, defined loosely by the difference between the transport operator in curved space and that of Minkowski spacetime, T g − T η , creates an error term in the commutator with the complete lifts and in turn obstructing terms in the resulting energy estimates.
The first crucial structure in the transport part of the massless system is the null structure of the perturbation terms. There are roughly three distincts sources of null structures. Two of them arise from the decomposition of the metric components and the momentum variables with respect to a null frame. The third arises from the identification of null forms for products involving (t, x)-derivatives of the metric components and vderivatives of the Vlasov field. These null structures are all discussed on Subsection 2.4.2.
It can be shown, as for the Vlasov-Maxwell system [8], that this structure is conserved under commutation with complete lifts. What is crucial in a subsequent step is to assure that this null structure can be exploited at all levels of regularity, which is not straightforward to validate. A particular difficulty occurs when well-behaved components of the metric perturbation need to be estimated in energy. In that case the bulk energies of Lindblad and Rodnianski are insufficient to close the estimates. We return to this issue below.

1.4.2.
A null structure in the energy-momentum tensor and its consequence for propagation of the metric perturbation. The energy momentum tensor for massless particles is tracefree. As a consequence of that, the 4-Ricci tensor is proportional to the energy-momentum tensor. From the aforementioned null structure in the momentum components, after decomposition on a standard null frame, we obtain a system of wave equations where certain matter source terms enjoy improved decay in comparison with a generic energymomentum tensor term. This structure is another characteristic feature of the massless system. To our knowledge, in the massive case, matter source terms are usually taken of the generic type and an underlying hierarchy was never exploited.
To derive suitable energy estimates for the frame components of the metric, we consider additional energy norms for the metric components. The resulting estimates are better than the generic ones due to the fast decaying matter source terms and improved null properties satisfied by the semi-linear terms of the Einstein equations. It is those energy norms that in turn can be used to estimate the good frame components of the metric perturbation when the source terms in the Vlasov equation are analysed at top order. Moreover, compared to the proof of Lindblad-Rodnianski [26], this allows us to avoid the use of Hörmander's L 1 − L ∞ estimate. 1.4.3. Strong (t − r)-decay for velocity averages. In order to close the energy estimates for the particle density, we have to deal with the weak decay rate of the pertubation part of the metric in the interior of the light cone. In the case of Vlasov fields with compact support, massless particles will follow straight lines parallel to the light cone, so that the support of the Vlasov field is located close to the light cone. We capture this effect in the non-compactly supported case using hierarchized weighted-energy norms for the Vlasov field, similar to those considered in [9]. The extra weights allows us to prove strong decay away from the wave zone, i.e. when t − r is large. 1.4.4. The Lie derivative. As in [27], we commute the Einstein equations with Lie derivatives. Following a strategy initially developped for the Vlasov-Maxwell system in [6], we also write the error terms arising in the commutation of the Vlasov equation in terms of Lie derivatives of the metric components. Compared to [15], this reduces the complexity of the error terms, and fully conserves the null structure of the system after commutation, which appears to be crucial in our proof. Moreover, it also allows to avoid many hierarchies considered in [26] in the commuted Einstein equations and in [15] in the commuted Vlasov equation. 1.4.5. Decay loss and v-derivatives. At the linear level, derivatives in v do not commute well with the massless transport operator, so that one should expect that the presence of terms of the form ∂ v i Z I f in the source term of the Vlasov equation to be problematic. In the massive case [15,27], the introduction of improved commutators seemed necessary to deal with the similar issue. Here, this issue can be resolved essentially by using the null structure of the system, the strong decay in t − r of the Vlasov field and a hierarchy of growth in t at the top order. 2. Strategy of the proof and outline of the paper 2.1. The Cauchy problem in wave coordinates and initial data. It is well known that the Einstein equations can be formulated as a Cauchy problem and in the case of the Einstein-Vlasov system, the well-posedness is guaranteed by a theorem of Choquet-Bruhat [11]. See also [38] for the massless case. A detailed formulation of the Cauchy problem for the Einstein Vlasov system can be found in [31].
Consider a smooth 3-dimensional manifold Σ 0 with a Riemannian metricg, a symmetric covariant 2-tensork and a functionf defined on T Σ (or equivalently on T ⋆ Σ), with all data assumed to be smooth and such that the constraint equations (see [31] for details) are satisfied. The Cauchy problem then consists in constructing a 4-dimensional manifold M with Lorentz metric g, a smooth function f defined on P, satisfying the Einstein-Vlasov system (1.1), and an embedding i : Σ → M such that i * g =g, i * k =k, f •pr −1 Σ =f , where k is the second fundamental form of i(Σ) in (M, g) and the function pr Σ : π −1 (Σ) → T ⋆ Σ, with π : P ⊂ T ⋆ M → M the canonical projection, is defined analogously to [31,Definition 13.30], i.e. prΣ projects p ∈ π −1 (Σ), for some hypersurfaceΣ ∈ M, to the part p ⊥ of p being perpendicular to the unit normal vector ofΣ.
In view of the decomposition (2.2), the equations (2.3) can be rewritten as a system for the components of h 1 , with extra source terms depending on h 0 . Thus, the unknowns of the reduced Einstein-Vlasov system are h 1 and the distribution function f . The initial data will be chosen small in the sense that the mass parameter M and certain energy norms of h 1 and f are bounded by a small constant ǫ > 0.
Moreover, there exists a global system of wave coordinates (t, x 1 , x 2 , x 3 ), and a constant 0 < δ(ǫ) < γ 20 , with δ(ǫ) → ǫ→0 0, in which the following energy bounds hold. For the Vlasov field, ∀ t ∈ R + , For the metric perturbation h 1 , ∀ t ∈ R + , and for all (t, Remark 2.3. At the top order, the strong growth on the energy norm of f leads to a strong growth on the L 2 norm of the pertubation of the metric. For a technical reason and in order to avoid a much stronger decay hypothesis on h 1 (0, ·), we, in some sense, include this strong growth through the weight (1 + t + r) −1 into the top order energy norm of h 1 .
The proof of the main theorem is based on vector field methods and a continuity argument so that it essentially consists in improving bootstrap assumptions on well-chosen energy norms of h 1 and f . The global-in-time existence then follows by standard arguments. As we use vector fields method, we then need to • commute the equations by high order derivatives composed by elements of K for the Einstein equations and P 0 for the Vlasov equations. • Perform energy estimates in order to propagate weighted L 2 norms of h 1 and weighted L 1 norms of f . • Obtain pointwise decay estimates on the solutions through Klainerman-Sobolev type inequalities. • Estimate all the error terms arising from the energy estimates using the decay estimates. As is usual for these type of problems, the main sources of difficulty arise from • the bad behaviour near the light cone and the weak decay rate of h 1 in the interior region t > r, • the bad commutation properties of the Vlasov equation, in particular, generating error terms containing ∂ v derivatives of f , • the top order estimates, where some of the structural properties of the equations cannot be used anymore. We present below some key technical ingredients of the proof that addresses in particular the above issues.
2.4. L 1 estimates for the Vlasov field.
2.4.1. Naive estimate. As Z, the complete lift of a Killing vector field 5 Z, commute with the flat relativistic transport operator T η := |v|∂ t + v i ∂ v i and since |g − η| is expected to be small, commuting T g (f ) = 0 with Z should create controllable error terms. However, a naive estimate leads to and, during the proof, we will have x f ||v|dvdxdτ +better terms. 5 The case of S, which is merely a conformal Killing vector field, is slightly different but do not create more complicated error terms.
Controling the left-hand side is necessary to close the energy estimates for f using a Grönwall type inequality. However, with the above naive estimate, there are two obstacles preventing us to do so (1) The decay rate degenerates near the light cone t = r. As mentionned earlier, we will deal with this issue by taking advantage of the null structure of the equations. (2) The decay rate is not integrable (and not even almost integrable). Even if we could transform the t − r decay into a t + r one, the overall t decay is too weak to derive an estimate such as Zf L 1 x,v ǫ(1 + t) η for any Z ∈ P 0 , with η ≪ 1.

2.4.2.
The null structure of the Vlasov equation. Let us denote g −1 −η −1 by H and v 0 +|v| by ∆v. Then, the deviation of T g from the flat relativistic transport operator is

Now, recall
• that the derivatives of H tangential to the light cone can be compared to those of h and have a better behavior than the others. More precisely, It will be important to notice that a similar property hold for |Lf |. • from [26,Section 8] and the wave gauge condition that the LT components of H enjoy improved decay estimates near the light cone, We will prove that ∇ e A (H) LL decay even more faster near the light, which will be crucial in our proof. • from [6, Proposition 2.9], that certain null components of v behaves better than others. In particular, in the flat case where v 0 = −|v|, one can control 6 (1 + |t − r|) 9 8 |v L |dvdxdτ by the initial energy of | Zf |, so that, in the presence of v L , we can exploit the decay in t − r in order to close the energy estimates. Moreover, the angular components satisfy, still in the flat case, |v A | |v||v L |, so that angular components also behave better than generic ones. [15,Subsection 4.2], that ∆v satisfied a kind of null condition. In our case, Now note that a naive estimate of (2.11) gives us The exponent 9 8 appearing in the denominator could be replaced by any number a > 1.
whereas, expanding all the error terms according to a null frame and taking advantage of the improved properties satisfied by the good null components of the solutions, we obtain This last estimate is much better since either the decay rate is almost integrable for t ≈ r or the Vlasov field is multiplied by |v||v L |, which allows to use part of the decay in t − r. This indicates how important the structure of the non-linearities is and how important it is to conserve them by commutation. By differentiating the metric by Lie derivatives, we will obtain that 7 which improves the commutation formula obtained in [15], where the quantities controlled, Z(h µν ), are not geometric, and where the full structure of the non-lineraties were not preserved. This will allow us to improve our naive estimate (2.10) in the following way x f ||v L |dvdxdτ + better terms, (2.14) so that we can expect to propagate the bound Zf (t, ·) L 1 x,v ǫ(1 + t) η , with η ≪ 1 independant of δ, provided that we can improve the decay in t − r of the velocity averages of f and its derivatives. Note that we will take η = δ 2 during the proof.

2.4.3.
Dealing with the non integrable decay rate. Even after exploiting the null structure as explained above, we are still left with error terms which are not time-integrable and therefore with energy norms a priori growing in time. We will circumvent this difficulty by following the strategy of [9] and we will then consider hierarchized weighted L 1 norms. It essentially relies on the following two properties.
(1) The translations ∂ µ , when applied to solutions of a wave equation, provide an extra decay far from the light cone compared to the other commutation vector fields. In view of (2.12)-(2.13), we can expect the following improved behavior for T g (∂ x µ f ), which would considerably improved the estimate (2.14) for Z = ∂ x µ . Since the worst source terms of T g ( Zf ), for any Z ∈ P 0 , contains only standard derivatives ∂ t,x f of the particle density, the system composed by the commuted Vlasov equations is in some sense triangular.
(2) The weight m := 1 can be used in order to obtain stronger decay on f . It essentially 8 arises from the contraction of the Morawetz 7 The commutations formula for the scaling and the Lorentz boosts contain more terms which can be handled in a similar way than those of (2.11). 8 The overall exponent 1/4 is here only for homogeneity, so that m ∼ t, for t ≫ r.
conformal Killing vector field K = (t 2 + r 2 )∂ t + 2tr∂ r with the flat velocity current. It satisfies in particular so that one can expect T g (m n f ) to be small and then propagate L 1 norms of f weighted by m n .
As a consequence of these two observations, we will then be able to prove an estimate such as m This will then allow us to improve the estimate (2.14) by dvdxdτ + better terms, and then prove Zf (t, ·) L 1 x,v ǫ(1 + t) η . Since we will have to consider higher order derivatives, in order to apply this strategy, we will rather consider energy norms of the form m Q− 2 3 I P Z I f (t, ·) L 1 x,v , with Q > 0 sufficiently large and where I P is the number of homogeneous vector fields composing Z I .

2.5.
Study of the metric perturbation h 1 . As already observed by Lindblad [24], differentiating the metric by Lie derivatives considerably simplifies the study of the Einstein equations. In particular for the two reasons presented here.
2.5.1. The wave gauge condition is preserved by commutation with L J Z , where Z J ∈ K |J| . More precisy, the wave gauge condition g x ν = 0 leads to and one can prove (see Subsection 4.2) that this property is preserved by differentiation by the Lie derivative, i.e.
This implies in particular, with ∇ := (∇ L , ∇ e 1 , ∇ e 2 ) containing the good derivatives of our null frame (those tangential to the light cone), that for any |J| ≤ N , In [26] (and in [15]), this property was obtained for ∇h but could not be directly obtained for its derivatives, since the quantities controlled, Z I (h µν ), were not geometric. For the purpose of this article, it is crucial to derive improved estimated on the null components of the higher order derivatives of h in order to close the energy estimates. Otherwise, certain error terms of the commuted Vlasov equations would lack too much t + r decay.
Remark 2.4. In [26], a lack of (t + r) δ -decay in the error terms of the commuted Einstein equations was circumvented by considering several hierarchies so that ∇Z I h 1 µν (t, ·) L 2 ǫ(1 + t) δ |I| , with δ |I| ≪ 1 growing with |I|. In our case the lack of decay seems to be much worse (recall the naive estimate (2.11)) and this prevents us to consider such hierarchies between the energy norms at top order.

Remark 2.5. Several analogies exists between the Eintein equations and the Maxwell equations
where the electromagnetic field F is a 2-form, * F is its Hodge dual and the source term J is a current. In particular, studying the Eintein equations in wave coordinates has to be compared to considering the Maxwell equations in the Lorentz gauge. This means that we work with a potential A satisfying dA = F and the Lorentz gauge condition ∇ µ A µ = 0, which has to be compared to the wave gauge condition since it gives |∇(A) L | |∇A|. Moreover, we noticed in [6] that so that commuting with L Z conserves the Maxwell equations as well as the Lorentz gauge condition.
2.5.2. The null structure of the Einstein equations. For the study of the Einstein equations (2.3), all the error terms arising after commutation will have enough decay outside from the wave zone. To control the error terms near the wave zone, one of course, needs to exploit the null structure and the weak null structure of the equations. Indeed, one cannot propagate L 2 estimate on h 1 by performing naive estimates. It was shown in [26] that F µν (h)(∇h, ∇h) is composed of cubic terms which decay strongly, of quadractic terms Q µν (∇h, ∇h), which are a linear combination of standard null forms, and other quadratic terms P (∇ µ h, ∇ ν h) which contains semi-linear terms satisfying Since the wave gauge condition holds, the problem comes from the term |∇h| 2 T U . To deal with it, the proof of [26] used the L 1 − L ∞ estimate of Hörmander which gave that |∇h| T U ǫ(1 + t) −1 . We provide in this paper an alternative way for treating this issue, which seems in fact necessary in order to deal with the top order energy estimates for the Vlasov field (see Subsection 2.6). The L 2 bound that we will have on h 1 is We then observe that for any (T, U ) ∈ T × U , P (∇ T h, ∇ U h) satisfies the null condition and that T [f ] T U , due to the presence of the good component v T in the integrand, decay much faster near the light cone than |T [f ]|. As a consequence, we will be able to prove that where κ ≪ 1 can be choosen independently of δ, allowing us to control sufficiently well the error term |∇h| 2 T U . During the proof, we will take κ = δ.
Remark 2.6. These estimates reflect that, even estimated in L 2 , |∇h 1 | T U has a better behavior than ∇h 1 for t ≈ r. As no improvement can be obtained far from the light cone, this property can only be captured if the L 2 norm of |∇h 1 | T U carries a weaker weight in t − r than the one of ∇h 1 .
Again, it is then important to prove that the structure of the source terms of the Einstein equations are conserved by commutation with L J Z . As noticed in [24], we have for a Killing vector field 9 Z, Moreover, the structure of the commutator is also preserved by the action of L J Z and the cubic terms as well as g h 0 µν can be easily handled. Similarly, one can prove that so that L Z (T [f ]) enjoys the same improved properties as T [f ] in the good null directions.
2.6. The top order estimates. After commuting the Vlasov equation by Z I , with |I| = N and where N is the maximal number of commutation, a specific difficulty appears with the error terms of the form where all the null structure is contained in the h 1 -factor. Since |I| = N , one cannot gain t + r decay by expressing the good derivatives ∇ in terms of the commutation vector fields anymore. Since the estimate Then, even the energy bound E 2γ,1+γ would not allow us to close the energy estimates at top order. Indeed, we would obtain Z I f (t, ·) L 1 . For a technical reason and even though |T [ Z I f ]| T U has a good behavior, this will prevent us to prove a better estimate than E 2γ,1+γ Since δ > 0, we could not improve all the bootstrap assumptions. The idea then is to remark that g (L I Z h 1 ) LL strongly decay near the light cone, so that one can propagate the bound where η 0 ≪ 1 can be choosen independantly of all the other bootstrap assumptions.

2.7.
Organization of the paper. In Section 3, we introduce the notations used in this article. Useful results for the analysis of the null structure of the equations concerning the commutation vector fields, the velocity current v and the weights preserved by the free transport operator are presented. We also introduce the energy norms used to study the solutions. In Section 4, we study the consequences of the wave gauge condition and the source terms of the commuted Einstein equations. Section 5 is devoted to the commutation formula of the Vlasov equation, as well as its analysis and in Section 6, we compute the derivatives of the energy momentum tensor T [f ]. The energy estimates used for the metric pertubation are proved in Section 7 and the one for the particle density is derived in Section 8. We set-up the bootstrap assumptions in Section 9. In Section 10, we prove pointwise decay estimates on the null components of h 1 and its derivatives and we use them to bound all the source terms of the Einstein equations but for the contribution of T [f ] in Section 11. In section 12 (respectively Section 13), we improve the bootstrap assumptions on h 1 (respectively f ). Finally, in Section 14, we prove the required estimates on the L 2 norm of T [f ] in order to close the energy estimates.

Preliminaries
In this section, we set-up the problem and introduce basic mathematical tools and notations.
3.1. Basic notations. We will use two sets of coordinates on R 1+3 , the Cartesian (t, x 1 , x 2 , x 3 ), in which the metric η of Minkowski spacetime satisfies η = diag(−1, 1, 1, 1), and null coordinates (u, u, ω 1 , ω 2 ), where and (ω 1 , ω 2 ) are spherical variables, which are spherical coordinates on the spheres (t, r) = constant. These coordinates are defined globally on R 1+3 apart from the usual degeneration of spherical coordinates and at r = 0. We will use the notation ∇ for the covariant differentiation in Minkowski spacetime. We denote by / ∇ the intrinsic covariant differentiation on the spheres (t, r) = constant and by (e 1 , e 2 ) an orthonormal basis of their tangent spaces. Capital Roman indices such as A or B will always correspond to spherical variables. The null derivatives are defined by With respect to the null frame {L, L, e 1 , e 2 }, the Minkowski metric has the following components We define further ∇ = (∇ L , ∇ e 1 , ∇ e 2 ), the derivatives tangential to the light cone, as well as U = {L, L, e 1 , e 2 }, T = {L, e 1 , e 2 } and L = {L}, which will be useful in order to study the behavior of certain tensor fields in null directions. For that purpose, we introduce for a (0, 2)-tensor field of cartesian components k αβ , If V = W = U , we will drop the subscript U U . For instance, |k| := |k| U U . As we study massless particles, the functions considered in this paper will not be defined for v = 0 so we introduce R 3 v := R 3 \ {0}. We will use the notation D 1 D 2 for an inequality such as D 1 ≤ CD 2 , where C > 0 is a positive constant independent of the solutions but which could depend on N ∈ N, the maximal order of commutation, and fixed parameters (δ, γ,...). We will raise and lower indices using the Minkowski metric η. For instance, x µ = x ν η νµ and, for a current p, The only exception is made for the metric g, where in this case, g µν will denote the (µ, ν) component of g −1 .
Finally, we extend the Kronecker symbol to vector fields, i.e. if X and Y are two vector fields, δ Y X = 0 if X = Y and δ Y X = 1 otherwise.

3.2.
Vlasov fields in the cotangent bundle formulation. Our framework for the study of the Vlasov equation and the Vlasov field is adapted from the one developed in [16] and is thus based on the co-tangent formulation of the Vlasov equation. The presentation below follows closely that of [16], but takes into account the fact that we consider here massless particles only. Let (M, g) be a smooth time-oriented, oriented, 4-dimensional Lorentzian manifold. We denote by P the following subset of the cotangent bundle T ⋆ M P := (x, v) ∈ T ⋆ M : g −1 x (v, v) = 0 and v future oriented . Note in particular that for v to be a future oriented covector, necessarily v = 0. P is a smooth 7-dimensional manifold, as the level set of a smooth function.
In the massive case, P is often referred to as the co−massshell. By an abuse of language, we will keep calling P the co-massshell, even in the present massless case. We will denote by π the canonical projection π : P → M. Given a coordinate system on M, (U, x α ) with U ⊂ M , we obtain a local coordinate system on T ⋆ M, by considering the coordinates v α conjugate to the x α such that for any We now assume that there exist local coordinates (x α ) such that x 0 = t is a smooth temporal function, i.e. its gradient is past directed and timelike. In that case, the algebraic equation It follows that (x α , v i ), 1 ≤ i ≤ 3 are smooth coordinates on P and for any x ∈ M, (v i ), 1 ≤ i ≤ 3 are smooth coordinates on π −1 (x). Note that the requirement that v = 0, implies that v i ∈ R 3 \ {0}. We thus define R 3 v := R 3 \ {0}. All integrations in v can be performed using the (v i ) coordinates in which case, the domain of integration will always be R 3 v . With respect to these coordinates, we introduce a volume form dµ π −1 (x) on π −1 (x) defined by For any sufficiently regular distribution function f : P → R, we define its energymomentum tensor as the tensor field For the above integral to be well-defined, one needs f (x, ·) to be locally integrable in v, to decay sufficiently fast in v as |v| → +∞, as well as |v|f to be integrable near 0, in view of the fact that the volume form dµ π −1 (x) becomes singular near v = 0. All distribution functions considered in this paper will always be such that these properties hold. Moreover, we will also require f to possess additional decay in x and v, so that we can perform the various intergration by parts needed. In any case, one can assume for simplicity for the computations to hold that all distribution functions are smooth, compactly supported, with a support away from v = 0, and then use the standard approximation arguments to obtain the results in the non-compactly supported case.
The Vlasov field f is required to solve the Vlasov equation, which can be written in the (x α , v i ) coordinate system as It follows from the Vlasov equation that the energy-momentum tensor is divergence free and more generally, for any sufficiently regular distribution function k : P → R, 3.3. The system of equations. We decompose the metric as is the Schwarzschild-part, and χ : R → R is a smooth cutoff function such that χ(s) = 0 if s ≤ 1 4 and χ(s) = 1 if s ≥ 1 2 . For the inverse metric we will use the decomposition The relation between h 1 and H 1 is made precise in Section 4.1. Define the reduced wave operator The massless Einstein-Vlasov system then reads and the co-mass shell condition where P (∇ µ h, ∇ ν h), Q µν (∇h, ∇h) and G µν (h)(∇h, ∇h) are (0, 2)-tensor fields, the indices (µ, ν) refers to their components in the wave coordinates system (t, x), and P, Q, G are defined as follows.
• P contains the source terms which do not satisfy the null condition and is given by • Q is a combination of the standard null forms and is given by • Finally, G µν (h)(∇h, ∇h) contains cubic terms and can be written as a linear combination of where all the indices are taken in 0, 3 . The null structure of the quadratic terms are of fundamental importance and is described in the following result.
Proof. According to (3.5) and since η LL = η LA = 0, we have for any (V, W ) ∈ U 2 , This implies all the inequalities which concern P (∇k, ∇q). Note now that, for any cartesian component (µ, ν), Q µν (∇k, ∇q) can be written as linear combination of where at least one of the λ i is equal to µ or to ν and are the standard null forms. They satisfy (see [34,Chapter 2] for a proof), for any α < β, 3.4. Commutation vector fields for wave equations. Let P be the generators of the Poincaré algebra, i.e. the set containing • the rotations • the hyperbolic rotations which are Killing vector fields of Minkowski spacetime. We also consider K := P ∪ {S}, where S = x µ ∂ µ is the scaling vector field which is merely a conformal Killing vector field. The elements of P are well known to commute with the flat wave operator η = −∂ 2 t + ∂ 2 1 + ∂ 2 2 + ∂ 2 3 and we also have [ η , S] = 2 η . We consider an ordering on K = {Z 1 , . . . , Z 11 } such that Z 11 = S and we define, for any multi-index J ∈ 1, 11 n of length n ∈ N * , Z J = Z J 1 . . . Z Jn . By convention, if |J| = 0, When commuting the system (3.4a)-(3.4b), we will use the Lie derivative to differentiate the metric g in order to preserve the structure of the equations. In coordinates, the Lie derivative L X (k) of a tensor field k α 1 ···αn β 1 ···βm with respect to a vector field X is given by . Note that that for n ∈ N, we have the equivalence relation The following standard lemma can be obtained using where C ij A are bounded smooth functions of (ω 1 , ω 2 ), and The purpose of the following result is to generalize Lemma 3.2 to tensor fields.
With this notation, we claim that for V ∈ {L, T , U } and V ∈ V, Using (3.16), we obtain, as 1 + t + r r on {r ≥ 1+t 2 },, where V, W ∈ {U , T , L}. It then only remains to bound |∇(k V W )| and ∇(k V W ) . Start by noticing that, by Lemma 3.2, so that, using (3.17) and that 1 + t + r r on {r ≥ 1+t 2 }, The following two results will be useful in order to commute the Einstein equations geometrically.
Lemma 3.4. Let k be a (0, 2) tensor fields, so that ∇k and ∇∇k are respectively (0, 3) and (0, 4) tensor fields of cartesian components For all Z ∈ K, we have Proof. Both relations follow from (3.8) and the fact that ∂ α Z β is constant for any (α, β) ∈ 0, 3 2 and Z ∈ K. Let us give more details for the first one. For cartesian components (α, µ, ν), we have To derive the equality ∇L Z k = L Z ∇k, it only remains to remark that ∂ σ ∂ ρ Z λ = 0 for all 0 ≤ σ, ρ, λ ≤ 3.

3.5.
Analysis on the co-tangent bundle. As in [16], we will commute the Vlasov equation using the complete lift Z of the Killing vector fields Z ∈ P of Minkowski spacetime. They are given by and they commute with the flat massless relativistic transport operator T η : [16,Section 2.7] for more details). Even if the complete lift S of S satisfies [T η , S] = 0, we will rather commute the Vlasov equation with S, which verifies [T η , S] = T η , for technical reason (see Lemma 3.9 below). We then introduce the ordered set where Z 11 = S and Z i = Z i if i ∈ 1, 10 , so that for any multi-index J ∈ 1, 11 n , Z J := Z J 1 . . . Z Jn . For simplicity, we will denote by Z an arbitrary element of P 0 , even if the scaling vector field S is not the complete lift of a vector field X µ ∂ x µ of the tangent bundle of Minskowski spacetime. Similarly, we will use the following convention, mostly to write concisely the commutation formula. For any Z ∈ P 0 , if Z = S, then Z will stand for the Killing vector field which has Z as complete lift and if Z = S, then we will take Z = S. The sets {Ω 12 , Ω 13 , Ω 23 , Ω 01 , Ω 02 , Ω 03 , S}, { Ω 12 , Ω 13 , Ω 23 , Ω 01 , Ω 02 , Ω 03 , S} contain all the homogeneous vector fields of K and P 0 . As suggested by Lemma 3.2, ∂ µ φ has a better behavior than Zφ for Z an arbitray element of K. It will then be important, in order to exploit several hierarchies in the commuted Vlasov equation, to count the number of homogeneous vector fields which hit the particle density f in the error terms. Given a multi-index J so that Z J ∈ K |J| and Z J ∈ P |J| 0 , we denote by J P (respectively J T ) the number of homogeneous vector fields (respectively translations) composing Z J and Z J . For instance, if Z J = ∂ t Ω 12 S∂ 2 ∂ 1 , J T = 3 and J P = 2. The following technical lemma will be in particular useful for commuting the energy momentum tensor T [f ] and then the Einstein equations. It illustrates the compatibility between the commutation vector fields of the wave equation and those of the relativistic transport equation.
x × R 3 v → R be a sufficiently regular function and Z ∈ P. Then, Proof. Let, for any Killing vector field Z ∈ P, Z w := Z − Z. We have, It then remains to note that, and, by integration by parts in v, In order to treat the curved part of the metric as pure perturbation, we define the one form Using that w U = w µ U µ = η(w, U ) for any vector field U , we directly obtain As [16], we introduce the set of weights and we consider, as suggested by [10,Remark 2.3], All the above weights are obtained by contracting the current w with the conformal Killing vector fields of Minkowski spacetime. They are preserved along the flow of T η and will be used in order to obtain strong improved decay estimates on the distribution function. In particular, m has to be compared with the Morawetz vector field (t+r) 2 2 L + (t−r) 2 2 L when used as a multiplier for the wave equation. Note that m ≤ 0, so that we often work with |m|.
We now define z as an overall positive weight, by Note also that T η (z) = 0 and moreover, since |w 0 | |v| = 1, z∈k 0 |z| |v|(1 + t + r) and |m| ≤ |v|(1 + t + r) 2 , we have The following lemma illustrates how the null components of w and the weight z interact.
Lemma 3.7. The following estimates hold from which it follows that | / w| w 0 z 1 + t + r and 1 z 1 + |t − r| .
Proof. Since w L ≤ 0 and w L ≤ 0, we have which proves the first two inequalities. For the third inequality, we use the mass shell relation for the flat spacetime from which it follows that The fourth estimate then ensues from the third and the second one. For the last inequality, we use w 0 |w L | + |w L || |w L |w 0 + |w L |w 0 and then apply the first two inequalities.
The following Lemma illustrates the good interactions between the weights z ∈ k 0 , m and the vector fields Z ∈ K.
Lemma 3.8. For all µ ∈ 0, 3 , 1 ≤ i < j ≤ 3 and k ∈ 1, 3 , we have (3.20) in order to get A straightforward computation reveals that for all z ∈ k 0 , Z ∈ P 0 , there holds Z(z) ∈ span{k 0 }, and consequently For the weight m, one can check that We then obtain the first three inequalities of the lemma by taking Y = ∂ µ , S and Ω ij in (3.22) and using (3.23)- (3.24). For the Lorentz boosts, we use the decomposition We then deduce so that, taking Y = x q r Ω 0q in (3.22) and using (3.20), (3.23) as well as (1 + |t − r|) z (see Lemma 3.7), we obtain (3.27) x q r Ω 0q (z) |m| |v|z We also obtain from (3.26) that Since |t − r| z and using that (x j w k − x k w j ) ∈ k 0 and (tw i − x i w 0 ) ∈ k 0 , we obtain from the last two equalities Combining this last inequality with (3.22), applied with Y = x j r Ω 0k − x k r Ω 0j , and (3.23), we get The estimate | Ω 0k (z)| z then directly ensues from (3.25), (3.27) and (3.29).
3.6. Decomposition of ∂ v . In this subsection, we state the decompositions and estimates that will allow us to deal with error terms of the form ∂ x i φ∂ v i ψ which appear in the commuted Vlasov equation (see Section 5), where φ is a function on M and ψ is a function on P. We start by introducing the notation The v derivatives are not part of the commutation vector fields and will be transformed using so that, for ψ a sufficiently regular solution to the free relativistic massless transport equation w µ ∂ µ ψ = 0, |∇ v ψ| essentially behaves as (t + r)|∇ t,x ψ|. In the following lemma, we prove that the radial component has a better behavior near the light cone.
Lemma 3.9. For the radial component of ∇ v the following estimates hold Let A denote a spherical frame field index. The angular part verifies the weaker estimates Proof. Since the assertion (3.31) follows by Lemma 3.8. For the first inequality of (3.32), recall that the vector field e A can be written as The second inequality of (3.32) is obtained by applying the last estimate to ψ = z and using again Lemma 3.8.
Similar to the case of the wave equation, we can then deduce that Lψ enjoys improved decay near the light cone. More precisely, This can be obtained by combining the previous Lemma with the relation

3.7.
The energy norms. We define here the energy norms both for the distribution function f and the metric perturbation h 1 . First, for any (a, b) ∈ R 2 , we introduce the weight function Then, define, for all sufficiently regular function ψ : For a, b ∈ R * + , an integer n ≥ 0 and a real number ℓ ≥ 2 3 n, we define the energies Remark 3.10. During the proof of Theorem 2.1, as we will take ℓ ≥ 1 8 and since 1 + |t − r| z according to Lemma 3.7, the energy norm E ℓ n [f ] will control Σt R 3 v Z I f dvdx for any |I| ≤ n.

Functional inequalities.
We end this section with some functional inequalities, starting with the following Hardy type inequality, which essentially follows from a similar one of [26].
Since ∇ ∂r V = ∇ ∂r W = 0, we have |∂ r (k V W )| = |∇ ∂r (k) V W | and the result follows from the definition of |∇k| VW .
The following technical result will be useful to prove boundedness for energy norms.
Proof. This follows from a integration by parts in the variable τ , Recall the decomposition (2.2), where χ is a smooth cutoff function such that χ = 0 on ] − ∞, 1 4 ] and χ = 1 on [ 1 2 , +∞[. It will be useful to control the derivatives of the cutt-off χ r t+1 which is the content of the next lemma.
Proposition 3.14. Let k be a sufficiently regular tensor field defined on Proof. It is sufficient to prove the proposition for scalar functions φ since we can apply the inequality to each cartesian component of k and then use that Recall the classical Klainerman-Sobolev inequality and that χ is a smooth cutoff function such that χ = 0 on ] − ∞, 1 4 ] and χ = 1 on gives, using Leibniz formula and Lemma 3.13, .
we obtain the result for the region considered. Consider now (t, x) ∈ [0, T [×R 3 such that |x| ≤ 1+t 4 and let us introduce τ − := (1 + |t − r| 2 ) 1 2 for regularity issues. Applying the classical Klainerman-Sobolev inequality (3.38) Note that , which can be obtained by following the proof of Lemma 3.13. In particular, we have r −1 (1 + t + r) −1 on the support of the two integrands on the right hand side of the previous inequality.
, so that t − r is bounded on the support of χ ′ (r − t) and χ ′ (t − r + 2), which implies the result.
Furthermore, we will need a slight improvement of the Klainerman-Sobolev inequality for massless Vlasov fields originally proved in [16].
Proof. As we do not have the inequality | Z I (z)| z at our disposal if |I| ≥ 2 and since ω b a is not C 3 class, one cannot apply a standard L 1 Klainerman-Sobolev inequality for velocity averages to z c f ω b a and derive the result. In fact, one just have to slightly modify one step of its proof. Remark . Hence, since | Z(z c )| z c according to Lemma 3.8, we obtain, applying Lemma 3.6, Following the proof of [8, Proposition 3.6], with f formally replaced by z c |v|f ω b a , and using (3.39) instead of Lemma 3.6, each time where this lemma is applied in [8, Proposition 3.6], we get the result.

Preliminary analysis for the study of the metric coefficients
In this section, we recall standard analytical properties of the metric coefficients in wave coordinates, independently of the Vlasov field. Most of the material of this section can be found in either [26] or [27]. In order to be self-contained, we present here not only the statements but also detailed proofs.
We fix, for all Sections 4-6, a sufficiently regular metric g and its decomposition as We assume that g is defined on [0, T [×R 3 , satisfies the wave gauge condition (3.3) and verifies the following regularities conditions. For an integer N ≥ 6 and 0 < ǫ ≤ 1 4 small enough, M ≤ √ ǫ and These conditions, which will be verified during the proof of Theorem 2.1 for a certain N ≥ 6 (see the bootstrap assumption (9.5) and the decay estimates of Propositions 10.1-10.2) and ǫ > 0, ensure that all the quantities considered in the next three sections are well-defined. In particular, the series of functions appearing below will be convergent in L 2 (Σ t ). Let us start by estimating pointwise the Schwarzschild part and its derivatives.
is the number of translations (respectively homogeneous vector fields) composing Z J 0 . By the Leibniz rule we have, By Lemma 3.13 and a straightforward computation, we have where P K P (t, r, x r ) is a certain polynomial in (t, r, x r ) which has degree K P in (t, r). Applying this to Z J 0 = Z J and using that 1 + t + r r on the support of h 0 as well as , we obtain the first estimate. For the second one, note that and apply (4.4)-(4.5) to Z J 0 = ∂ µ Z J for all µ ∈ 0, 3 .

4.1.
Difference between H and h. In this subsection, we study the difference between H µν := g µν − η µν and h µν := h αβ η αµ η βν . For this, let us define Using the expansion in Taylor series of the inverse matrix function, we then obtain The goal now is to compare H with h and H 1 with h 1 . In order to unify the treatment of these two cases, we consider ( Recall now, as the elements of K \ {S} are Killing vector fields and since S is a conformal Killing vector field of factor 2, that, when acting on the contravariant tensor η µν , As the lie derivative commute with contraction, this implies Iterating the previous arguments, we then obtain Moreover, using (4.6), we also obtain that (4.10) Similarly, one can prove that We then immediately obtain the following result. Here More precise inequalities will be required during the proof of Proposition 5.14 in the case where Z J contains at least one translation, i.e. J T ≥ 1. Since M T = J T in the sums on the right hand sides of (4.7)-(4.9) and that 1≤i≤n J T i = J T in the one of (4.10), we have

4.2.
Wave gauge condition. Using the wave gauge condition, one can estimate the bad derivative L of good components LT of the metric by good derivatives of the metric and cubic terms. We emphasize that the result also holds for L J Z (H) since, crucially, we are differentiating the metric geometrically.
Proposition 4.4. Let N ≥ 6 be such that (4.2) holds and assume that the wave gauge condition is satisfied. Then, for all |I| ≤ N , we have Remark 4.5. From the wave gauge condition, one can also derive It can be obtained by expressing (4.14) in terms of H instead of h and by following the rest of the upcoming proof. Note that a slightly weaker estimate could be obtained by combining Propositions 4.2 and 4.4.
Proof. Remark first that we only need to prove these inequalities for ∇ L L I Z (h) LT and . In order to lighten the notations, we will use O µν (|h| 2 ) in order to denote a tensor field of the form • For all |J| ≤ N , +∞ n=2 L J Z (P n (h)) and +∞ n=2 ∇L J Z (P n (h)) are absolutely convergent in L 2 (Σ t ) and we have This will be implied by the fact that g satisfies the condition (4.2). • The tensor field O µν (|h| 2 ) can be different from one line to another. Recall from (3.3) that the wave gauge condition implies Expanding the determinant of g (the first order term is the trace), we have where P(|h| 2 ) is a polynomial in the variables (h αβ ) 0≤α,β≤3 of degree at most 4 and of valuation at least 2. Hence, using H µν = −h µν + O µν (|h| 2 ) and the expansion in Taylor series of the square root function, we get 11 (4.14) Now, observe by a straightforward calculation that for a general tensor field F µν , we have for all Z ∈ K and since the Lie derivative commute with contractions, The identities (4.14), (4.15) and (4.16) yield, by an easy induction, to For a vector field U and a tensor field F µν , there holds the formula Applying this identity to U = T ∈ T , F = L I Z (h) and then F = tr(L I Z h)η, one has, since η LT = 0, Combining (4.17) with (4.13), (4.19) and (4.20), we obtain The first estimate (4.11) then follows from We now turn to the second one. Note first that As h = h 0 + h 1 and δ µν − η µν = 2δ 0µ δ 0ν , the condition (4.14) leads to As the support of χ ′ is included in [ 1 4 , 1 2 ], we obtain, since Z J is a combination of translations and homogeneous vector fields 12 , Using (4.15) and (4.16), we then get for all |J| ≤ N and ν ∈ 0, 3 , Since (4.19) and (4.20) also hold if h is replaced by h 1 , the inequality (4.12) ensues from (4.13) and (4.22).

4.3.
Commutation formula for the Einstein equations. In this section, we compute the source terms of the wave equation satisfied by the cartesian components of L J Z (h 1 ). In order to do it in a geometric way, we define, for any sufficiently regular (0, 2)-tensor field k, the (0, 2)-tensor field g (k) whose components in wave coordinates satisfy since ∇ is the covariant differentiation of Minkowski spacetime whose Christoffel symbols vanish in the coordinates system (t, x). Our goal now is to compute, for any Z J ∈ K |J| , g (L J Z h 1 ). The first step consist in determining the commutator g (L J Z h 1 ) − L J Z ( g h 1 ) and then we will describe L J Z ( g h 1 ). We start by the following technical result. Lemma 4.6. Let K be a (2, 0)-tensor field and k a (0, 2)-tensor field, both sufficiently regular. Then, for all Z ∈ K, we have Proof. We will use here that K αβ ∇ α ∇ β k is obtained by contracting K with the (0, 4)tensor field ∇∇k. Since the Lie derivative commute with contraction, we have for any 0 ≤ µ, ν ≤ 3 and for all Z ∈ K, It then remains to apply Lemma 3.4, which gives ( We are now able to compute the commutator. 12 We refer to the proof of Lemma 3.13 for a more detailed estimate of a similar quantity.
For all multi-index |I| ≤ N , there exist integers C I K , C I J,K ∈ Z such that Proof. Let Z ∈ K and recall that g (h 1 ) = g αβ ∇ β ∇ α h 1 . Then, applying Lemma 4.6, we get It only remains to use For the higher order commutation formula, we proceed by induction on |I| (note that the result is straightforward if |I| = 0). Let n ∈ N * and assume that the result holds for all multi-indices |I 0 | = n. We then consider a multi-index I of length n + 1 and we introduce Z ∈ K and |I 0 | = n such that Z I = ZZ I 0 . Then, According to the first order commutation formula applied to L I 0 Z h 1 , All the terms on the right hand side of this equality have the required form since 1 ≤ |I 0 | < |I|. Using the induction hypothesis, we can write L Z g L I 0 Z h 1 − L I 0 Z g h 1 as linear combination of terms of the form It remains to apply Lemma 4.6 in order to deal with the first ones and the first order commutation formula for the last ones (note that |J| + |K| + 1 ≤ |I 0 | + 1 = |I| and |K| + 1 < |I|).
We now focus on L J Z g h 1 .
Lemma 4.8. Let k and q be two sufficiently regular (0, 2)-tensor fields. Then, for all Z ∈ K, Iterating these relations, we obtain that for all |I| ≤ N , there exist integers C I J,K such that Proof. This directly follows from the definition of P (∇k, ∇q) and Q(∇k, ∇q) (3.5)-(3.6) as well as Lemma 3.5.
We then deduce the commutation formula for the Einstein equations (3.4a). such that, for any (µ, ν) ∈ 0, 3 2 , The derivatives of T [f ] and g h 0 will be computed in Section 6 and Proposition 11.2. For the cubic terms, we have under the assumption (4.2), Proof. The commutation formula for the Einstein equations (3.4a) follows from an induction on |I| relying on Corollary 4.7 and Lemma 4.8. For the estimate on the cubic terms, we obtain from (3.7) and the definition of the Lie derivative (3.8) that L I Z (G(h)(∇h, ∇h)) µν can be bounded by a linear combination of terms of the form where all the multi-indices are in 0, 3 and |J 0 | + |J 1 | + |J 2 | + |J 3 | ≤ |I|. Note now, using (3.9) and Lemma 3.4 that Finally, without loss of generality, we can assume that |J 0 | ≤ N − 3, so that, using Proposition 4.2 and the assumption (4.2), Z J 0 H α 0 β 0 1. This concludes the proof.

Commutation of the Vlasov equation
The purpose of this section is to compute the commutator [T g , Z I ], for Z I ∈ P |I| 0 . The commutation formula obtained here is more geometric than the one used by [16]. In the spirit of [8] for the Vlasov-Maxwell system (see in particular Subsection 2.5), we express the error terms using Lie derivatives of the metric instead of derivatives of its Cartesian components. We recall the following notations and we consider for all this section a sufficiently regular symmetric tensor field H µν and a sufficiently regular function ψ : We define the vertical parts S w and Z w , for Z ∈ P a Killing, respectively conformal Killing, vector field, by Recall also that, in order to simplify the presentation of the commutation formula, we use the following convention. For any Z ∈ P 0 , if Z = S, then we denote by Z the Killing vector field which has Z as its complete lift and if Z = S, then we set Z = S. Finally, we extend the Kronecker symbol to vector fields (X, Y ), i.e. δ Y X = 1 if X = Y and δ Y X = 0 otherwise.

Geometric notations.
In order to clearly identify the structure of the error terms in the commuted equations, let us rewrite the two parts composing the operator T g . For this, we will denote the differential in the spacetime variables (t, x) of ψ by dψ and we recall that ∇H denotes the covariant derivatives of H with respect to the Minkowski metric. We then have With these notations, Note that the transport operator can then be rewritten as and where T η = |v|∂ t + v i ∂ i ψ = w µ ∂ µ is the massless relativistic transport operator with respect to the Minkowski metric. Let us mention that the quantity (5.3) will appear as an error term in the commutator [T g , Ω 0k ]. We now prove a technical lemma which contains useful identities.
Lemma 5.1. Let θ = θ µ dx µ and θ = θ µ dx µ be two 1-forms and Z ∈ P 0 . Then, Proof. As the Cartesian components of w do not depend on (t, and then that . In order to compute (5.7) and (5.8), let us introduce

5.2.
Commutation formula for T g . We start by deriving a commutation formula for the first part T g of the transport operator. To this end, we first decompose it as The following lemma is a prerequisite for Lemma 5.3.
Applying the identity (5.6) of Lemma 5.1, we get which leads in particular to and then concludes the first part of the proof. The second formula follows from We then derive the commutation formula for the operator T g .
If Z I ∈ P |I| 0 , there exists integers C I Q , C I J,K and C I µ,J 1 ,J 2 ,K such that where the multi-indices J, J 1 , J 2 and K in the last two sums satisfy one of the following two conditions, Combining the first order commutation formula with the identity (5.20), written below, one can check that Z K and Z Q (respectively Z J , Z J 2 and Z J 1 ) is built by at most |I| − 1 (respectively at most |J|, at most |J 2 | and at most |J 1 |) of the vector fields composing Z I , so that K P ≤ I P and Q P ≤ I P . If K P = I P , this means that there is at least one translation in Z I which is part of Z J and either Z J 2 or Z J 1 , i.e. J T ≥ 1 and Proof. Let Z ∈ P 0 and recall from Subsection 3.5 that Applying the first equality of Lemma 5.2 to H = H and the second one to H = g −1 and µ = 0, we get The first order commutation formula directly follows from (5.17), (5.18) and (5.19). The higher order formula can be proved similarly by performing an induction on |I|, using (5.20) [ and applying the first equality (respectively the second equality) of Lemma 5.2 to Z K ψ and H = L J Z (H) (respectively H = L J 2 Z (g −1 ) ), for well-chosen multi-indices J, J 2 and K.
Remark 5.5. Expressing the error terms in the commutation formula using v instead of w, we find, since

5.3.
Commutation formula for the transport operator. In view of Lemma 5.3 it remains to study the action of Z I on the term The following identities will then be useful in order to determine [T g , Z I ].
We are now able to compute the first order commutation formula. In fact we will state it in two different ways. The second one has the advantage of being more concise whereas the first one will be more adapted to the problem studied in this paper and for the purpose of deriving the higher order formula.
Alternatively, expressing the error terms using v instead of w, we get Proof. The first commutation formula follows from Lemma 5.3 and Lemma 5.6 applied to H = H and (µ, ν) = (0, 0). The second formula can be obtained from the first one using that v = w + ∆vdt and Remark 5.8. Even if the second commutation formula might seem to be more convenient, we will work with the first one for two reasons.
• The second and higher order formulas are not more concise when expressed in terms of v instead of w. • Working with w instead of v is more adapted to our method since no inequality analogous to |w L | w 0 z 2 (1+t+r) 2 holds for the component v L . Indeed, according to Lemma 5.12 proved below and | / w| |v||w L | (see Lemma 3.7), we have, if g satisfies (4.2) and for ǫ small enough, Although we will have, during the proof of Theorem 2.1, |w L ||H|+ |v||w L ||H| LT |v| z 2 (1+t+r) 2 , the term |v||H LL | will not behave sufficiently well near the light cone. Because of the Schwarzschild part, |H LL | cannot decay faster than (1 + t + r) −1 and no decay can be extracted from the weight z if t ≈ r without a good component of the flat velocity vector w L or / w.
Due to the new error terms generated by the Lorentz boosts, the following additional identities are required in order to compute the higher order commutation formula. Lemma 5.9. Let Z ∈ P 0 , (λ, ν) ∈ 0, 3 and q ∈ 1, 3 . Then, Proof. Note first that Then use the identities (5.6) and (5. We are now ready to describe the error terms of the higher order commutator [T g , Z I ] in full detail. Proposition 5.10. Let Z I ∈ P |I| 0 . Then, [T g , Z I ](ψ) can be written as a linear combination with polynomial coefficients in w ξ w 0 , 0 ≤ ξ ≤ 3, of the following terms, Moreover K, J, Q and M 1 satisfy the following condition (1) either K P < I P , (2) or K P = I P and then J T ≥ 1, The conditions on the multi-indices are easy to check when |I| = 1 (see Proposition 5.7). In that case there holds |K| = K P = 0. So, if Z I = Z is a homogeneous vector field, we have K P < I P = 1. Otherwise, Z I is a translation ∂ x µ and each source term contains either the factor L ∂ x µ (H) or ∂ x µ (∆v). Moreover, K P < I P always holds for the terms of the form (5.27) since they do not appear when Z I = ∂ x µ . One can check during the induction, and more precisely when we apply Lemmas 5.6 and 5.9, that these conditions hold for all I (the general principle is explained in Remark 5.4).
Remark 5.11. As mentioned in Subsection 2.4.3, we would not be able to close the energy estimates on the Vlasov field without taking advantage on the conditions on K P and I P given in Proposition 5.10.
We also point out that the condition K P < I P for the terms (5.27) is of fundamental importance. We would not be able to handle such terms if the case K P = I P was possible, even if we had at the same time J T ≥ 1.

5.4.
Null structure of the error terms in the commuted Vlasov equation. The aim of this subsection is to describe the null structure of the terms given by Proposition 5.10. We start by estimating Z M (∆v), which will be useful in order to deal with (5.28)-(5.32). with |M | ≤ N and assume that the metric g satisfies the wave gauge condition and (4.2). Then, if ǫ is sufficiently small, we have Proof. According to Proposition 4.2 and (4.2), we have , which implies, since w 0 = −|v| and if ǫ is sufficiently small, Consequently, As |H| √ ǫ, we obtain, if ǫ is sufficiently small, that |∆v| ≤ 2 |H(w,w)| |v| . Now, recall from Lemma 3.7 that w A w A |v||w L |, which implies and the result holds for |M | = 0. The next step consists in proving an inequality which will allow us to prove the result by induction in |M |. The starting point is the decomposition dx k and (5.6), we get It then follows that Iterating the process, one can prove that, for all Z M ∈ P .
Consider now N 0 ≤ N − 1 and suppose that (5.33) holds for all |I| ≤ N 0 . Then, let M be a multi-index satisfying |M | = N 0 + 1.
Following the computations made in (5.36), we then get In order to bound the second sum in the right hand side of (5.37), start by noticing that, Now, by the induction hypothesis, The claim then follows from (5.37), (5.38), the last two inequalities and which is a direct consequence of the induction hypothesis.
In the next lemma, we deal with the remaining error terms given by (5.25), (5.26) and (5.27) by expanding them with respect to the null frame (L, L, e 1 , e 2 ).
Lemma 5.13. The following estimates hold, Proof. The first inequality follows from and from Lemma 3.7 as well as (3.35), which give Remark now that for a symmetric tensor G µν , Consequently, using again that |w A | |v||w L |, we get Recall from Lemma 3.9 that The last two estimates then result from (5.39), (5.40), (5.41) and

5.5.
Final classification of the error terms. In this section, we list of all the error terms that appear in the commuted equations in such a way that we will able to easily estimate them when we try to improve all the bootstrap assumptions on the energy norms of the Vlasov field.
Proposition 5.14. Let N ≥ 6 be such that the metric g satisfies (4.2), assume that the wave gauge condition holds and consider Z I ∈ P |I| 0 with |I| ≤ N . Then, [T g , Z I ](ψ) can be bounded by a linear combination of terms taken in the following families. The terms arising from the source terms The terms arising from the Schwarzschild part, where, Z ∈ P 0 , • |Q| + |J| + |K| ≤ |I|, |K| ≤ |I| − 1, K P ≤ I P . The quadratic terms, where, Z ∈ P 0 , • |J| + |K| ≤ |I|, |K| ≤ |I| − 1. • K and J satisfy one of the following conditions.
(1) Either K P < I P , (2) or K P = I P and J T ≥ 1.
(1) Either K P < I P , (2) or K P = I P and M T + J T ≥ 1.
Remark 5.15. To clarify the analysis, we have denoted by S or E, the error terms that contains factors of Z Z K ψ , and by S or E, error terms containing ∇ Z K ψ , so that we know that the last derivative hitting ψ is a translation.
Proof. Since g verifies (4.2) and in view of Proposition 4.2, we will use throughout this proof that Consider a multi-index I such that |I| ≤ N . In order to clarify the analysis, let us introduce a notation. Fix q ∈ 4, 11 and multi-indices (J, K) satisfying the conditions presented in the proposition which are associated to E J,K I,q . Then, for a sufficiently regular tensor field k, denote by E J,K I,q [k] the quantity corresponding to E J,K I,q , but where h 1 is replaced by k. For instance, We [h], where p ∈ 14, 17 and (J 0 , K), (M 0 , J 0 , K) as well as (Q 0 , M 0 , J 0 , K) satisfy the conditions presented in the proposition. This follows from Remark 4.3, so that, for instance, [h].
Similar relations can be obtained, using also ( [H]. (2) For all n ∈ 1, 3 and q ∈ 4, 11 , we have This ensues from the decomposition h = h 1 + h 0 and Proposition 4.1, which gives that, for all |J|, Then, the last term is bounded by E J,K I,1 if K P < I P . Otherwise K P = I P and M T + J T ≥ 1, so that it can be bounded by E J,K I,3 if M T ≥ 1 and by E J,K I,1 if J T ≥ 1. The remainder of the proof then consists in bounding the terms written in Proposition 5.10 by (5.42) and those of (5.51)-(5.68), with h 1 replaced by H. For that purpose, we will use several times Lemmas 5.12 and 5.13. Until the end of this section, each time that we will refer to one of the terms (5.51)-(5.68), h 1 has to be replaced by H. • with |M | + |Q| + |J| + |K| ≤ |I|, |K| ≤ |I| − 1 and K P < I P or K P = I P and which comes from (5.41), we finally get quartic terms of the form (5.68) and, using (5.69), cubic terms (5.63). • Finally, since for two functions φ and ψ, there holds we can bound, using (5.41), the terms (5.29) and (5.31) by with |M 1 | + |J| + |K| ≤ |I|, |K| ≤ |I| − 1 and K P < I P or K P = I P and which follows from Lemma 5.12, leads to terms of the form (5.63) and (5.65)-(5.68).
It will be convenient to introduce the following notations.
Definition 5. 16. Given one of the error terms E J,K I,i , i ∈ 4, 11 , listed in Proposition 5.14, we define A J,K I,i as the quantity which contains everything of E J,K I,i but the ψ-part |∇ Z K ψ|. We define similarly, for n ∈ 1, 3 and p ∈ 14, 17 , , so that

Commutation of the Vlasov energy momentum tensor
To evaluate the commuted Einstein equations (see Proposition 4.9), we will require the null components of the tensor field L I Z (T [f ]). In order to simplify the presentation of the following results as well as their proofs, we denote by T [ψ] the energy-momentum tensor of the Vlasov field in the flat case, i.e.
This field is considered in the following. Proof. The result for the Killing vector fields Z ∈ P holds in a more general setting. More precisely, if X is Killing for a metric h and T [ψ] is the energy-momentum tensor of a Vlasov field ψ for the metric h, then L X T [ψ] = T [ Xψ], with X the complete lift of X, as can easily be verified by choosing a local coordinate system such that X coincides with one of the coordinate derivatives. For the scaling vector field 13 , S = x µ ∂ µ , we have We now turn on the real energy momentum tensor T [ψ]. 13 The types of formula can be in fact generalized to any conformal Killing fields on a general Lorentzian manifold.
Proposition 6.2. Let I be a multi-index and Z I ∈ K |I| . Then, there exist integers C I J,K , C I,λ J,K,M ;µν and C I J,K,L,M ;µν such that Proof. The formula is satisfied for |I| = 0 since w 0 = |v| and The result for arbitrary multi-indices I follows by induction, applying several times Lemmas 3.6 and 6.1.
Recall that the metric g satisfies the decomposition (4.1) and the condition (4.2).
Proposition 6.3. Let N ≥ 6 and g be a metric such that (4.2) holds. Then, for all Z I ∈ K |I| such that |I| ≤ N and V, W ∈ U , we have, if ǫ small enough, Proof. Note first that according to Proposition 4.2 and the assumptions (4.2), Hence, using Lemma 5.12, we have Suppose that holds. Then, from Proposition 6.2 and (6.3)-(6.4), it holds The result then follow from which holds for any |J| ≤ N and follows from the decomposition h = h 0 + h 1 and Proposition 4.1. It then only remains to prove (6.4). For this, note first that, using v = w + ∆vdt, g −1 = η −1 + H, (6.3) and (6.2), Similarly, using that det(g −1 ) is a polynomial of degree 4 in g µν , 0 ≤ µ, ν ≤ 3, we get Using |H| √ ǫ, |∆v| √ ǫ, v = w +∆vdt, (6.3), and that the determinant is a multilinear mapping, we obtain, for ǫ small enough, The inequality (6.4) then follows from the Leibniz rule, | Z Q (w 0 )| ≤ C Q |v| and the last four estimates.
Remark 6.4. Note that a better estimate could be obtained for the good components of L I Z (T [f ]) in Propositions 6.2 and 6.3 but the result stated in this section will be sufficient in order to close the energy estimates.

Energy estimates for the wave equation
The aim of this section is to prove energy inequalities for solutions to wave equations in a curved background whose metric g is close and converges to the Minkowski metric η. These results can be found in Section 6 of [26] and we give here, for completness, an slightly different proof. More precisely, the goal is to control, for some (a, b) ∈ R 2 + and a sufficiently regular function φ, energy norms will allow us to take advantage of the decay in t − r. Without an a priori good estimate on it, we would merely obtain that Note however that the bulk integral provides only a control on the derivatives tangential to the light cone, i.e. L and / ∇, and constitutes an important tool in order to exploit the null structure of the massless Einstein-Vlasov system. The problem when a = 0 or b = 0 is that the energy estimate derived below (see Proposition 7.5) will not allow us to control K. Moreover, if a > 0, the norm r≤t |∇ t,x φ| 2 ω b a dx is strictly weaker than r≤t |∇ t,x φ| 2 dx, which explains why we introduce E a,b [φ].
We introduce the energy normE a,b [φ] in order to avoid a strong growth at the top order which would force us to assume more decay on the initial data in order to close the energy estimates.
We fix, for the remaining of this section, T > 0 as well as a function φ and a metric g, both defined on [0, T [×R 3 and sufficiently regular. We also introduce H := g −1 − η −1 . In order to derive energy inequalities, we introduce the (1, 1)-tensor field The tensor field T [φ] is the energy momentum tensor of φ, written as a (1, 1) tensor. However, we point out that since we lower indices with respect to the Minkowski metric, T [φ] µν = ∂ µ φ∂ ν φ − 1 2 g µν g αβ ∂ α φ∂ β φ. The (1, 1) tensor field T [φ] appears to be well adapted in order to prove energy estimates for the norms that we are interested in.
Let us now compute the divergence of T [φ]. For this, it will be convenient to use the notation Remark 7.4. In general, T µν [φ] is not symmetric.
We are now ready to provide an alternative proof of Proposition 6.2 of [26].
Proposition 7.5. Let a, b ∈ R * + , C H > 0 and suppose that H satisfies Then, there exists a constant C := C 0 1+a+b min (1,a,b) , where C 0 > 0 is an absolute constant, such that, if ǫ is sufficiently small 14 Finally, there also holds Proof. In order to lighten the proof, we will not keep track of the constant C H , which appears merely when √ ǫ does. The (euclidian) divergence theorem yields As |H| √ ǫ, we have, if ǫ is sufficiently small enough, 14 One can check that ǫ needs to satisfy a condition of the form C1CH √ ǫ(1 + a + b) ≤ 1 4 min (1, a, b), for a certain constant C1 > 0.
The first inequality (7.1) then follows, if ǫ is sufficiently small 15 , from the second equality of Lemma 7.3 as well as t 0 Στ In order to prove (7.6), start by noticing that This, together with
We now turn on the second inequality (7.2), which can be obtained by taking the sum of (7.1) and 16 To prove this estimate, apply the euclidian divergence theorem to T [φ] µ 0 and follow the proof of (7.1). The identity (7.4) does not depend of (a, b) and (7.5)-(7.6) are trivial for (a, b) = (0, 0) as ω 0 0 = 0. It then remains to bound sufficiently well the left hand side of (7.7) when (a, b) = (0, 0). For this note that (7.8), (7.9) and (7.10) still hold in that context and that Finally, (7.3) can be proved similarly as (7.1) by applying the divergence theorem to 1+t+r (see Lemma 7.3). Apart from the fact that each integral contains an extra |1 + t + r| −1 (or |1 + τ + r| −1 ) weight, the only significant difference is that we need to control In view of sign considerations and since |H| √ ǫ, we can bound it by which concludes the proof.

L 1 -Energy estimates for Vlasov fields
Let ψ be a sufficiently regular function defined on the co-mass shell P and recall the Vlasov L 1 -energy In this section, we prove the following L 1 -energy estimate for Vlasov fields.

Proposition 8.1. Assume the bounds
For any parameters a, b > 0 and 0 ≤ t 1 ≤ t 2 < ∞ and any sufficiently regular function ψ : where C and C are two constants such that C depends only on (a, b).
The boundary terms at t = t i are given by Thus, using (8.4) and the assumptions on H, Consider now the last term on the RHS of (8.2), for which we have which we rewrite as In view of the bounds on H, it follows that, so that using (8.4), we have It follows that the contribution of the last term on the RHS of (8.2), for some constant C > 0, and, using (8.3)-(8.5), that The LHS of this last inequality will provide the spacetime term of E a,b [ψ](t 2 ) when we sum all the terms at the end of the analysis. Note that it will arise with the same sign as the boundary term at t = t 2 . Finally, we consider the contribution of the terms and we use Lemma 5.12 in order to get Using the assumptions on H, we derive, since |v A v B | |v||w L | by Lemma 3.7, where we note that the contribution of the first term on the RHS can be absorbed if ǫ is small enough into the spacetime positive term containing |w L | obtained above, while the contribution of the second term can be simply estimated in terms of the energy.

Bootstrap assumptions
We consider the following bootstrap assumptions on certain energy norms which has been defined in Subsection 3.7. Let N ≥ 13, ℓ = 2 3 N + 6 and consider the parameters 0 < 20δ < γ < 1 20 . We have • bootstrap assumptions for the Vlasov field: For all t ∈ [0, T [, where C f , C, C T U and C LL are constants larger than 1 which will be fixed during the proof in Section 12. As is usual for this type of proof, the above bootstrap assumptions are satisfied with strict inequality for t = 0 by our assumptions on the initial data and provided that C f , C, C T U and C LL are large enough. By standard well-posedness theory, it follows that they are satisfied on some non-empty interval of time [0, T [, with T > 0. Theorem 2.1 then holds provided we can improve each of the above bootstrap assumptions.
The growth on the bootstrap assumptions (9.1), (9.2) and (9.8) are independant from all the other ones and could be choosen to be of the form (1 + t) η , with η arbitrary small.
The following result will be useful in order to improve the boostrap assumptions (9.6)-(9.8). The rough idea is that the L 2 -norm of |∇L J Z (h 1 )(V, W )| and |∇ L J Z h 1 (V, W ) | are equivalent.
Lemma 9.2. There exists a constant C > 0 independant of C, C T U and C LL such that, for all t ∈ [0, T [, Proof. For the purpose of keeping track of certain quantities, all the constants hidden in will be independant of C, C T U and C LL . This convention will only hold during this proof. In order to lighten the notations, we introduce k J := L J Z (h 1 ) for any |J| ≤ N . Then, observe that according to the triangle inequality, the lemma would follow if we could prove the first inequality (respectively the last two inequalities) with N − 1 (respectively N ) replaced by 0 and h 1 by k J for any |J| ≤ N − 1 (respectively |J| ≤ N ).
We start by an intermediary result. Let us fix (V, W) ∈ {U , T , L} 2 , 0 ≤ a ≤ 1 + 2γ and 0 ≤ b ≤ 1 + γ. Since Note that since the domain of integration of the four integrals on the right hand side of the previous inequality are located far from the light cone, we do not keep track of 17 V and W. Our goal now is to bound them sufficiently well for well choosen values of |J| and (a, b) in order to obtain Cǫ. (9.14) 17 It is only near the light cone that certain null components of the metric enjoy improved decay estimates.
For the purpose of controlling the four integrals on the right hand side of (9.11), we will use many times the inequality 1 + τ + r 1 + |τ − r| which holds on their domain of integration. We start by dealing with the case |J| ≤ N − 1 and (a, b) = (2γ, 1 + γ).
Applying the Hardy inequality of Lemma 3.11 and making similar computations, one gets We now assume that |J| ≤ N and we introduce η ∈ {0, γ} in order to unify the treatment of the remaining two cases. We have, Applying the Hardy inequality of Lemma 3.11, one obtains Now recall from the bootstrap assumptions (9.4) and (9.5) that Using also that 2δ < γ, we can deduce (9.12)-(9.14) from the last estimates. We now turn on the second part of the proof. Note that According to the Hardy type inequality of Lemma 3.11 and the bootstrap assumptions (9.5) and (9.7), we have 18 , since 2δ < γ, The third inequality of the Lemma then ensues from (9.14), (9.15) and these last two estimates. By similar considerations, one can obtain for |J| ≤ N − 1, and, for |J| ≤ N , 18 Note that we could avoid the use of the bootstrap assumption (9.7) by taking advantage of the wave gauge condition. The consequence is that the right hand side of the third inequality of Lemma 9.2 could be independant of CT U .
All these integrals will be estimated using the Hardy inequality of Lemma 3.11. For those of (9.16), we have Using the bootstrap assumptions (9.4) and 2δ < γ, we have The first inequality of the Lemma follows from these last three estimates, (9.12) and (9.16).
For the integrals on the right hand side of (9.17), one has, according to the bootstrap assumption (9.5), The second inequality of the Lemma then ensues from the last two estimates, (9.13) and (9.17).

Pointwise decay estimates on the metric
We prove here pointwise decay estimates on h 1 and its (lower order) derivatives using the bootstrap assumptions (9.4) and (9.6). The Schwarzschild part h 0 can always be estimated pointwise using its explicit form. This will then allow us to obtain asymptotic properties on h = h 1 + h 0 .
Proof. The first inequality directly follows from the bootstrap assumption (9.4) and the Klainerman-Sobolev inequality of Proposition 3.14, applied with a = 0 and b = 1 + 2γ. Let |J| ≤ N − 3, θ ∈ S 2 , (µ, ν) ∈ 0, 3 and so that L J Z (h 1 )(t, rθ) = ϕ(t + r, t − r). We start by considering the exterior of the light cone, i.e. we fix (t, r) ∈ [0, T [×R * + such that r ≥ t. Hence, We can now treat the remaining region and we then fix (t, r) ∈ [0, T [×R * + such that r ≤ t. We have For the third estimate, we use the inequality (3.11) of Lemma 3.3 and the estimate (10.2).
In order to obtain the decay rate of L J Z (h), for |J| ≤ N − 3, it remains to study h 0 and its derivatives. The following result is a direct consequence of Proposition 4.1 and M ≤ √ ǫ.
Proposition 10.2. For all Z J ∈ K |I| , there exists C J > 0 such that for all (t, x) ∈ R + ×R 3 , Remark 10.3. In the interior of the light cone, the behaviour of L J Z (h) is clearly given by L J Z (h 1 ). In the exterior region, note that L J Z (h 0 ) has a weaker decay rate than L J Z (h 1 ) when r > 2t but a stronger one when t ∼ r.
We can improve the decay estimates satisfied by certain null components of h 1 through the wave gauge condition. According to Proposition 4.4 as well as the pointwise decay estimates given by Propositions 10.1 and 10.2 (recall that h = h 0 + h 1 ), we obtain the following results.
ǫ (1 + t + r) 6 (10.5) Remark 10.5. This inequality will be used several times in this article. Apart from its application during the proof of Propositions 12.8 and 13.4 below, we will always bound the term ∇L J Z h 1 2 T U by ∇L J Z h 1 2 .
It then remains to bound the right hand side of the previous inequality. Let us fix µ ∈ 0, 3 and (T, U ) ∈ T × U . Using Lemma 3.13 we get, for any |I| ≤ 2, .
Using the notation [Z 1 Z 2 , X] in order to denote [Z 1 , [Z 2 , X]] for any vector fields Z 1 , Z 2 and X, we can bound the right hand side of the previous inequality by .
• Following the proof of (3.17) and using one can prove that for all r ≥ 1+t 2 and |L| ≤ 2, where |d X | 1+|t−r| 1+t+r and |b W | + |b Y | 1 since 1 + t + r r on this region.
We then deduce, since 1+|t−r| The pointwise decay estimate (10.6) then follows from the bootstrap assumptions (9.4) and (9.6) as well as 2δ < γ. Now consider the LT components and assume that |J| ≤ N − 4. The first estimate can be obtained from the wave gauge condition (10.5) and the three inequalities of Proposition 10.1. For the second one, fix θ ∈ S 2 and consider, for T ∈ T , the function Let now (t, r) ∈ [0, T [×R * + such that r ≥ t. Using the estimate (10.7) and the good decay properties of the initial data, we obtain On the other hand, if r ≤ t, we have so that we will be able to apply the energy estimates of Propositions 7.5 and 8.1 for wellchosen parameters a and b.
The estimate |∇H| LL √ ǫ 1+|t−r| (1+t+r) 2 , which can be obtained in a similar way, will also be useful.
When h 1 is differentiated by at least one translation, we can improve the pointwise decay estimates given by Propositions 10.1 and 10.6. Note that certain of the following decay rates could be improved, in particular in the exterior of the light cone.

Bounds on the source terms of the Einstein equations
The aim of this subsection is to bound the source terms of the commuted Einstein equations which are given in Section 4.3. We will control them sufficiently well in order to close the energy estimates but more decay in t − r could be proved for certain terms. We start by the semi linear terms L I Z (F (h)(∇h, ∇h)) µν = L I Z (P (∇h, ∇h)) µν +L I Z (Q(∇h, ∇h)) µν +L I Z (G(h)(∇h, ∇h)) µν .
Proposition 11.1. Let I be a multi-index with |I| ≤ N . Then Proof. Let |I| ≤ N and recall from Lemma 4.8 that there exists integers C I J,K such that Moreover, according to Proposition 4.9 and the split the split h = h 0 + h 1 , We start by dealing with the cubic terms and we define for j, k, q ∈ {0, 1} and multi-indices J, K, Q such that |J| + |K| + |Q| ≤ |I|, Using the pointwise decay estimates given by proposition 10.2 on h 0 and its derivatives, we have (1 + t + r) 5 Finally, using also the pointwise decay estimates given by Proposition 10.1 on h 1 and its derivatives (at most one of the multi-indices J, K and Q has a length larger than N − 3), it follows I 0,1,1 J,K,Q + I 1,0,1 J,K,Q + I 1,1,0 The inequalities (11.1)-(11.3) provide a sufficiently good bound on the cubic terms for the purpose of proving the three estimates of the Proposition. Consider now the semi-linear terms Q and P . Start by decomposing h into h 0 + h 1 so that, using the pointwise decay estimates on h 0 given in Proposition 10.2, we get for any null components (V, W ) ∈ U 2 , It then remains us to study the last term of the previous two inequalities for (V, W ) ∈ U U (respectively (V, W ) ∈ T U and (V, W ) = (L, L)) in order to derive the first (respectively the second and the third) estimate of the Proposition. For the quadratic terms P , recall from Lemma 3.1 that if V = W = L, the null condition is not satisfied. More precisely, Hence, using the pointwise decay estimates given by Propositions 10.1, 10.2 and 10.6 as well as the wave gauge condition (4.12), we find that for any null components (V, W ) ∈ U 2 , Since (1 + |u|) γ ≤ (1 + |u|) 1 2 −γ and according to (11.1)-(11.3), this bound is sufficient in order to prove the first estimate of the proposition. Now we deal with the T U components of P and the U U components of Q together. According to Lemma 3.1 and the pointwise decay estimates of Proposition 10.1, we have for any (T, U ) ∈ T × U and (V, W ) ∈ U 2 , Note that this inequality need to be improved in order to prove the third estimate of the Proposition, i.e. for the case T = U = V = W = L, but is sufficient for the first two estimates. Finally, applying again Proposition 10.1 and Lemma 3.1, we obtain This implies the last estimate of the Proposition and concludes the proof.
Next we consider the Schwarzschild part h 0 .
Proposition 11.2. Let I be a multi-index such that |I| ≤ N and (µ, ν) ∈ 0, 3 2 . Then, Proof. Recall from Subsection 4.3 the definition of the tensor field g h 0 and start by decomposing g in η + H σθ ∇ σ ∇ θ . Then, as η 1 r = 0, we have, for all 0 ≤ µ, ν ≤ 3, According to (3.9), there holds 0≤µ,ν≤3 Fix then |Q| ≤ |I|. One can easily check, by similar computations as those made in the proof of Proposition 4.1 and in view of the support of χ ′ , that 2 } . Similarly, since 1 + t + r r on the support of χ( r t+1 ) and using (3.9), we have We now estimate the error terms arising from the commutator g L J Z h 1 − L J Z g h 1 .
Proposition 11.3. Let n ≤ N and J, K be multi-indices such that |J| + |K| ≤ n and |K| ≤ n − 1. For V, W ∈ {U , T , L}, there holds For the LL component, we have the improved estimate Proof. Start by noticing that for V, W ∈ {U , T , L}, Applying Lemma 3.3 and using that [Z, ∂ λ ] ∈ {0}∪{±∂ ν / 0 ≤ ν ≤ 3} as well as L ∂ν = ∇ ∂ν yields Applying Proposition 4.2, which makes the transition from H to h precise, and then using the split h = h 1 + h 0 as well as the pointwise decay estimates given by Propositions 10.2, for the Schwarzschild part h 0 , and 10.1, for h 1 , one obtains We then deduce that Note that one factor of each of the quadratic terms in h 1 can be estimated pointwise since N ≥ n ≥ 12. Hence, using the decay estimates given by Propositions 10.1 and 10.6, we obtain the following bound In order to estimate the first term on the right hand side of the previous inequality, we use the pointwise decay estimates of Propositions 10.1 and 10.6 which provides The asserted bounds now follow (note that we use δ ≤ 1 2 and that we do not keep all the decay given by the last estimates.) Finally we bound the error terms coming from the commutation of g with the contraction with the frame fields T U or LL and the commutation of g with the multiplication by the characteristic function χ r 1+t .
Lemma 11.4. Let k µν be a (2, 0) tensor field and (T, U ) ∈ T × U . Then Proof. We will use in the upcoming computations that and that, for any U ∈ U , there exist bounded functions a U,V and b U,V such that These last relations can be proved similarly as (3.16). As a consequence, we immediatly deduce that for any (T, U ) ∈ T × U , and, using also Proposition 4.2 combined with the decay estimates of Proposition 10.1, These two estimates are sufficiently good in order to prove the two inequalities of the Lemma (recall that (L, L) ∈ T × U ). It then remains us to study the commutation of the frame fields with ∇ A ∇ A . If (T, U ) ∈ T × U , one has, since The first inequality of the Lemma can then be obtained using (11.4) and |∇ A k| ≤ |∇k|.
Remark 11.6. Note that the error terms given by Lemmas 11.4 and 11.5 are of size √ ǫ whereas the source terms of the Einstein equations are of size ǫ. For this reason, we will have to consider a hierarchy between the different energy norms considered for h 1 . In particular, when we will improve the bootstrap assumption on E 1+γ,1+γ , the terms given by the previous two lemmas will have to be bounded indenpendantly of C T U and C LL (respectively C LL ).
12. Improved energy estimates for the metric perturbations We start by the first one. For this, recall from Remark 7.5 that we can apply the second energy estimate of Proposition 7.5 to L J Z (h 1 ) for (a, b) = (γ, 1 + 2γ) and for any |J| ≤ N − 1. Consequently, by the Cauchy-Schwarz inequality and the bootstrap assumption (9.4), we obtain for all t ∈ [0, T [, where C > 0 is a constant which does not depend on C. We are now ready to prove the following result.
Then, if C is choosen sufficiently large and if ǫ is small enough, we have Proof. In view of the commutation formula of Proposition 4.9, the analysis of the source terms of the wave equation satisfied by L J Z (h 1 ) µν , which has been carried out in Section 11, and the inequality (12.1), we are led to bound sufficiently well the following integrals, defined for all multi-indices |J| ≤ N − 1.
Let us precise that, • Proposition 11.2 gives the terms I 0 and I J 3 • Proposition 11.1 gives the terms I 0 , I J 1 , I J 2 and I J 3 . • Proposition 11.3 gives I J 3 , I J 4 and I J 5 . • I J 6 is the source term related to the Vlasov field, it is estimated in Proposition 14.15. According to (12.1), the result follows if we prove, for any |J| ≤ N − 1 and all q ∈ 1, 6 , For later use, it will be useful to bound I 0 by an auxiliary quantity I 0 . Since 1 + 2γ ≤ 2, one easily finds that We fix |J| ≤ N − 1. Using the bootstrap assumption (9.6), we get By the crude estimate (1 + |u|) γ ≤ (1 + τ + r) 1−2δ and then bootstrap assumption (9.4), one obtains The Hardy type inequality of Lemma 3.11 yields We then deduce, using the bootstrap assumption (9.4) and 6δ ≤ 1 2 , that The next term can be estimated easily, using agin the bootstrap assumption (9.4), For I J 5 , the first step consists in applying the Hardy inequality of Lemma 3.11. For this reason, we cannot exploit all the decay in u = t − r in the exterior region (for simplicity, we do not keep all the decay in t − r that we have at our disposal in the interior region as well). We have Now, recall from (10.5) that Then, remark that, since 1 + |u| ≤ 1 + τ + r, so that, according to the bootstrap assumption (9.4) and the previous computations, Finally, the required bound on I J 6 is given by the assumptions of the proposition. This concludes the proof.
In order to improve the bootstrap assumption (9.4), one then only has to combine the previous result with Proposition 14.15, which will be proved in Subsection 14.3.
We now turn onE γ,2+2γ In the same way that we derive (12.1), one can prove using the third energy estimate of Proposition 7.5, the Cauchy-Schwarz inequality and the bootstrap assumption (9.5), that, for all t ∈ [0, T [, where C > 0 is a constant which does not depend on C. This last estimate, combined with Proposition 14.15 and the following result improves the bootstrap assumption (9.5) if ǫ is small enough and provided that C is choosen large enough.
Then, if C is choosen sufficiently large and if ǫ is small enough, we have Proof. The proof is similar to the one of Proposition 12.1. In view of the commutation formula of Proposition 4.9 and the estimates obtained on the error terms in Propositions 11.1-11.3, the result would follow if we bound by ǫ 2 (1 + t) 2δ the following integrals, defined for all multi-indices |J| ≤ N .
(1 + τ + r) 8 dxdτ, Note first that, using (12.2),I 0 ≤ I 0 ǫ 2 . We fix |J| ≤ N for the remainder of the proof. Using the bootstrap assumption (9.5), we directly obtain By the bootstrap assumption (9.7) and γ ≥ 3δ, we get Since 1 − 2δ ≥ 0, the bootstrap assumption (9.5) gives Using first the Hardy type inequality of Lemma 3.11 as well as the inequality 1 + |u| ≤ 1 + τ + r and then the bootstrap assumption (9.5) as well as 7δ ≤ 1, we obtain Applying the Hardy inequality of Lemma 3.11, we get Using (10.5) and ω 1+2γ 1+γ = ω 2+2γ γ 1+|u| , we obtain, using the previous computations, Finally, by the assumptions of the Proposition and Lemma 3.12, As a consequence, the constant C can be choosen independantly of C T U and C LL , provided that ǫ is small enough.
Then, there exist a constant C 0 independant of ǫ, C T U and C LL and a constant C independant of ǫ, such that, for all t ∈ [0, T [, Remark 12.5. Note that C T U has to be fixed sufficiently large compared to C but there is no restriction related to C LL .
In order to simplify the presentation of the following computations, all the constants hidden by will not depend on C T U nor on C LL . This convention will hold in and only in this subsection. We mention that all the energy norms which will be used here are defined in Subsection 3.7. We start by the following result. Proposition 12.6. There exists a constant C 0 independant of ǫ, C T U and C LL such that, for all t ∈ [0, T [, Proof. As these two estimates can be obtained in a very similar way, we only prove the second one. In order to lighten the notations, let us introduce φ J T U := χ( r t+1 )L J Z (h 1 ) T U for any |J| ≤ N and (T, U ) ∈ T × U . We can obtain from the first energy inequality of Proposition 7.5 and the Cauchy-Schwarz inequality that, According to Lemma 9.2, the smallness assumption on h 1 (t = 0) and the bootstrap assumption (9.7), we obtain, using also C T U ≥ 1, It then remains to combine these last four estimates.
Proposition 12.4 then ensues from the following two results.
Proof. According to the commutation formula of Proposition 4.9 and the result of Section 11, the proposition would follow if we could bound sufficiently well the quantities J J k defined below, for any multi-index J satisfying |J| ≤ N − 1 and any null components (T, U ) ∈ T × U . Those arising from the commutation of the wave operator with the cut off function (see Lemma 11.5), For J J 8 , we start by applying the Hardy inequality of Lemma 3.11. For this reason, we cannot use all the decay in t − r in the exterior region. We have where I 0 is defined and bounded by ǫ 2 in (12.2) and Since 1 − 4δ + γ > 0, we get from the bootstrap assumption (9.4) that Finally, Lemma 3.12, combined with the bootstrap assumption (9.4) and γ ≥ 3δ, gives Proof. The proof is similar to the one of Proposition 12.7. According to the commutation formula of Proposition 4.9, Propositions 11.1-11.3 and Lemma 11.4-11.5, it is suffient to bound by C 0 ǫ + Cǫ 2 (1 + t) 2δ the following integrals, defined for any |J| ≤ N and For J J 8 , as previsouly for similar integrals, we cannot keep all the decay in t − r when we apply the Hardy inequality of Lemma 3.11 (the problem comes from the exterior region). We have, since 1 ≥ 2δ, on the domain of integration of all these integrals. Since |∇(L J Z (h 1 ) T U )| |∇L J Z (h 1 )| + 1 r |L J Z (h 1 )| and 1 + τ + r 1 + |τ − r| for all r ≤ τ +1 2 , we have We also have Applying the Hardy type inequality of Lemma 3.11 and using the bootstrap assumption (9.5), we get Finally, the bootstrap assumption (9.5) gives 12.3. LL-energy. The purpose of this subsection is to prove the following result which, combined with Proposition 14.15, improves the bootstrap assumption (9.8) provided that ǫ is small enough and C LL choosen large enough.
Proposition 12.9. Assume that the following estimate holds Then, there exist a constant C 0 independant of ǫ and C LL and a constant C independant of ǫ, such that, Remark 12.10. For the conclusion of the previous proposition, it was crucial that C and C T U were fixed independantly of C LL (see Remarks 12.3 and 12.5).
In order to simplify the presentation of the following computations, all the constants hidden by will not depend on C LL . This convention will hold in and only in this subsection. The following result is the first step of the proof.
Proposition 12.11. There exists a constant C 0 independant of ǫ and C LL , such that, for all t ∈ [0, T [, Proof. In order to lighten the notations, let us introduce φ J := χ( r t+1 )L J Z (h 1 ) LL for any |J| ≤ N . We can obtain from the second energy inequality of Proposition 7.5 and the Cauchy-Schwarz inequality that, According to Lemma 9.2, the smallness assumption on h 1 (t = 0) and the bootstrap assumption (9.8), we obtain It then remains to combine these last four estimates.
We are then led to prove the following proposition.
Proof. Let us point out that C LL will only appear when we will use the bootstrap assumption (9.8). In order to prove this result, we are led to bound sufficiently well the following spacetime integrals, where the multi-index J will satisfy |J| ≤ N . Those coming from the commutation of the wave operator with the cut off function (see Lemma 11.5), Those coming from the commutation of the contraction with LL and the wave operator (see Lemma 11.4), Those coming from the contraction of g L I Z (h 1 ) µν with L µ L ν , More precisely, • Proposition 11.1 gives us the terms L 5 , L J 6 and L J 7 . • Proposition 11.2 leads us to control L 5 and L J 6 . • Proposition 11.3 gives the terms L J 6 and L J 8 . • L J 9 is the source term related to the Vlasov field, it is estimated in Proposition 14.15.
We start by the easiest ones, L 5 , L J 6 , L J 7 , L J 8 and L J 9 . First, according to (12.2), the hypothesis (12.7) and (12.5), We obtain from Lemma 3.12, the bootstrap assumption (9.7) and 2δ < 1 − 2δ, that According to the bootstrap assumption (9.8), we have We now focus on L J 1 , L J 2 , L J 3 and L J 4 . Since these integrals are of size ǫ (and not ǫ 2 ), we cannot use the bootstrap assumption (9.8) in order to control them as it would give us a bound larger than C LL ǫ(1 + t) δ . We will use several times the inequality 1 + τ + r ≤ 5r, which holds for all r ≥ τ +1 4 (and then on the domain of integration of each of these integrals). Using the inequality |∇(L J Note also that Consequently, applying the Hardy type inequality of Lemma 3.11 and using the bootstrap assumption (9.5), we get, since 1 − 2δ ≥ γ and 2δ < γ, Finally, as (1 + |u|) 1−γ ≤ (1 + τ + r) 1−γ , we obtain, using Lemma 3.12, the bootstrap assumption (9.7) and 2δ < γ, that The proof of Proposition 12.9 follows directly from Propositions 12.11 and 12.12, which concludes this section. 13. Improvement of the bootstrap assumptions on the particle density 13.1. General scheme. In this section we prove the following proposition.
Proposition 13.1. There exist an absolute constant C 0 > 0 and a constant 19 This improves in particular the bootstrap assumptions (9.1)-(9.3) if ǫ is small enough and provided that C f is choosen large enough.
Remark 13.2. One can check during the upcoming computations that the initial decay hypotheses on f could be lowered. We made the choice to simplify the presentation and then to work with energy norms weighted by z a , where the exponent a is as simple as possible.
In order to unify the proof of these three inequalities, we introduce for any multi-index |I| ≤ N the quantity (13.4) ℓ |I| := ℓ + 3 = 2 3 N + 9, |I| ≤ N − 5, According to the energy estimate of Proposition 8.1, we have where C is an absolute constant, which in particular does not depend on C f . In view of • the bootstrap assumptions (9.1)-(9.3), which give Proposition 13.1 is implied by the following two results.
Proposition 13.3. Let I be a multi-index of length |I| ≤ N . Then, Proposition 13.4. Let I be a multi-index of length |I| ≤ N . Then, 13.2. Proof of Proposition 13.3. Since the weights z are preserved by the flat relativistic transport operator T η , i.e. η αβ w α ∂ β (z) = 0, we have, using the notations introduced in Subsection 5.1, By a direct application of Lemmas 3.7 and 3.8 we have and recall from Proposition 10.7 that 20 We can then bound the first term of (13.6) by using (5.36) and the last two ones by applying Lemma 5.13, so that we obtain, since |w L | ≤ |v||w L |, We then deduce that The result ensues from the bootstrap assumptions (9.2) and (9.3).
13.3. Proof of Proposition 13.4. The starting point consists in bounding the commutator T g , Z I (f ) by a linear combination of the terms listed in Proposition 5.14. Then, in order to close the energy estimates and to deal with the weak decay rate of the metric, we will have to pay attention to the hierarchies related to the weights z which have been built into the Vlasov energy norms E ℓ+3 Before performing the proof, let us explain the strategy, which will be illustrated by the treatment in full details of the integral arising from the two families of error terms where Z ∈ P 0 , |J| + |K| ≤ |I|, |K| ≤ |I| − 1 and • either K P < I P • or K P = I P and J T ≥ 1, so that Z J contains at least one translation ∂ µ . We will then have to bound sufficiently well Apart for the error terms S K I,1 and S K I,2 , there is two cases to consider.
Step 1: if all the metric factors 21 can be estimated pointwise, e.g. ∇L J Z (h 1 ) for E J,K I,1 and ∇L J Z (h 1 ) LL for E J,K I,10 , so that |J| ≤ N − 5 according to Propositions 10.1 and 10.6. Then, the particle density is estimated in L 1 through the following result.
Lemma 13.5. Consider Z ∈ P 0 and let I and K be two multi-indices such that |I| ≤ N , |K| ≤ |I| − 1 and K P ≤ I P . Then, • Otherwise K P = I P and we have E Proof. This directly ensues from the fact that ∇ Z K (respectively Z Z K ) contains K P (respectively at most K P + 1) homogeneous vector fields and that ℓ |I| ≤ ℓ |K|+1 since |I| ≥ |K| + 1.
We need to consider two subcases for the most problematic terms, the quadratic and some of the cubic ones (see Proposition 5.14), in order to deal with a non integrable decay rate.
• If Z K contains less homogeneous vector fields than Z I , i.e. K P < I P , then the terms containing the factor Z Z K f are good since we control the energy norm of z ℓ |I| − 2 3 I P Z Z K f and the pointwise decay estimates on the metric provide an integrable decay rate. For I, we obtain from the pointwise decay estimates of Proposition 10.1, Lemma 13.5 and the bootstrap assumptions (9.1)-(9.3), For the remaining quadratic and cubic terms, which contain the factor ∇ Z K f , the pointwise decay estimates on the metric do not provide an integrable decay rate. The idea is to take advantage of the fact that we control the L 1 norm of z ℓ |I| − 2 3 I P + 2 3 ∇ Z K f and to gain decay through the extra weight z − 2 3 and Lemma 3.7. For J , we use Proposition 10.6, the inequality z − 2

3
(1 + |t − r|) − 2 3 which comes from Lemma 3.7, that δ ≤ γ < 1 6 , Lemma 13.5 and the bootstrap assumptions (9.1)-(9.3). We have In summary, we have proved first that . For I and J , this means that Z J contains a translation ∂ µ and that we can use the improved pointwise decay estimates of Proposition 10.8. We then get, using also Lemma 13.5 and the bootstrap assumptions (9.1)-(9.3), For I, as we merely control the energy norm of z ℓ |I| − 2 which comes from (3.21), so that In summary, we have proved first that √ ǫ|v| 1 + τ + r and then we have applied Lemma 13.5.
Step 2: if one of the metric factor cannot be estimated pointwise. In that case, the error term considered contain a factor where h 1 has been differentiated too many times so that we cannot apply Propositions 10.1 and 10.6 anymore. For J , this means that |J| ≥ N − 4. For I, we could have dealed with the cases |J| ∈ {N −4, N −3} during the first step but for simplicity we treat them here. Since |J| + |K| ≤ |I| ≤ N , we necessarily have |I| ≥ N − 4 and |K| ≤ 4 ≤ N − 9, so that the Vlasov field can be estimated pointwise. Note also that if |J| = N then |I| = N . Moreover, since ℓ |I| + 3 = ℓ |K|+1 , we will be able to gain decay through the weight z and Lemma 3.7 using |w L | |v|z 2 (1+t+r) 2 or 1 z 1+|t−r| . For I, we get, applying the Cauchy-Schwarz inequality in (τ, x) and since |w L | ≤ |w L ||v|, For J , we have (τ + r) Remark 13.6. We point out that E J,K I,10 is the most problematic term and that its treatment is more complicated than the ones of the other error terms. In particular, it is this term which prevents us to prove that E We are then led to prove the following lemma, which will also be useful for all the other error terms.
Since Z Z K contains at most I P + 1 homogeneous vector fields, |K| ≤ 5 ≤ N − 8 and ℓ |I| + 3 = ℓ + 3 = ℓ |K|+1 , we obtain from (9.9) and the bootstrap assumption (9.1) that which give us The bound on A K I can be obtained in the same way using this time that ∇ Z K contains at most I P homogeneous vector fields.
We can then bound I using the bootstrap assumptions (9.5). For any |J| ≤ N , The treatment of E J,K I, 10 , and then J , is much different for the case |J| = N than for N − 4 ≤ |J| ≤ N − 1. In both cases, we need to use an energy norm related to special components of h 1 in order to close the energy estimates. Assume first that |J| = N , which implies |I| = N . Then, using sup r∈R + 1+τ +r 1+|τ −r| 1 + τ , γ ≤ 1 16 and the bootstrap assumption (9.8), we obtain We now turn on the case N − 4 ≤ |J| ≤ N − 1. Apply first the inequality (3.14), so that .
Then, we bound A K I by using Lemma 13.7 and we apply the Hardy inequality of Lemma 3.11. Note that once again we need to be careful since we cannot use all the decay in u = τ − r in the exterior region. We obtain Fix now |J 0 | ≤ N and use the estimate (10.5), which was obtained using the wave gauge condition, in order to get t 0 Στ Using first that 1 + τ ≤ 1 + τ + r, δ ≤ γ, γ ≤ 1 + 1 8 and then the Hardy inequality of Lemma 3.11, we get We then deduce from the bootstrap assumption (9.5) and 4δ < 1 that I ǫ. Finally, as γ ≤ 1 8 , Lemma 3.12 combined with the bootstrap assumption (9.7) and γ > 3δ give t 0 Στ We then deduce from the previous estimates that J ǫ We now analyse the other error terms.
13.3.1. The terms arising from the source terms. Since T g (f ) = 0 we have Z I 0 (T g (f )) = 0 for any |I 0 | < |I| and all the error terms of the form (5.42) are equal to 0. 13.3.2. The terms which do not contain h 1 . We start by dealing with the error terms z ℓ |I| − 2 3 I P S K I,0 and z ℓ |I| − 2 3 I P S K I,00 since their treatment is different from the other ones.
Lemma 13.8. Let K be a multi-index satisfying |K| ≤ |I| − 1 and K P ≤ I P . Then, for any Z ∈ P 0 , Proof. As the Schwarzschild mass satisfies M √ ǫ, we have It remains to use the bootstrap assumption (9.1), (9.2) or (9.3).

13.3.3.
A sufficient condition for Proposition 13.4 to hold.
The two examples treated just before suggest us to prove the following three results, where we use the notations introduced in Definition 5.16. The first two ones concern the case where all the metric factor can be estimated pointwise. In the last result, we deal with the case where one of the h 1 factor has to be estimated in L 2 . Let us start by the easiest terms.
Lemma 13.9. Let Q, M , J and K be multi-indices satisfying |Q|+|M |+|J|+|K| ≤ N −5, |K| ≤ |I| − 1 and K P ≤ I P . Fix also Z ∈ P 0 . If for all (τ,   Proof. This follows from the definition of the quantities considered here and from the inequality z 2 3 ≤ (1 + τ + r) 2 3 , so that z ℓ |I| − 2 3 I P S J,K I,1 + S J,K I,2 + E Q,J,K I,12 + E Q,J,K I,13 Recall now the definition (3.36) of the norm E 1 8 , 1 8 [·], so that, using Lemma 13.5, the integrals considered in the statement of the lemma can be bounded by and it remains to use the bootstrap assumptions (9.1)-(9.3).
We now focus on the problematic terms. Those for which we need to use our hierarchy related to the weight z and the number of homogeneous vector fields composing Z I and Z K .
Fix also Z ∈ P 0 and define (1 + τ + r) Then, Proof. The proof is similar to the one of the previous lemma. Note that if K P < I P , G + A J,K I,11 · z ℓ |I| − 2 3 I P + 3 2 |∇ Z K f |. 22 Recall that we cannot have K P = I P in the error term E J,K I, 11 .
We now prove a similar result for the error terms containing a high order derivative of h 1 . Similarly, we have that . It then remains to remark that we necessarily have |K| ≤ 4 and to apply Lemma 13.7.
We consider now the first three terms. If K P < I P , we have which give the required bounds. If K P = I P , then J T ≥ 1 so that we can use the improved decay estimates given by Proposition 10.8. This leads to (1 + t + r) 2 3 A J,K I,2 + A J,K I,3 √ ǫ|v| We now treat the remaining terms, using again the pointwise decay estimates of Propositions 10.1 and 10.6 as well as the ones of Proposition 10.8 when J T ≥ 1. We have, using the inequality (1 + |τ − r|) 2 3 z 2 3 , which comes from Lemma 3.7, and then 2ab ≤ a 2 + b 2 , .
Otherwise we have J T ≥ 1 so that A J,K I,6 + A J,K I,9 = √ ǫ |v||w L | and we have then obtained the expected bounds when K P < I P . Similarly, one obtains This leads to the required bounds since z − 2
It remains to prove that the hypotheses of Lemma 13.11 hold.
Proposition 13.14. Let K be a multi-index such that |K| ≤ |I| − 1 and K P ≤ I P . Proof. Recall that we already dealt with the term associated to A J,K I,10 when we have bounded J (see (13.8)). We also already treated the integral associated to A J,K I,1 but we will repeat the proof here. We will oftenly use that 1 + |u| ≤ 1 + τ + r as well as the inequalities (13.9) 1 which come Lemma 3.7. We start by the terms of degree 1 in h 1 , i.e. the quadratic terms and some of the terms arising from the Schwarzschild part. We obtain by using (13.9) Similarly, we have Finally, using the wave gauge condition (10.5), there holds, as 1 + |τ − r| ≤ 1 + τ + r We now study the remaining terms. Note that without loss of generality, we can assume that |M | ≤ N − 5. Since |Q| ≤ N − 5 or |M | ≤ N − 5, we have, using the pointwise decay estimates of Proposition 10.1 and (13.9), Combining all the previous estimates, we are then led to prove that for all |I 0 | ≤ N , (1 + τ + r) 1−2δ (1 + |u|) 4 ω When we will apply the Hardy inequality of Lemma 3.11 in the upcoming computations, we will not be able to exploit all the decay in u = τ − r in the exterior region. Using first the Hardy inequality and then the wave gauge condition (10.5), we have (1 + |u|) 2 ω 1+2γ 1+γ dxdτ.
Moreover I j , J, Q and M 1 satisfy the following condition (1) either I P j < I P i , (2) or I P j = I P i and then J T ≥ 1, Q T + M T 1 ≥ 1. For the term ∇ λ L J Z H (w, w) · w λ w 0 ∂ vq F [ Z I j f ], J and I j satisfy the improved condition |J| + |I j | ≤ |I i | − 1 and I P j < I P i .

Remark 14.3.
Notice that if |I i | = N − 5, then A q i = 0 for all 1 ≤ q ≤ |M|. Proof. One only has to apply the commutation formula of Proposition 5.10 to Z I i f and to replace each derivatives of the Vlasov field Z K f , for |K| = N − 5, by the corresponding component of F or W . If |K| = N − 5, we replace it by the corresponding component of F for the following reason. In the terms listed on the Lemma, a derivative is applied to the components W k . Hence, if |K k | ≤ N − 6, we are able to rewrite ∂ x µ W k and ∂ v i W k as a combination of components of W , which will be important later.
The goal is to obtain an L 2 -estimate on F . For this, let us split it in F := F hom + F inh , where T g (F hom ) + A · F hom = 0, F hom (0, ·, ·) = F (0, ·, ·), T g (F inh ) + A · F inh = B · W, F inh (0, ·, ·) = 0. and then prove L 2 estimates on the velocity average of F hom and F inh . To do it, we will schematically establish that F inh = KW , with K a matrix such that E[KKW ] do not growth too fast, and then use the pointwise decay estimates on v |W |dv given by (9.10) to obtain the expected decay rate on v |F inh |dv L 2 x . For v |F hom |dv L 2 x , we will make crucial use of the Klainerman-Sobolev inequality of Proposition 3.15 so that we will need to commute the transport equation satisfied by F hom and prove L 1 -bounds such as we did in Section 13.
It will be convenient to denote, as for F , the components F hom i and F inh i of F hom and F inh as follows, Remark 14.4. Contrary to [16], we kept, as in [8], the v derivatives in the statement of Lemma 14.2 in order to take advantage of the good behavior of radial component of ∇ v F . If we had already transformed the v derivatives, we would have obtained terms such as x j r (t − r)∂ x j F from (∇ v F ) r (see Lemma 3.9). We would then have to deal with factors such as t 3 |x| 3 during the treatment of the homogeneous part F hom (apply three boost to x k |x| ), when we will have to commute at least three times T g (F hom ) + A · F hom = 0.
However, this creates two new technical difficulties compared to the strategy of [16] and we will circumvent it by following [8]. The first one concerns F hom and will lead us to consider a new hierarchy (see Subsection 14.1). The other one concerns certain source terms of the transport equation satisfied by F inh , which contain derivatives of F inh . Because of the presence of top order derivatives of h 1 , we will not commute this equation and these derivatives have to be rewritten as a combination of components F inh and controlled terms, which will be derivatives of F hom .
14.1. The homogeneous system. In order to obtain L ∞ , and then L 2 , estimates on v |F inh |dv, we will have to commute at least three times the transport equation satisfied by each component of F inh . However, if for instance |I i | = N − 4, we need to control the L 1 norm of Z K F hom [ Z I j f ], with |K| = 4 and |I j | = N − 5, to bound Z I F hom [ Z I i f ] L 1 x,v , with |I| = 3. We then consider the following energy norm (recall that ℓ = 2 3 N + 6), We have the following commutation formula.
Moreover K, J j , J, Q and M 1 satisfy the following condition (1) either K P + I P j < I P + I P i , (2) or K P + I P j = I P + I P i and then J T ≥ 1, Q T + M T 1 ≥ 1. For the term (14.2), J and K satisfy the improved condition K P + I P j < I P + I P i . Proof. Let i ∈ 1, |M| and |I| ≤ N + 3 − |I i |. The starting point is the relation According to Proposition 5.10, the error terms arising from the commutator T g , Z I F hom [ Z I i f ] are • such as those listed in the lemma, with I j = I i . Note that the conditions on |J| and |M 1 |+|M 2 |+|Q| follows from |J|+|K|, |M 1 |+|M 2 |+|Q|+|K| ≤ |I| ≤ N +3−|I i | ≤ 8 and N ≥ 13.
• Or such as Z I 0 T g (F hom [ Z I i f ]) , with |I 0 | < |I| and I P 0 < I P .
The analysis of the other source terms is similar to the one made in order to derive the commutation formula of Proposition 5. 10  dxdτ.
Using |T g (z)| ≤ √ ǫ|v|z 1+t+r + √ ǫ|w L |z 1+|t−r| (see (13.7)) and (3.36), we obtain Then, Definition (14.1) of E F hom and the bootstrap assumption on it lead to The integral Z I,I i can be bounded similarly using Corollary 14.6 instead of (13.7). We then deduce from (14.1) and the last estimates that there exists a constant C 0 independant of C F such that which improves the bootstrap assumption if ǫ is small enough and C F choosen large enough. This implies that T 0 = T . The pointwise decay estimates can then be obtained from the Klainerman-Sobolev inequality of Proposition 3.15 and the fact E F hom gives a control on the derivatives up to third order of z ℓ−2− 2 3 (I P i +I P ) Z I F hom Z I i f , for any |I| + |I i | ≤ N .
14.2. The inhomogenous system. To derive an L 2 estimate on F inh , we cannot commute the transport equation because B contains top order derivatives of h 1 . We then need to rewrite the derivatives of F inh , kept in the matrix A in order to use the full null structure of the system, in terms quantities that we can control. More precisely, we will use the following result.
Lemma 14.8. Let i ∈ 1, |M| such that |I i | ≤ N − 1 and 0 ≤ µ ≤ 3. Then, For the v derivatives, there holds Proof. Recall that F = F hom +F inh and note that for any Z ∈ P 0 and N −5 ≤ |I i | ≤ N −1, In order to rewrite the transport equation satisfied by F inh , we will then need to consider a bigger vector valued field than W . Moreover, in order to take advantage of the hierarchies that we identified in the commuted Vlasov equation, we will work with a slightly different quantity than F inh .  We define Y as a the vector valued field of length l Y containing the following quantities We are now ready to prove the following two results.
Moreover, A and B are such that, if i ∈ 1, |M| , T F (F inh z,i ) can be bounded by a linear combination of terms of the form √ ǫ|v| |Y |.
Using the Cauchy-Schwarz inequality in v, we then obtain, as i ≤ |M n |, (1 + τ + r) As z ≤ z 2 , we obtain that I + I ǫ dx, we obtain from the pointwise decay estimate on v z 4 |Y ||v|dv and the bootstrap assumption on E n F inh that We then deduce that I B ǫ (1 + t) 1+ 3 2 δ otherwise, so that If ǫ is small enough, this improves the bootstrap assumptions on E N −1 F inh and E N F inh .
14.3. The L 2 estimates. We start by estimating the L 2 norm of R 3 v z| Z K f |dv. Proof. Assume first that |I| ≤ N − 4. Then, using the Cauchy-Schwarz inequality in v and then the pointwise decay estimate (9.10) as well as the bootstrap assumption (9.2), we get Otherwise |I| ≥ N − 3 and there exists i ∈ M such that We deduce that K ≤ K hom + K inh , where, using Proposition 14.7, Recall Definition 14.9 and that K · Y = F inh z . Hence, Using first the Cauchy-Schwarz inequality in v and then the pointwise decay estimates (14.8), I i = I as well as Proposition 14.13, we obtain We are now able to prove the following result.
The proof follows from (14.9)-(14.10) and the estimates obtained on J 1 , J 2 and J 3 .