Spatially homogeneous Friedmann–Lemaître–Robertson–Walker (FLRW) solutions constitute an infinite dimensional family of explicit solutions of the Einstein–massless Vlasov system with vanishing cosmological constant. Each member expands toward the future at a decelerated rate. These solutions are shown to be nonlinearly future stable to compactly supported spherically symmetric perturbations, in the case that the spatial topology is that of R3. The decay rates of the energy momentum tensor components, with respect to an appropriately normalised double null frame, are compared to those around Minkowski space. When measured with respect to their respective t coordinates, certain components decay faster around Minkowski space, while others decay faster around FLRW.

Standard homogeneous isotropic cosmological models in general relativity are described by the Friedmann–Lemaître–Robertson–Walker (FLRW) spacetimes
M=I×Σ,g=dt2+a(t)2gΣ,
(1.1)
where IR is an open interval, (Σ, gΣ) is a constant curvature manifold, and a: I → (0, ∞) is an appropriate scale factor. See Sec. 5.3 of Ref. 1 for more on FLRW spacetimes. This article concerns radiation filled FLRW cosmologies in which the constant curvature manifold is Euclidean space, (Σ,gΣ)=(R3,gEucl), with radiation described by spatially homogenous solutions of the massless Vlasov equation, and their stability properties.
Consider a 3 + 1 dimensional Lorentzian manifold (M,g) and let
P={(t,x,p)TMg(t,x)(p,p)=0}TM,
denote the mass shell of (M,g). Consider some local coordinates {t, x1, x2, x3} on M, where t is a time function, and let {t, xi, pμ} denote the corresponding conjugate coordinate system for TM, i.e., (t, xi, pμ) describes the point pμxμ|(t,x)TM. The massless Vlasov equation on (M,g), for a particle density f: P → [0, ∞), takes the form
p0tf+pixifpμpνΓμνipif=0,
(1.2)
where p0 is defined by the mass shell relation
gμνpμpν=0.
(1.3)
The Einstein–massless Vlasov system consists of Eq. (1.2) coupled to the Einstein equations
Ric(g)μν12R(g)gμν=Tμν,
(1.4)
where the energy momentum tensor takes the form
Tμν(t,x)=P(t,x)f(t,x,p)pμpνdetgp0dp1dp2dp3,
(1.5)
where indices are raised and lowered with respect to the metric g (so that, for example, p0 = g0μpμ). Here t is also denoted x0, Greek indices, such as μ, ν, range over 0, 1, 2, 3, and lower case Latin indices, such as i, j, k range over 1, 2, 3.

The spatially homogeneous FLRW family constitutes an infinite dimensional family of explicit solutions of the Einstein–massless Vlasov system (1.2)(1.5).

Define the manifold
M=(0,)×R3,
(1.6)
and consider the Euclidean metric
gEucl=(dx1)2+(dx2)2+(dx3)2,
on R3. Given xR3 and p=(p1,p2,p3)TxR3, define
|p|2=|p|gEucl2=(gEucl)ijpipj=(p1)2+(p2)2+(p3)2.
For any smooth, sufficiently decaying function μ: [0, ∞) → [0, ∞) such that μ ≢ 0, the metric g and particle density f defined by
g=dt2+a(t)2(dx1)2+(dx2)2+(dx3)2,f(t,x,p)=μ(a(t)4|p|2),
(1.7)
where
a(t)=t124ϱ314,ϱ=|p|μ(|p|2)dp,
(1.8)
define a solution of (1.2)(1.5) on (1.6). Such solutions are known as spatially homogeneous FLRW solutions.

Remark 1.1

(Decelerated expansion). Spacetimes with metric of the form (1.1) are said to undergo accelerated expansion if the scale factor a(t) satisfies a″(t) > 0 for all t, linear expansion if a″(t) = 0 for all t, or decelerated expansion if a″(t) < 0 for all t. Note that the solutions (1.7) and (1.8) undergo decelerated expansion. Previous stability works in cosmological settings have typically considered perturbations of spacetimes undergoing linear or accelerated expansion. See Sec. I D below.

Remark 1.2
(Anisotropic spatially homogeneous FLRW solutions). The spatially homogeneous FLRW family of solutions (1.7) and (1.8) lie inside a more general family of anisotropic spatially homogeneous FLRW solutions. Given any smooth, sufficiently decaying function F:R3[0,) such that F ≢ 0 and
R3F(p)pidp=0,i=1,2,3,andR3F(p)pipj|p|dp=ϱ3δij,i,j=1,2,3,
the metric and particle density
g=dt2+a(t)2(dx1)2+(dx2)2+(dx3)2,f(t,x,p)=F(a(t)2p1,a(t)2p2,a(t)2p3),
where
a(t)=t124ϱ314,ϱ=|p|F(p)dp,
also defines a solution of (1.2)(1.5) on (1.6). Since these solutions are, in general, not spherically symmetric (see Sec. II A), consideration here is restricted to the solutions of the form (1.7) and (1.8).

Remark 1.3

(T3 spatial topology). Each (1.7) and (1.8), and more generally the solutions of Remark 1.2, also define a spatially homogeneous FLRW solution of the Einstein–massless Vlasov system on the manifold (0,)×T3. This article will only be concerned with solutions on (0,)×R3. (See also Remark 3.5 for a further comment on the solutions on (0,)×T3.)

The components of the energy momentum tensor of (1.7) and (1.8), with respect to the Cartesian (t, x1, x2, x3) coordinate system, take the form
(T°)tt(t,x)=34t2,(T°)ij(t,x)=14t2δij,(T°)it(t,x)=0,
(1.9)
for i, j = 1, 2, 3, where (T°)ij(t,x)=gik(t,x)(T°)jk(t,x).

The spacetimes (M,g) expand from a spacelike singularity at t = 0 (with Kretchmann scalar |Rm|g2=32t4) and are future geodesically complete. See the Penrose diagram of Fig. 1. Such solutions feature in the Ehlers–Geren–Sachs Theorem,2 which in particular ensures that all solutions of the Einstein–massless Vlasov system for which f is isotropic and irrotational are either stationary or described by an FLRW metric as above.

FIG. 1.

Penrose diagram of the FLRW spacetime. Solutions of the massless Vlasov equation arising from compactly supported initial data on {t = 1} are supported in the darker shaded region. The solutions, of the coupled system, of Theorem 1.4 admit a similar Penrose diagram, where now the perturbation is supported in the darker shaded region.

FIG. 1.

Penrose diagram of the FLRW spacetime. Solutions of the massless Vlasov equation arising from compactly supported initial data on {t = 1} are supported in the darker shaded region. The solutions, of the coupled system, of Theorem 1.4 admit a similar Penrose diagram, where now the perturbation is supported in the darker shaded region.

Close modal

The main result of the present work concerns the future stability of the family (1.7) and (1.8).

Theorem 1.4

(Future stability of FLRW). Each spatially homogeneous FLRW solution is future nonlinearly stable, as a solution of the Einstein–massless Vlasov system (1.2)(1.5), to compactly supported spherically symmetric perturbations.

More precisely, consider an FLRW solution of the form (1.7) and (1.8) for some smooth, non-identically vanishing μ (decaying suitably so that ϱ is defined). For all compactly supported spherically symmetric initial data sufficiently close to that of (1.7) and (1.8) on Σ1 = {t = 1}, the resulting solution is future geodesically complete, is isometric to (1.7) and (1.8) after some retarded time and, in an appropriately normalised double null gauge, the (appropriately normalised) components of the energy momentum tensor decay to zero as (1.9) to leading order, and each metric component and Christoffel symbol either remains close or decays to its FLRW value to the past of this retarded time, with quantitative polynomial rates.

A more precise version of Theorem 1.4 is stated below in Sec. IV. See also the Penrose diagram in Fig. 1.

The proof of Theorem 1.4 is based on an understanding of decay properties of solutions of the massless Vlasov equation. As such, in addition to Theorem 1.4, decay properties of solutions of the massless Vlasov equation on a fixed FLRW background of the form (1.7) and (1.8) are given. See Theorem 3.2.

The assumption of compact support can be relaxed, and is made for convenience in order to localise the proof. Indeed, following,3 the perturbation of f, in the coupled problem, vanishes close to the centre of spherical symmetry at late times. The absence of singularities in future evolution therefore follows from the aforementioned quantitative decay properties, along with comparatively soft arguments including an extension principle around non-central points (see Ref. 4 or Theorem 2.3 below).

Remark 1.5

(Asymptotic stability of FLRW). Theorem 1.4, from the spacetime (as opposed to mass shell) perspective, can be viewed as the asymptotic stability of the FLRW solution (1.7) and (1.8) to spherically symmetric perturbations. There is residual double null gauge freedom (see Remark 2.1) which can be appropriately renormalised so that each metric component, Christoffel symbol and energy momentum tensor component is equal to its (possibly decaying) FLRW value, plus a remainder which decays strictly faster than its corresponding FLRW value. Such a renormalisation is not required in the Proof of Theorem 1.4 however, and thus is not further considered here.

Remark 1.6

(Comparison with decay rates on Minkowski space). The decay rates of solutions of the massless Vlasov equation on FLRW, obtained in Theorem 3.2, are later compared with those of the massless Vlasov equation on Minkowski space. When measured with respect to their respective t coordinates, certain components of the energy momentum tensor, with respect to an appropriately normalised double null frame, decay faster in Minkowski space, while others decay faster in FLRW. See Remark 3.4. Note also that solutions of the linear wave equation with localised initial data decay at the same rate on each of these two spacetimes.5 

Remark 1.7

(Birkhoff-type theorem). The Proof of Theorem 1.4 in particular contains a Birkhoff-type theorem for the system (1.2)(1.5) which ensures that any spherically symmetric solution with a regular centre, for which f is equal to its FLRW value f [defined by (1.7)], is locally isometric to the FLRW spacetime (M,g). See the final step in the proof of Theorem 5.1 (from which a general theorem can be extracted). Recall also the Ehlers–Geren–Sachs Theorem.2 

Remark 1.8
(Einstein–Euler with radiation equation of state). The metric g, defined by the former of (1.7) with a(t)=t12, also arises as a solution of the Einstein–Euler system, i.e., Eq. (1.4) with
Tμν=(ρ+p)uμuν+pgμν,
for unknowns g, four velocity uμ — satisfying gμνuμuν = −1 — pressure p, and density ρ, with p and ρ related by the radiation equation of state
p=ρ3.
The fluid variables for the FLRW metric take the form, in Cartesian coordinates,
u=t,ρ=34t2.
(1.10)

It is known that the solution (1.10) is unstable to formation of shock waves for the Euler equationsμTμν = 0 on the fixed FLRW background (M,g) (see Chap. 9 of Ref. 6).

Note also the linear Jeans instability for the Einstein–Euler system linearised around FLRW7 (see also the discussion in Ref. 8).

There have been previous works on the Einstein–Vlasov system in cosmological settings for Einstein equations with a positive cosmological constant,9 and also for perturbations of the vacuum Milne solution.10–12 Note also related works13–17 on the Einstein–Euler system. See also Refs. 18–21. Each of these works considers perturbations of spacetimes with metric of the form (1.1) undergoing accelerated expansion or linear expansion (see Remark 1.1). Contrast with the solutions (1.7) and (1.8) of the present work in which undergo slow decelerated expansion.

In the asymptotically flat setting, Minkowski space has been shown to be nonlinearly stable for both the massless and massive Einstein–Vlasov systems, first in spherical symmetry,3,22 and later to general perturbations.23–25 The present work is based on the approach of Ref. 3. Anti-de Sitter space has also been shown to be nonlinearly unstable as a solution of the Einstein–massless Vlasov system with a negative cosmological constant.26 

Though there is no direct non-relativistic analogue of the present problem, there is an infinite dimensional family of explicit, spatially homogenous, stationary solutions of the Vlasov–Poisson system. The works27–29 on Landau damping, and also the celebrated works30–32 on the torus T3, concern the stability of these families of solutions. Note however that one must apply the Jeans swindle to make sense of the problems considered in these works in the gravitational setting. See also Ref. 33.

Section II contains preliminaries on the Einstein–massless Vlasov system in double null gauge. Section III concerns further properties of the FLRW spacetimes, introduced in Sec. I B. A double null gauge in FLRW is defined in Sec. III A and, as a precursor to the proof of Theorem 1.4, properties of the massless Vlasov equation (1.2) and (1.3) on the FLRW background spacetimes are presented in Sec. III B. In Sec. IV A more precise version of Theorem 1.4 is formulated, and in Sec. V its proof is given.

This section concerns facts about the Einstein–massless Vlasov system in spherical symmetry. In Sec. II A the spherical symmetry assumption is introduced. In Sec. II B spherically symmetric double null gauges are introduced, and the residual spherically symmetric double null freedom is discussed. In Sec. II C the Einstein–massless Vlasov system in a spherically symmetric double null gauge is presented. In Sec. II D functions t and r are introduced. Finally, in Sec. II E the Cauchy problem is discussed and an extension principle, which will be used in the Proof of Theorem 1.4, is stated.

A solution (M,g,f) of (1.2)(1.5) is spherically symmetric if SO(3) acts by isometry on (M,g) and preserves f. More precisely, it is assumed that there is a smooth isometric action O:M×SO(3)M on (M,g) such that, for each pM, the orbit of the action satisfies either Orb(p) ≃ S2, or Orb(p) = {p} and, moreover, for each of the generators Ω1, Ω2, Ω3 of the SO(3) action, the corresponding flows Φi:(0,2π)×MM satisfy
f(x,p)=f(Φsi(x),(Φsi)*p),for all s(0,2π),i=1,2,3,
(2.1)
for all (x,p)P. It is assumed that the centre
Γ={xMO(x,s)=x for all sSO(3)},
consists of a single timelike geodesic, and that M\Γ splits diffeomorphically as
M\ΓQ×S2,
for a smooth two-manifold Q=(M\Γ)/SO(3). The quotient space M/SO(3), also denoted Q, is viewed as a manifold with boundary, with boundary identified with the centre Γ.

A double null gauge consists of functions u,v:QR such that the level hypersurfaces of u foliate Q by outgoing lines which are null with respect to the induced metric on Q, and the level hypersurfaces of v foliate Q by ingoing null lines. The level hypersurfaces of u and v lift to outgoing and incoming null cones of M and such a double null gauge can be complemented with local coordinates (θ1, θ2) on S2 to local coordinates (u, v, θ1, θ2) for M.

For a given double null gauge, the metric g can be written in double null form
g=Ω2dudv+R2γ,
(2.2)
where γ is the unit round metric on S2, where Ω is a function on Q and R:QR is the area radius function
R(u,v)=Area(Su,v)/4π,
where Su,v={(u,v)}×S2M. The area radius R extends regularly to 0 on the centre Γ. Define
λ=vR,ν=uR.

Remark 2.1
(Residual gauge freedom). There is residual double null gauge freedom present in (2.2). If ũ,ṽ:IR are increasing functions, with IR a suitable interval, then the change uũ(u), vṽ(v) preserves the double null form (2.2). Under such a change, the metric takes the form
g=Ω̃2(ũ,ṽ)dũdṽ+R(ũ,ṽ)2γ=Ω̃2(ũ(u),ṽ(v))dũdu(u)dṽdv(v)dudv+R(ũ(u),ṽ(v))2γ.
It thus follows that
Ω2(u,v)=Ω̃2(ũ(u),ṽ(v))ũ(u)ṽ(v),
(2.3)
R(u,v)=R̃(ũ(u),ṽ(v)),
(2.4)
ν(u,v)=ν̃(ũ(u),ṽ(v))ũ(u),
(2.5)
λ(u,v)=λ̃(ũ(u),ṽ(v))ṽ(v).
(2.6)

The inverse metric of (2.2) takes the form
g1=2Ω2(uv+vu)+R2γABθAθB,
where A and B range over 1 and 2. The nonvanishing Christoffel symbols of g are
ΓABu=2λRΩ2γAB,ΓABv=2νRΩ2γAB,ΓBvA=λRδBA,ΓBuA=νRδBA,Γuuu=ulogΩ2,Γvvv=ulogΩ2,
and
ΓBCA=ΓBCA,
where ΓBCA are the Christoffel symbols of (S2, γ).
For a given spherically symmetric double null gauge (u, v, θ1, θ2) for (M,g), one defines a corresponding conjugate coordinate system (u, v, θ1, θ2, pu, pv, p1, p2) for TM, whereby (u, v, θ1, θ2, pu, pv, p1, p2) defines the point
puu+pvv+pAθA(u,v,θ1,θ2)TM.
Such a coordinate system induces a coordinate system (u, v, θ1, θ2, pv, p1, p2) on the mass shell P, with pu defined by the mass shell relation (1.3), which takes the form
pu=R2γABpApBΩ2pv=L2Ω2R2pv,
(2.7)
where the angular momentum L: P → [0, ∞) is defined by
L(x,p)2=R4γABpApB.
(2.8)
The volume form on
Px={pTxMgx(p,p)=0}TxM,
in this induced coordinate system takes the form
R2detγpvdpvdp1dp2.
The condition (2.1) that f is spherically symmetric means that, in a given double null gauge, f can be written as a function — which, abusing notation slightly, is also denoted f:Q×[0,)×[0,)[0,) — of u, v, pv, and L=(R4γABpApB)12,
f(x,p)=f(u,v,pv,L).
(2.9)
The energy momentum tensor T on M takes the form
T=Tuu(u,v)du2+2Tuv(u,v)dudv+Tvv(u,v)dv2+TAB(u,v)dθAdθB.
In double null gauge, the components take the form
Tuu(u,v)=Ω4R2π200f(u,v,pv,L)pvLdLdpv,
(2.10)
Tuv(u,v)=Ω4R2π200f(u,v,pv,L)puLdLdpv,
(2.11)
Tvv(u,v)=Ω4R2π200f(u,v,pv,L)(pu)2pvLdLdpv,
(2.12)
with pu = pu(u, v, pv, L) defined by (2.7). It moreover follows from the mass shell relation (2.7) that
γABTAB=R24Ω2Tuv.
The Einstein equations (1.4) in spherical symmetry take the form of the following system of equations for (Ω, R, f),
uv(R2)=Ω22+R2Tuv,
(2.13)
uvlogΩ2=Ω22R21+4Ω2uRvRTuvΩ24R2γABTAB,
(2.14)
u(Ω2uR)=12RΩ2Tuu,
(2.15)
v(Ω2vR)=12RΩ2Tvv,
(2.16)
where g=R2γ. The mass shell relation (1.3) takes the form
Ω2pupv+R2γABpApB=0.
(2.17)
In particular, Ω2Tuv = R2γABTAB, or 4Ω−2Tuv = R−2γABTAB, and so Eq. (2.14) can be rewritten
uvlogΩ2=Ω22R21+4Ω2uRvR2Tuv,
(2.18)
The Vlasov equation (1.2) in spherical symmetry takes the form
X(f)puuf+pvvfR2γABpApB2νRΩ2+(pv)2vlogΩ2pvf2λpv+νpuRLLf=0.
(2.19)
The angular momentum L: P → [0, ∞) defined by (2.8) is preserved by the geodesic flow
X(L)=0.
(2.20)

The system (2.13)(2.19) is equivalent to the Einstein–massless Vlasov system (1.2)(1.5) in the sense that any spherically symmetric solution (M,g,f) of (2.13)(2.19) defines a solution (Q,Ω2,R,f) of the reduced system (2.13)(2.19) and, conversely, any solution (Q,Ω2,R,f) of the reduced system (2.13)(2.19) defines a solution of (1.2)(1.5) with M=Q×S2 and g defined by (2.2).

See Ref. 34 for more on the spherically symmetric Einstein–massless Vlasov system.

Define functions t and r by
t(u,v)=14(v+u)2,r(u,v)=vu.
When the FLRW spacetime (1.7) and (1.8) is expressed in the double null gauge introduced in Sec. III A below, the functions t and r coincide with those of Sec. I B.
The Einstein–massless Vlasov system (1.2)(1.5) admits the following local existence theorem. Initial data consists of a Riemannian 3 manifold (Σ,ḡ), along with a symmetric (0,2) tensor K on Σ, and a function f1: TΣ → [0, ∞), satisfying constraint equations
R̄|K|ḡ2+(tr̄K)2=2detḡf1|p|ḡdp1dp2dp3,div̄Kiitr̄K=detḡf1pidp1dp2dp3,
where R̄ is the scalar curvature and ̄ is the Levi–Cività connection of ḡ. A development of initial data (Σ,ḡ,K,f1) is a globally hyperbolic spacetime (M,g) which admits Σ as a Cauchy hypersurface, with induced first and second fundamental form ḡ and K respectively, such that the restriction of f to the mass shell over Σ, P|Σ, coincides with f1, when TΣ is appropriately identified with P|Σ. See Ref. 9 for more on the Cauchy problem for the Einstein–Vlasov system.

Theorem 2.2

(Local well posedness of the Cauchy problem for the Einstein–massless Vlasov system35,36). For any smooth initial data set (Σ,ḡ,K,f1) for the Einstein–massless Vlasov system, as above, there exists a unique maximal development (M,g,f) solving the system (1.2)(1.5). Moreover, if (Σ,ḡ,K,f1) is spherically symmetric then the maximal development (M,g,f) is spherically symmetric.

See also Ref. 34 for a version of Theorem 2.2 in a more general spherically symmetric setting.

The following extension principle for spherically symmetric solutions of the system (1.2)(1.5), concerning non-central points, will be used in what follows. It is assumed that (M,g,f) is the maximal development of a smooth spherically symmetric initial data set (Σ,ḡ,K,f1) for the Einstein–massless Vlasov system (1.2)(1.5), and that Q is the quotient manifold of M by the action of the SO(3) isometry, as in Sec. II A above.

Theorem 2.3
(Extension principle, around non-central points, for the spherically symmetric Einstein–massless Vlasov system4). Let Q be as above. Let (u*, v*) be such that there exists U < u* and V < v* such that the characteristic diamond DU,Vu*,v*={Uu<u*,Vv<v*} is contained in Q,
DU,Vu*,v*Q.
(2.21)
Assume also that there exists 0 < R0 < R1 such that
Uu*Vv*Ω2(u,v)dudv<,andR0R(u,v)R1,for all (u,v)DU,Vu*,v*,
and that f(u, v, ⋅, ⋅) is compactly supported for all (u,v)DU,Vu*,v*. Then (u*,v*)Q.

For a proof of Theorem 2.3 see Sec. IV C of Ref. 4 [where one has to replace the massive mass shell relation with its massless analogue (2.17), which does not affect the proof]. Theorem 2.3 is preferred to the softer extension principle of Ref. 37 in view of the presence of the anti-trapped surfaces in the spacetimes under consideration (see Remark 3.1 below).

Recall the FLRW spacetimes introduced in Sec. I B. In Sec. III A 1 double null gauge is introduced in each of these spacetimes. Section III B concerns properties of the massless Vlasov equation on these spacetimes. Though the main theorem of Sec. III B, Theorem 3.2, is not used in the proof of Theorem 1.4, its proof is presented as a simple precursor to that of Theorem 1.4. There is no symmetry assumption required, however, for Theorem 3.2.

Recall that the FLRW spacetime (M,g) takes the form
M=(0,)×R3,g=dt2+t4ϱ312(dx1)2+(dx2)2+(dx3)2.
For t ∈ (0, ∞), define
Σt={t}×R3M.
Define double null coordinates
u=t12r2,v=t12+r2,where r=4ϱ314(x1)2+(x2)2+(x3)2,
and note that, since t ≥ 0 and r ≥ 0 in M,
v0,vu,vu,in M.
The FLRW metric g in the above double null gauge takes the form
g=4tdudv+tr2γ,t=14(v+u)2,r=vu,
(3.1)
defined on the quotient manifold
Q={(u,v)R2v0,vu,vu}.
The metric g can be written
g=Ω2dudv+R2γ,whereΩ2=4t,R=t12r,
and moreover
λ=vR=v=r2+t12,ν=uR=u=r2t12,detg=2t2r2detγ.

Remark 3.1

(Anti-trapped spheres). Note that the FLRW spacetime (M,g) contains anti-trapped spheres of symmetry: if (u,v)Q is such that u < 0, then λ(u, v) > 0 and ν(u, v) > 0.

The non-vanishing Christoffel symbols of (3.1) take the form
Γ̊ABu=rv2t12γAB=1212t12+1rr2γAB,Γ̊ABv=ru2t12γAB=1212t121rr2γAB,
Γ̊uuu=Γ̊vvv=t12,Γ̊BvA=12t12+1rδAB,Γ̊BuA=12t121rδAB,Γ̊BCA=ΓBCA,
where ΓBCA are the Christoffel symbols of (S2, γ).
It follows from (1.9) that the null components of the energy momentum tensor of f satisfy
T°uu=1t,T°uv=12t,T°vv=1t.
(3.2)
The mass shell relation (2.17) takes the form
4pupv+r2γABpApB=0.
(3.3)
Recall the expression (2.7) for pu, and define
p̊u=r2γABpApB4pv.
Where there is no ambiguity, for example in Sec. III B, p̊u will also be denoted pu.
Given a solution (Q,Ω2,R,f) of the spherically symmetric Einstein–massless Vlasov system in double null gauge (2.13)(2.19), one identifies with the spatially homogeneous FLRW solutions (1.7) and (1.8) using values of the coordinates (u, v) coordinates for Q, and the values of the coordinates (u, v, pv, L) for Q×[0,)×[0,), along with the above double null gauge for the spatially homogeneous FLRW solutions. For example
(Ω2Ω2)(u,v)=Ω2(u,v)(v+u)2.
Note in particular that, if f = f, then
Ω4R2Tuu=Ω4R2T°uu,Ω2R4Tuv=Ω2R4T°uv,R6Tvv=R6T°vv.
The asymmetry in u and v arises from the choice of parameterising the mass shell by pv, rather than pu.

The Proof of Theorem 1.4 is based on the following proof of decay of components of solutions of the massless Vlasov equation on a fixed FLRW backround.

Theorem 3.2
(Solutions of massless Vlasov equation on a fixed FLRW background with R3 spatial topology). Let f be a solution of the massless Vlasov equation (1.2) on ((0,)×R3,g), where g is the FLRW metric (3.1), such that f1 = f|{t=1} is compactly supported. The components of the energy momentum tensor satisfy
Tuuf1Lt2,Tuvf1Lt3,Tvvf1Lt4,
(3.4)
for t ≥ 1. Moreover, there exist U0U1 such that
π(supp(f)){U0uU1},
(3.5)
where π:PM denotes the natural projection.

The support of such solutions are depicted in Fig. 1.

The proof of Theorem 3.2 is based on properties of null geodesics in FLRW. Given a null geodesic γ:[1,)M, write the tangent vector to γ as
γ̇(s)=pu(s)u+pv(s)v+pA(s)θA.
The geodesic equations for pu(s) and pv(s) take the form
ṗu(s)+t12pu(s)2+1212t12+1rr2γABpApB=0,
(3.6)
ṗv(s)+t12pv(s)2+1212t121rr2γABpApB=0.
(3.7)
One has equations for t(s) = t(γ(s)) and r(s) = r(γ(s)) along the geodesic
ṫ(s)=t(s)12pv(s)+pu(s),ṙ(s)=pv(s)pu(s).
(3.8)
Note in particular that
ddst(pvpu)=trr2γABpApB0.
(3.9)

Proposition 3.3
(Null geodesics in FLRW). Let BP|Σ1 be a compact subset of the mass shell over Σ1 and let γ:[1,)M be a future directed null geodesic such that (γ(1),γ̇(1))B and let v0 be sufficiently large. There exist 0 < c < 1 < C, L0 ≥ 0, and s0 ∈ [1, ∞) such that v(γ(s0)) = v0 and the components of the tangent vector to γ satisfy, for all ss0,
γABpA(s)pB(s)12L0t(s)r(s)2,0pu(s)Cpv(s0)t(s)r(s)2,cpv(s0)t(s)pv(s)Cpv(s0)t(s).
(3.10)
Moreover, there exists retarded times U0 < 0 < U1 such that the u component of γ satisfies
U0u(s)U1,for all s1.
(3.11)

Proof.

Provided v0 is suitably large, the existence of s0 ∈ [1, ∞) such that v(γ(s0)) = v0 follows from the compactness of B.

Consider now the bounds (3.10). For simplicity, suppose that r(γ(1)) ≠ 0 (otherwise replace 1 in the following with 1 + ϵ for some ϵ > 0). By the conservation of angular momentum (2.20), one has the conservation law
L(s)2t(s)2r(s)4γABpA(s)pB(s)=t(1)2r(1)4γABpA(1)pB(1)for all s1,
(3.12)
from which the first of (3.10) trivially follows.
For the remaining inequalities of (3.10), suppose first that γ is radial, so that the conserved quantity L(s) vanishes. It follows from the mass shell relation (3.3) that either pu(1) = 0 or pv(1) = 0. The geodesic Eqs. (3.6) and (3.7) take the form
ddstpu=0,ddstpv=0,
and so, as long as the geodesic remains away from the centre where the coordinate system breaks down, one has t(s)pu(s) = pu(1) and t(s)pv(s) = pv(1). If pu(1) = 0 then, by (3.8), ṙ(s)>0 for all s ≥ 1 and so the solution remains away from the centre. The bounds (3.10) then trivially hold. If pv(1) = 0 then Eq. (3.8) imply
t(s)=32(s1)pu(1)+123,r(s)=r(1)1spu(1)32(s1)pu(1)+123ds.
It follows that the there is a time s* ≥ 1 at which the geodesic hits the centre, i.e., a time at which r(s*) = 0. The geodesic then becomes outgoing, with
pu(s)=0 for s>s*,andlimss*pv(s)=limss*pu(s)=pu(1)t(s*).
The bounds (3.10) then again trivially hold, as before.
Suppose now that γ is not radial, so that L ≠ 0, where L is defined by (3.12). Recall [see (3.9)] that
ddst(pvpu)=L2tr30.
If pv(1) ≥ pu(1) then it follows that pv(s) ≥ pu(s) for all s ≥ 1 and, by the latter of (3.8), r(s) is non-decreasing in s and the geodesic remains away from the centre. If pv(1) < pu(1) then, since 0ϵr3dr= for all ϵ > 0, it follows that there exists a time s* ≥ 1 such that pv(s*) = pu(s*), and moreover r(s) > 0 for all 1 ≤ ss*. As above r(s) is then non-decreasing for ss* and so the geodesic again remains away from the centre. Using now the mass shell relation (3.3) and the Eqs. (3.6) and (3.7),
ddsr2pupv=4r(pu)2pv0,
and so (defining s* = 1 if pv(1) ≥ pu(1) and integrating from s = s*),
r(s)2pu(s)r(s*)2pv(s),for all ss*.
(3.13)
Consider the latter of (3.10). It follows from (3.7) and (3.8) that
ddstpv=t2rr2γABpApB=2tr3t2r4γABpApB4=2t(s)r(s)3t(s0)2r(s0)2pu(s0)pv(s0),
(3.14)
The lower bound follows from the sign of the right hand side of (3.14). For the upper bound, (3.13) guarantees that, if v0 is sufficiently large (and hence s0s* is sufficiently large), then pv(s0) ≥ 2pu(s0). Equation (3.9) implies that
t(s)pv(s)pu(s)t(s0)pv(s0)pu(s0)for all ss0,
(3.15)
and so Eq. (3.14) implies that
t(s)pv(s)t(s0)pv(s0)s0s2t(s)r(s)3dst(s0)2r(s0)2pu(s0)pv(s0)r(s0)r(s)2r3drt(s0)r(s0)2pu(s0)pv(s0)pv(s0)pu(s0)4t(s0)pv(s0),
and the upper bound follows. The estimate for pu(s) in (3.10) finally follows from returning to (3.13).
Consider now (3.11). The existence of U0 follows trivially from the fact that
u̇(s)=pu(s)0.
The existence of U1 follows from the fact that,
u(s)u(s0)=s0spu(s)dsr(s0)r(s)Ctr2pv(s0)pvpudr2Ct(s0)r(s0),
by (3.10) and (3.15).□

The proof of Theorem 3.2 can now be given.

Proof of Theorem 3.2.
Consider Proposition 3.3 with B=supp(f1). The inclusion (3.5) follows immediately from (3.11). Let v0 be sufficiently large so that R(U1, v0) > 0. The inclusion (3.5) in particular implies that,
crt12Cr,in supp(f){vv0},
and so (3.10) imply that
γABpApB12L0t2,0puCpvt,0pvCt,ct2R2Ct2in supp(f){vv0}.
Hence, in view of the expression (2.10),
Tuu=π2Ω4R20L00Ct1fpvdpvLdLf1Lt2.
Similarly, in view of (2.11) and (2.12),
Tuv=π2Ω4R20L00Ct1fpudpvLdLf1Lt3,Tvv=π2Ω4R20L00Ct1f(pu)2pvdpvLdLf1Lt4.

Remark 3.4
(Massless Vlasov equation on Minkowski space). Consider Minkowski space (R3+1,m), where m=dt2+(dx1)2+(dx2)2+(dx3)2. In standard double null coordinates u = tr, v = t + r, r=(x1)2+(x2)2+(x3)2, with respect to which the metric takes the form
m=dudv+r2γ,
the components of the energy momentum tensor (1.5) for solutions of the massless Vlasov equation arising from compactly supported initial data decay with the rates
Tuuf1Lt2,Tuvf1Lt4,Tvvf1Lt6,
for t ≥ 1.25  Compare with the rates (3.4) which, with respect to the double null frame e3=t12u, e4=t12v normalised so that g(e3, e4) = −2, take the form
T33f1Lt3,T34f1Lt4,T44f1Lt5.

Remark 3.5
(Massless Vlasov equation on FLRW with T3 topology). The FLRW metric g, along with the function f, defined by (1.7) and (1.8), also define a solution of the Einstein–massless Vlasov system on the manifold (0,)×T3. The components of the energy momentum tensor (1.5) for solutions of the massless Vlasov equation arising from compactly supported initial data on such a background decay with the rates
Tttf|Σ1t2,i=13Titf|Σ1t32,i,j=13Tijf|Σ1t,
for all t, for a suitable norm ‖ ⋅‖. Compare again with the rates (3.4).

Theorem 1.4 can be more precisely stated as follows. Recall that initial data for the Einstein–massless Vlasov system consists of (Σ,ḡ,K,f1) (see Sec. II E). Let ‖ ⋅‖ be a norm in which local well posedness and Cauchy stability holds for the Einstein–massless Vlasov system (1.2)(1.5) in spherical symmetry. The norm on (Σ1,g|Σ1),
(ḡ,K,f1)=ḡCk+1+KCk+f1Ck,
suffices, for example, for k suitably large.

Theorem 4.1
(Future stability of FLRW). Let μ: [0, ∞) → [0, ∞) be a smooth function, decaying suitably so that R3|p|μ(|p|2)dp is finite, such that μ ≢ 0. Let g and f be defined by (1.7), (1.8). Let (Σ1,ḡ,K,f1) be a spherically symmetric initial data set for the Einstein–massless Vlasov system on Σ1={t=1}M such that f1f|Σ1 is compactly supported on P|Σ1 and ḡg|Σ1 and K12tg|Σ1 are compactly supported on Σ1. There exists ɛ0 > 0 such that, if
ḡg|Σ1,K12tg|Σ1,f1f|Σ1ε0,
then the resulting maximal developement of Theorem 2.2 is future geodesically complete and there exists U0 < 0 < U1 such that, in an appropriately normalised double null gauge, the associated (Q,Ω2,R,f) satisfies the estimates
Ω4R2TuuΩ4R2T°uuε0t6,Ω2R4TuvΩ2R4T°uvε0t8,R6TvvR6T°vvε0t10,
(4.1)
|Ω24t|ε0t,|Rt12r|ε0,|λv|ε0,|ν+u|ε0,
(4.2)
in the region {U0 < u < U1}, where t=14(v+u)2 and r = vu. Moreover,
π(supp(ff)){U0uU1},
where π:PM is the natural projection, and the solution is isometric to the spatially homogeneous FLRW solution (1.7) and (1.8) in the region {uU1}.

Note that, since the estimates (4.2) imply that Ω4R2t4, Ω2R4t5 and R6t6 in supp(ff) ∩ {vv0}, the rates (4.1) should be viewed as exactly the rates of Theorem 3.2.

The remainder of the paper is concerned with the Proof of Theorem 4.1.

This section concerns the Proof of Theorem 4.1. In Sec. V A the proof is reduced so that the main region of consideration arises from an appropriate characteristic initial value problem. In Sec. V B 1 more convenient renormalisation of the equations of Sec. II C are given. Section V C concerns Theorem 5.1, which forms the main content of the Proof of Theorem 4.1. In Sec. V D the Proof of Theorem 4.1 is completed.

Recall the setting of Theorem 4.1. Let (M,g) denote the unique maximal development of (Σ1,ḡ,K,f1), and let (Q,Ω2,R,f) denote the associated maximal solution of the reduced system (2.13)(2.19) (see Sec. II C), in which the residual gauge freedom, discussed in Remark 2.1, is normalised so that λ|{t=1} = λ|{t=1} and ν|{t=1}∪Γ = ν|{t=1}∪Γ.

Consider some v0 > U1, where U1 > 0 is to be determined later. Given δ0 > 0, L0 ≥ 0, define the set
Uδ0,L0,v0={(x,p)Pv(x)=v0,2v0u(x)v0δ0,0pu2pv,L(x,p)L0}P,
so that, in particular, rδ0 and t ≥ 1 in Uδ0,L0,v0.
By the domain of dependence property, in the region {uU0}, for some appropriate U0 < 0, the solution coincides with the FLRW solution (1.7) and (1.8). In particular,
Ω2=4t,vlogΩ2=t12,R=t12r,λ=v,ν=u,on {u=U0}.
(5.1)
Let t0=t(U1,v0)=14(v0+U1)2. Let L0 > 0 be such that
R4γABpApB(L0)2insupp(f1f|Σ1),
and let δ0 > 0 be such that r(U1, v0) = v0U1 > δ0. By Cauchy stability (for both the Einstein–massless Vlasov system (1.2)(1.5) and the geodesic equations) the solution exists up to {t = t0} and, by Proposition 3.3, for any future-maximal null geodesic γ:[1,S)M such that (γ(1),γ̇(1))supp(f1), there exists s0 ∈ [1, ∞) such that
(γ(s0),γ̇(s0))Uδ0,L0,v0,
if v0 is sufficiently large and ɛ0 > 0 is sufficiently small. Moreover
|Ω24t|ε0,|ulogΩ2t12|ε0,|Rt12r|ε0,|λv|ε0,|ν+u|ε0,on {v=v0}.
(5.2)

The remainder of the proof will then be concerned with the region {vv0} ∩ {uU0}.

Rather than considering the Eqs. (2.13)(2.16) directly, in the Proof of Theorem 4.1 it is more convenient to consider the following renormalisation of the equations with their respective FLRW quantities,
uv(R2R2)=Ω2Ω22+(R2R2)Tuv+R2(TuvT°uv),
(5.3)
uvlogΩ2Ω2=Ω22R2Ω22R2+2R2νλ2R2νλ2(TuvT°uv),
(5.4)
u(Ω2νΩ2ν)=12RΩ2T°uu12RΩ2Tuu,
(5.5)
v(Ω2λΩ2λ)=12RΩ2T°vv12RΩ2Tvv,
(5.6)
where the solution is compared with the FLRW solution using the double null gauge of Sec. III A, as discussed at the end of Sec. III A.
For any null geodesic γ:[1,S)M, write the tangent vector to γ as
γ̇(s)=pu(s)u+pv(s)v+pA(s)θA.
The geodesic equations for pu(s) and pv(s) take the form
ṗu(s)+ulogΩ2pu(s)2+2λRΩ2γABpApB=0,
(5.7)
ṗv(s)+vlogΩ2pv(s)2+2νRΩ2γABpApB=0.
(5.8)
The following renormalisation of (5.8),
dds(tpv)=tt12vlogΩ2pv2+(RR)2t12(νν)t12RΩ2γABpApB+2t32RΩ2γABpApB,
(5.9)
will also be used. Equation (5.9) is obtained from (5.8) using the mass shell relation (2.17) and the fact that ṫ=t12(pv+pu).

The main content of the Proof of Theorem 4.1 is contained in the following bootstrap theorem.

Given T > t0 and U1 > 0, define the region
QU1,T={(u,v)QU0uU1,vv0,andt0t(u,v)<T}.
(5.10)

Theorem 5.1
(Bootstrap theorem). Suppose that T > t0 is such that the maximal development exists up to time T — in the sense that (u(t,r),v(t,r))Q for all t0t < T, r ≥ 0 — the difference ff is supported in the region
supp(π(ff)){U0uU1},
(5.11)
and, for some fixed ɛ > 0 and U1 > 0, the solution moreover satisfies, for all (u,v)QU1,T,
Ω4R2TuuΩ4R2T°uu<εt6,Ω2R4TuvΩ2R4T°uv<εt8,R6TvvR6T°vv<εt10,
(5.12)
|Ω24t|<εt,vlogΩ21t12<εt,ulogΩ21t12<ε,
(5.13)
|Rt12r|<ε,|R2tr2|<εt,|λv|<ε,|ν+u|<ε.
(5.14)
Then, if U1 is suitably large and ɛ is suitably small, there exists a constant C > 0 (independent of ɛ and U1) such that the inequalities (5.12)(5.14) hold in QU1,T with ɛ replaced by 0 and (5.11) holds with U1 replaced by U1/2. Moreover, the solution is isometric to the spatially homogeneous FLRW solution (1.7) and (1.8) in the region {t0t(u, v) < T} ∩ {uU1}.

Proof. The proof is divided into several steps. First, the size of the support of ff is estimated. These estimates are then used to obtain estimates on the components of the energy momentum tensor. The metric quantities are then estimated in the region {U0uU1}. Finally, the solution is shown to be isometric to the spatially homogeneous FLRW solution (1.7) and (1.8) in the region {uU1}.

Throughout the proof the notation AB is used when there exists a constant K, which may depend on U1U0, such that AKB. Constants 0 < c ≪ 1 ≪ C are always independent of U1U0.

Estimates for the support of ff: Let δ0, L0 and v0 be as in Sec. V A, and recall that, for each null geodesic γ:[1,S)M emanating from supp(f1), if S is sufficiently large there exists a time s0 ∈ [1, S) such that (γ(s0),γ̇(s0))Uδ0,L0,v0.

For such a null geodesic, write
γ̇(s)=pu(s)u+pv(s)v+pA(s)θA.
By the bootstrap assumption (5.11), the u coordinate of γ satisfies
U0u(s)U1,
for all s0s < S. Recall moreover that v0 > U1 and so R(v0, U1) > 0. It follows in particular that
t12rR12v,along γ.
(5.15)
There exist constants 0 < c ≪ 1 ≪ C such that the components of the tangent vector γ̇ of any such null geodesic γ:[1,S)M emanating from supp(f1) satisfy
γABpA(s)pB(s)122L0t(s)r(s)2,0pu(s)Cpv(s0)t(s)r(s)2,cpv(s0)t(s)pv(s)Cpv(s0)t(s),
(5.16)
for all ss0. Indeed, let γ:[s0,S)M be a null geodesic such that (γ(s0),γ̇(s0))Uδ0,L0,v0. The first of (5.16) follows from the conservation of angular momentum (2.20), which implies that
R(s)4γABpA(s)pB(s)=R(s0)4γABpA(s0)pB(s0)(L0)2,
(5.17)
along with the fact that t2r4 ≤ 2R4, which follows from the bootstrap assumption (5.14).
It follows from the geodesic equation (5.9) that
ddstpvexps0sh(s)pv(s)ds=exps0sh(s)pv(s)ds2t32RΩ2γABpApB,
(5.18)
where
h(s)=t12vlogΩ2+(RR)2t12(νν)t12RΩ2γABpApB(pv)2.
In order to obtain the lower bound of the latter of (5.16), consider the set
E=s[s0,S)(tpv)(s)12(tpv)(s0) for all s0ss[s0,S).
The set E is manifestly a closed, connected, non-empty subset of [s0, S). Moreover, if sE, then
s0s|h(s)|pv(s)ds=v0v(s)|h|dvCv0v(s)εv2dvCε,
(5.19)
by the fact that
R(s)4γABpA(s)pB(s)=R(s0)4γABpA(s0)pB(s0)=R(s0)2Ω(s0)2pu(s0)pv(s0),
(5.20)
along with the bootstrap assumptions (5.13) and (5.14) and the equivalence (5.15). Since
ddstpvexps0sh(s)pv(s)ds0,
it follows that
(tpv)(s0)eCε(tpv)(s),
and thus, if ɛ is sufficiently small, the set E[s0,S) is open, and hence equal to [s0, S).
For the upper bound of the latter of (5.16), recall Eq. (5.18) and note that
s0s2t32RΩ2γABpApBds=v0v(s)2t52R3Ω2R4γABpApBtpvdvCv0v(s)pv(s0)v3dvCpv(s0),
by (5.15) and (5.20) and the lower bound for (tpv)(s). The upper bound for (tpv)(s) then follows from the estimate (5.19), after integrating Eq. (5.18).
The estimate for pu(s) in (5.16) finally follows from the mass shell relation (2.17), which gives
Ω(s)2t(s)pu(s)=R4γABpApB(s)t(s)pv(s),
along with the former and latter of (5.16).
It in particular follows from (5.16) that
γABpApB12CL0t2,0puCpvt,0pvCt,in supp(ff).
(5.21)
Consider now the support property (5.11). For any null geodesic γ:[s0,S)M as above, it follows from (5.16) that,
u(s)u(s0)=s0spu(s)dsv(s0)v(s)Ctv2pv(s0)pvdvCt(s0)v0.
Thus, since t(s0), v0 ≥ 1,
supp(π(ff)){U0uC},
(5.22)
and so (5.11) in fact holds with U1 replaced by U1/2, provided that U1 ≥ 2C.
Estimates for the energy momentum tensor components: Recall the expressions (2.10)(2.12) for the components of the energy momentum tensor, along with the expression (3.2) for the components of the energy momentum tensor of f. By (5.22), for any (u, v, pv, L),
|f(u,v,pv,L)f(u,v,pv,L)|ε01{U0uU1/2}.
(5.23)
Define
L2=R4γABpApB,
and note that R4L2=R4L2. It follows from the support property (5.21) that
Ω4R2TuuΩ4R2T°uuπ20L00Ct1|ff|pvdpvR4LdLε0t61{U0uU1/2}.
Noting that Ω2R2pu=Ω2R2p̊u, it similarly follows that
Ω2R4TuvΩ2R4T°uvπ20L00Ct1|ff|Ω2puR2dpvR4LdLε0t81{U0uU1/2},
and
R6TvvR6T°vvπ20L00Ct1|ff|Ω4(pu)2R4pvdpvR4LdLε0t101{U0uU1/2}.
Thus
Tuu1tε0t21{U0uU1/2}+t2|RR|+t2|Ω2Ω2|,
(5.24)
Tuv12tε0t31{U0uU1/2}+t2|RR|+t2|Ω2Ω2|,
(5.25)
Tvv1tε0t41{U0uU1/2}+t2|RR|.
(5.26)
Estimates for the metric quantities in the region uU1: For the metric quantities, consider the region U0uU1 and recall that
t12rR12v,for U0uU1,vv0.
Let u and v be such that U0uU1, vv0 and t0t(u, v) < T.
Note that
12v(R2R2)=Rλλ+(RR)λ.
Equation (5.3) and the estimate (5.25) then imply that
uR(λλ+(RR)λ|Ω2Ω2|+|RR|+ε0v2.
Thus, integrating and using the boundary condition (5.1),
|λλ|(u,v)v1|RR|(u,v)+v2U0u|Ω2Ω2|(u,v)+|RR|(u,v)du+ε0v4.
(5.27)
Similarly, since
νν=Ω2Ω2νΩ2νν(1Ω2Ω2),
it follows from Eq. (5.5) and the estimate (5.25) and the boundary condition (5.1) that
|νν|(u,v)v2|Ω2Ω2|(u,v)+v2U0u|Ω2Ω2|(u,v)+|RR|(u,v)du+ε0v2.
(5.28)
Moreover,
u(RR)=νν,
and so (5.28) and the Grönwall inequality imply that
|RR|(u,v)v2U0u|Ω2Ω2|(u,v)du+ε0v2.
(5.29)
Returning to (5.27) and (5.28), one then has
|λλ|(u,v)v2U0u|Ω2Ω2|(u,v)du+ε0v3,
(5.30)
and
|νν|(u,v)v2|Ω2Ω2|(u,v)+v2U0u|Ω2Ω2|(u,v)du+ε0v2.
(5.31)
Now, by Eq. (5.4) and the estimate (5.25)
uvlogΩ2Ω2v4|Ω2Ω2|+v4|RR|+v3|νν|+v4|λλ|+ε0v4,
(5.32)
Integrating in u and v, using the boundary condition (5.1) and (5.2), and inserting (5.29)(5.31), it follows that
logΩ2Ω2ε0+v0vU0u(v)4|Ω2Ω2|(v,u)dudv
(5.33)
Now since, for ɛ suitably small,
|Ω2Ω2|Ω2logΩ2Ω2v2logΩ2Ω2for U0uU1,
it follows from the Grönwall inequality, after integrating (5.33) once more in u, that
|Ω2Ω2|ε0v2.
Returning to (5.29)(5.31), it then follows that
|RR|ε0,|λλ|ε0,|νν|ε0.
Similarly, returning to (5.32),
vlogΩ21t12ε0v2,ulogΩ21t12ε0.
The region uU1: In order to see that the solution is isometric to the spatially homogeneous FLRW solution (1.7) and (1.8) in the region uU1, recall first (5.23), from which it follows that
R6Tvv(U1,v)=R6T°vv(U1,v)=t(U1,v)4r(U1,v)6for all vV1,
where V1 is such that R(U1, V1) = 0. Recall the residual gauge freedom of Remark 2.1 and consider the following change of gauge. Define increasing functions ṽ(v) and ũ(u) as solutions of
ṽ(V1)=U1,ṽ(v)=λ(U1,v)λ1(U1,ṽ(v)),for vV1,
and
ũ(u)=ufor uU1,ũ(u)=1ṽ(V1)Ω2(u,V1)Ω2(ũ(u),V1)for u>U1,
respectively. It follows that {u=U1}={ũ=U1}, and moreover the relations (2.3) and (2.6) imply that the solution in the resulting gauge satisfies
R(U1,U1)=0,Ω2(U1,U1)=Ω2(U1,U1)vR(U1,v)=vR(U1,v),for all vU1.
The Raychaudhuri Eq. (2.16), restricted to u = U1, takes the form
v(Ω2vR)(U1,v)=12R7Ω2R6T°vv(U1,v).
It then follows that
R(U1,v)=R(U1,v),Ω2(U1,v)=Ω2(U1,v)for all vU1.
By uniqueness for the characteristic initial value problem, it follows that R = R and Ω2=Ω2 for all uU1. □

Recall the set QU1,TQ defined by (5.10). Consider the set A[t0,) of times T ∈ [t0, ∞) such that the solution exists up to time T — in the sense that (u(t,r),v(t,r))Q for all t0t < T, r ≥ 0 — and moreover satisfies the estimates (5.12)(5.14) for all (u,v)QU1,T for some ɛ > 0 sufficiently small so that the conclusion of Theorem 5.1 holds. The set A[t0,) is manifestly connected and open, and is moreover non-empty by local existence and Cauchy stability, provided ɛ0 is sufficiently small. By Theorem 2.3, the bounds (5.13)(5.14) imply that (u(T,r),v(T,r))Q for all r ≥ 0 such that u(T, r) ≤ U1. By Theorem 5.1, the solution is isometric to the spatially homogeneous FLRW solution (1.7) and (1.8) in the region {t0t(u, v) < T}∩{uU1} and thus (u(T,r),v(T,r))Q for all r ≥ 0. Theorem 5.1 moreover guarantees that the bounds (5.12)(5.14) cannot be saturated, and thus the set A is moreover closed, and hence equal to [t0, ∞).

The future geodesic completeness follows from the bounds on the Christoffel symbols arising from (5.13) and (5.14) (as, for example in Sec. 16 of Ref. 38).

I acknowledge support through Royal Society Tata University Research Fellowship Grant No. URF∖R1∖191409. I am grateful to G. Fournodavlos, H. Masaood, and J. Speck for helpful discussions.

The author has no conflicts to disclose.

Martin Taylor: Conceptualization (equal); Investigation (equal); Methodology (equal); Writing – original draft (equal).

Data sharing is not applicable to this article as no new data were created or analyzed in this study.

1.
S. W.
Hawking
and
G. F.
Ellis
,
The Large Scale Structure of Space-Time
(
Cambridge University Press
,
2023
).
2.
J.
Ehlers
,
P.
Geren
, and
R. K.
Sachs
, “
Isotropic solutions of the Einstein–Liouville equations
,”
J. Math. Phys.
9
(
9
),
1344
1349
(
1968
).
3.
M.
Dafermos
, “
A note on the collapse of small data self-gravitating massless collisionless matter
,”
J. Hyperbolic Differ. Equ.
03
(
04
),
589
598
(
2006
).
4.
M.
Dafermos
and
A. D.
Rendall
, “
Strong cosmic censorship for surface-symmetric cosmological spacetimes with collisionless matter
,”
Commun. Pure Appl. Math.
69
(
5
),
815
908
(
2016
).
5.
J.
Natário
and
F.
Rossetti
, “
Explicit formulas and decay rates for the solution of the wave equation in cosmological spacetimes
,”
J. Math. Phys.
64
(
3
),
2023
(
2023
).
6.
J.
Speck
, “
The stabilizing effect of spacetime expansion on relativistic fluids with sharp results for the radiation equation of state
,”
Arch. Ration. Mech. Anal.
210
,
535
579
(
2013
).
7.
E. M.
Lifshitz
, “
On the gravitational stability of the expanding universe
,”
J. Phys.
10
,
116
129
(
1946
).
8.
G.
Börner
,
The Early Universe: Facts and Fiction
(
Springer Science & Business Media
,
2013
).
9.
H.
Ringström
,
On the Topology and Future Stability of the Universe
(
Oxford University Press
,
Oxford
,
2013
).
10.
L.
Andersson
and
D.
Fajman
, “
Nonlinear stability of the Milne model with matter
,”
Commun. Math. Phys.
378
(
1
),
261
298
(
2020
).
11.
H.
Barzegar
and
D.
Fajman
, “
Stable cosmologies with collisionless charged matter
,”
J. Hyperbolic Differ. Equ.
19
(
04
),
587
634
(
2022
).
12.
D.
Fajman
, “
The nonvacuum Einstein flow on surfaces of negative curvature and nonlinear stability
,”
Commun. Math. Phys.
353
(
2
),
905
961
(
2017
).
13.
D.
Fajman
,
M.
Ofner
, and
Z.
Wyatt
, “
Slowly expanding stable dust spacetimes
,” arXiv:2107.00457 (
2021
).
14.
M.
Hadžić
and
J.
Speck
, “
The global future stability of the FLRW solutions to the dust-Einstein system with a positive cosmological constant
,”
J. Hyperbolic Differ. Equ.
12
(
01
),
87
188
(
2015
).
15.
T. A.
Oliynyk
, “
Future stability of the FLRW fluid solutions in the presence of a positive cosmological constant
,”
Commun. Math. Phys.
346
,
293
312
(
2016
).
16.
I.
Rodnianski
and
J.
Speck
, “
The nonlinear future stability of the FLRW family of solutions to the irrotational Euler–Einstein system with a positive cosmological constant
,”
J. Eur. Math. Soc.
15
(
6
),
2369
2462
(
2013
).
17.
J.
Speck
, “
The nonlinear future stability of the FLRW family of solutions to the Euler–Einstein system with a positive cosmological constant
,”
Selecta Math.
18
(
3
),
633
715
(
2012
).
18.
L.
Andersson
and
V.
Moncrief
, “
Einstein spaces as attractors for the Einstein flow
,”
J. Differ. Geom.
89
(
1
),
1
47
(
2011
).
19.
D.
Fajman
and
Z.
Wyatt
, “
Attractors of the Einstein–Klein–Gordon system
,”
Commun. Partial Differ. Equ.
46
(
1
),
1
30
(
2021
).
20.
G.
Fournodavlos
, “
Future dynamics of FLRW for the massless-scalar field system with positive cosmological constant
,”
J. Math. Phys.
63
(
3
),
032502
(
2022
).
21.
P.
Mondal
, “
The nonlinear stability of n + 1 dimensional FLRW spacetimes
,” arXiv:2203.04785 (
2022
).
22.
G.
Rein
and
A. D.
Rendall
, “
Global existence of solutions of the spherically symmetric Vlasov–Einstein system with small initial data
,”
Commun. Math. Phys.
150
,
561
583
(
1992
).
23.
D.
Fajman
,
J.
Joudioux
, and
J.
Smulevici
, “
The stability of the Minkowski space for theEinstein–Vlasov system
,”
Anal. PDE
14
(
2
),
425
531
(
2021
).
24.
H.
Lindblad
and
M.
Taylor
, “
Global stability of Minkowski space for the Einstein–Vlasov system in the harmonic gauge
,”
Arch. Ration. Mech. Anal.
235
(
1
),
517
633
(
2020
).
25.
M.
Taylor
, “
The global nonlinear stability of Minkowski space for the massless Einstein–Vlasov system
,”
Ann. PDE
3
(
1
),
9
(
2017
).
26.
G.
Moschidis
, “
A proof of the instability of AdS for the Einstein-massless Vlasov system
,”
Invent. Math.
231
,
467
672
(
2022
).
27.
J.
Bedrossian
,
N.
Masmoudi
, and
C.
Mouhot
, “
Linearized wave-damping structure of Vlasov-Poisson in R3
,”
SIAM J. Math. Anal.
54
(
4
),
4379
4406
(
2022
).
28.
D.
Han-Kwan
,
T. T.
Nguyen
, and
F.
Rousset
, “
On the linearized Vlasov–Poisson system on the whole space around stable homogeneous equilibria
,”
Commun. Math. Phys.
387
,
1405
1440
(
2021
).
29.
A.
Ionescu
,
P.
Benoit
,
X.
Wang
, and
K.
Widmayer
, “
Nonlinear Landau damping for the Vlasov-Poisson system in R3: The Poisson equilibrium
,”
Ann. PDE
10
,
2
(
2024
).
30.
J.
Bedrossian
,
N.
Masmoudi
, and
C.
Mouhot
, “
Landau damping: Paraproducts and Gevrey regularity
,”
Ann. PDE
2
(
1
),
4
(
2016
).
31.
E.
Grenier
,
T. T.
Nguyen
, and
I.
Rodnianski
, “
Landau damping for analytic and Gevrey data
,” arXiv:2004.05979 (
2020
).
32.
C.
Mouhot
and
C.
Villani
, “
On Landau damping
,”
Acta Math.
207
,
29
201
(
2011
).
33.
M. K.-H.
Kiessling
, “
The ‘Jeans swindle’: A true story—Mathematically speaking
,”
Adv. Appl. Math.
31
(
1
),
132
149
(
2003
).
34.
G.
Moschidis
, “
The characteristic initial-boundary value problem for the Einstein–massless Vlasov system in spherical symmetry
,” arXiv:1812.04274 (
2018
).
35.
Y.
Choquet-Bruhat
, “
Problème de Cauchy pour le système intégro-différentiel d’Einstein–Liouville
,”
Ann.Inst. Fourier
21
,
181
201
(
1971
).
36.
Y.
Choquet-Bruhat
and
R.
Geroch
, “
Global aspects of the Cauchy problem in general relativity
,”
Commun. Math. Phys.
14
,
329
335
(
1969
).
37.
M.
Dafermos
and
A. D.
Rendall
, “
An extension principle for the Einstein–Vlasov system in spherical symmetry
,”
Ann. Henri Poincaré
6
,
1137
1155
(
2005
).
38.
H.
Lindblad
and
I.
Rodnianski
, “
Global existence for the Einstein vacuum equations in wave coordinates
,”
Commun. Math. Phys.
256
,
43
110
(
2005
).
Published open access through an agreement with JISC Collections