# Future stability of expanding spatially homogeneous FLRW solutions of the
spherically symmetric Einstein–massless Vlasov system with spatial topology
$\mathbb{R}^{3}$

Martin Taylor  Imperial College London, Department of Mathematics, South
Kensington Campus, London SW7 2AZ, United Kingdom<EMAIL_ADDRESS>
(June 30, 2023)

###### Abstract

Spatially homogeneous FLRW solutions constitute an infinite dimensional family
of explicit solutions of the Einstein–massless Vlasov system with vanishing
cosmological constant. Each member expands towards 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 $\mathbb{R}^{3}$. 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.

###### Contents

  1. 1 Introduction
    1. 1.1 The Einstein–massless Vlasov system
    2. 1.2 The spatially homogeneous FLRW family of solutions
    3. 1.3 First version of the main theorem
    4. 1.4 Previous works
    5. 1.5 Outline of the paper
  2. 2 The spherically symmetric Einstein–massless Vlasov system
    1. 2.1 Spherical symmetry
    2. 2.2 Double null gauge
    3. 2.3 The spherically symmetric Einstein–massless Vlasov system in double null gauge
    4. 2.4 The functions $t$ and $r$
    5. 2.5 The Cauchy problem
  3. 3 The FLRW spacetimes
    1. 3.1 The FLRW metrics in double null gauge
    2. 3.2 The massless Vlasov equation on an FLRW background
  4. 4 The main theorem
  5. 5 The proof of the main theorem
    1. 5.1 Cauchy stability and the domain of dependence property
    2. 5.2 The renormalised equations
    3. 5.3 The bootstrap theorem
    4. 5.4 The conclusion of the proof of the main theorem

##  1 Introduction

Standard homogeneous isotropic cosmological models in general relativity are
described by the Friedmann–Lemaître–Robertson–Walker (FLRW) spacetimes

$\mathcal{M}=I\times\Sigma,\qquad g=-dt^{2}+a(t)^{2}g_{\Sigma},$ (1.1)

where $I\subset\mathbb{R}$ is an open interval, $(\Sigma,g_{\Sigma})$ is a
constant curvature manifold, and $a\colon I\to(0,\infty)$ is an appropriate
_scale factor_. See Section 5.3 of [21] for more on FLRW spacetimes. This
article concerns radiation filled FLRW cosmologies in which the constant
curvature manifold is Euclidean space,
$(\Sigma,g_{\Sigma})=(\mathbb{R}^{3},g_{\mathrm{Eucl}})$, with radiation
described by spatially homogenous solutions of the massless Vlasov equation,
and their stability properties.

###  1.1 The Einstein–massless Vlasov system

Consider a $3+1$ dimensional Lorentzian manifold $(\mathcal{M},g)$ and let

$P=\\{(t,x,p)\in T\mathcal{M}\mid g_{(t,x)}(p,p)=0\\}\subset T\mathcal{M},$

denote the _mass shell_ of $(\mathcal{M},g)$. Consider some local coordinates
$\\{t,x^{1},x^{2},x^{3}\\}$ on $\mathcal{M}$, and let $\\{t,x^{i},p^{\mu}\\}$
denote the corresponding conjugate coordinate system for $T\mathcal{M}$, i.e.
$(t,x^{i},p^{\mu})$ describes the point $p^{\mu}\partial_{x^{\mu}}|_{(t,x)}\in
T\mathcal{M}$. The massless Vlasov equation on $(\mathcal{M},g)$ takes the
form

$p^{0}\partial_{t}f+p^{i}\partial_{x^{i}}f-p^{\mu}p^{\nu}\Gamma_{\mu\nu}^{i}\partial_{p^{i}}f=0,$
(1.2)

where $p^{0}$ is defined by the mass shell relation

$g_{\mu\nu}p^{\mu}p^{\nu}=0.$ (1.3)

The Einstein–massless Vlasov system consists of equation (1.2) coupled to the
Einstein equations

$Ric(g)_{\mu\nu}-\frac{1}{2}R(g)g_{\mu\nu}=T_{\mu\nu},$ (1.4)

where the energy momentum tensor takes the form

$T^{\mu\nu}(t,x)=\int_{P_{(t,x)}}f(t,x,p)p^{\mu}p^{\nu}\frac{\sqrt{-\det
g}}{-p_{0}}dp^{1}dp^{2}dp^{3},$ (1.5)

where indices are raised and lowered with respect to the metric $g$ (so that,
for example, $p_{0}=g_{0\mu}p^{\mu}$). Here $t$ is also denoted $x^{0}$, Greek
indices, such as $\mu,\nu$, range over $0,1,2,4$, and lower case Latin
indices, such as $i,j,k$ range over $1,2,3$.

###  1.2 The spatially homogeneous FLRW family of solutions

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

$\mathcal{M}_{\circ}=(0,\infty)\times\mathbb{R}^{3},$ (1.6)

and consider the Euclidean metric

$g_{\mathrm{Eucl}}=(dx^{1})^{2}+(dx^{2})^{2}+(dx^{3})^{2},$

on $\mathbb{R}^{3}$. Given $x\in\mathbb{R}^{3}$ and $p=(p^{1},p^{2},p^{3})\in
T_{x}\mathbb{R}^{3}$, define

$|p|^{2}=|p|_{g_{\mathrm{Eucl}}}^{2}=(g_{\mathrm{Eucl}})_{ij}p^{i}p^{j}=(p^{1})^{2}+(p^{2})^{2}+(p^{3})^{2}.$

For any smooth, sufficiently decaying function $\mu:[0,\infty)\to[0,\infty)$
such that $\mu\not\equiv 0$, the metric $g_{\circ}$ and particle density
$f_{\circ}$ defined by

$g_{\circ}=-dt^{2}+a(t)^{2}\big{(}(dx^{1})^{2}+(dx^{2})^{2}+(dx^{3})^{2}\big{)},\qquad
f_{\circ}(t,x,p)=\mu(a(t)^{4}|p|^{2}),$ (1.7)

where

$a(t)=t^{\frac{1}{2}}\left(\frac{4\varrho}{3}\right)^{\frac{1}{4}},\qquad\varrho=\int|p|\mu(|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^{\prime\prime}(t)>0$ for
all $t$, _linear expansion_ if $a^{\prime\prime}(t)=0$ for all $t$, or
_decelerated expansion_ if $a^{\prime\prime}(t)<0$ for all $t$. Note that the
solutions (1.7)–(1.8) undergo decelerated expansion. Previous stability works
in cosmological settings have typically considered perturbations of spacetimes
undergoing linear or accelerated expansion. See Section 1.4 below.

######  Remark 1.2 (Anisotropic spatially homogeneous FLRW solutions).

The spatially homogeneous FLRW family of solutions (1.7)–(1.8) lie inside a
more general family of _anisotropic_ spatially homogeneous FLRW family. Given
any smooth, sufficiently decaying function $F:\mathbb{R}^{3}\to[0,\infty)$
such that $F\not\equiv 0$ and

$\int_{\mathbb{R}^{3}}F(p)p^{i}dp=0,\quad
i=1,2,3,\qquad\text{and}\qquad\int_{\mathbb{R}^{3}}F(p)\frac{p^{i}p^{j}}{|p|}dp=\frac{\varrho}{3}\delta^{ij},\quad
i,j=1,2,3,$

the metric and particle density

$g_{\circ}=-dt^{2}+a(t)^{2}\big{(}(dx^{1})^{2}+(dx^{2})^{2}+(dx^{3})^{2}\big{)},\qquad
f_{\circ}(t,x,p)=F(a(t)^{2}p^{1},a(t)^{2}p^{2},a(t)^{2}p^{3}),$

where

$a(t)=t^{\frac{1}{2}}\left(\frac{4\varrho}{3}\right)^{\frac{1}{4}},\qquad\varrho=\int|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 Section 2.1), consideration here is
restricted to the solutions of the form (1.7)–(1.8).

######  Remark 1.3 ($\mathbb{T}^{3}$ spatial topology).

Each (1.7)–(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,\infty)\times\mathbb{T}^{3}$. This article will only be
concerned with solutions on $(0,\infty)\times\mathbb{R}^{3}$. (See also Remark
3.5 for a further comment on the solutions on
$(0,\infty)\times\mathbb{T}^{3}$.)

The components of the energy momentum tensor of (1.7)–(1.8), with respect to
the Cartesian $(t,x^{1},x^{2},x^{3})$ coordinate system, take the form

$(T^{\circ})^{tt}(t,x)=\frac{3}{4t^{2}},\qquad(T^{\circ})_{i}^{j}(t,x)=\frac{1}{4t^{2}}\delta_{i}^{j},\qquad(T^{\circ})_{it}(t,x)=0,$
(1.9)

for $i,j=1,2,3$, where
$(T^{\circ})_{i}^{j}(t,x)=g_{ik}(t,x)(T^{\circ})^{jk}(t,x)$.

The spacetimes $(\mathcal{M}_{\circ},g_{\circ})$ expand from a spacelike
singularity at $t=0$ (with Kretchmann scalar
$|Rm_{\circ}|^{2}_{g_{\circ}}=\frac{3}{2t^{4}}$) and are future geodesically
complete. See the Penrose diagram of Figure 1. Such solutions feature in the
Ehlers–Geren–Sachs Theorem [12], 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.

###  1.3 First version of the main theorem

The main result of the present work concerns the future stability of the
family (1.7)–(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)–(1.8) for some
smooth, non-identically vanishing $\mu$ (decaying suitably so that $\varrho$
is defined). For all compactly supported spherically symmetric initial data
sufficiently close to that of (1.7)–(1.8) on $\Sigma_{1}=\\{t=1\\}$, the
resulting solution is future geodesically complete, is isometric to
(1.7)–(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 Section 4. See also
the Penrose diagram in Figure 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)–(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 [9], 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 [11] or Theorem 2.3 below).

######  Remark 1.5 (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.

######  Remark 1.6 (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_{\circ}$
(defined by (1.7)), is locally isometric to the FLRW spacetime
$(\mathcal{M}_{\circ},g_{\circ})$. 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 [12].

######  Remark 1.7 (Einstein–Euler with radiation equation of state).

The metric $g_{\circ}$, defined by the former of (1.7) with
$a(t)=t^{\frac{1}{2}}$, also arises as a solution of the Einstein–Euler
system, i.e. equation (1.4) with

$T^{\mu\nu}=(\rho+p)u^{\mu}u^{\nu}+pg^{\mu\nu},$

for unknowns $g$, four velocity $u^{\mu}$ — satisfying
$g_{\mu\nu}u^{\mu}u^{\nu}=-1$ — pressure $p$, and density $\rho$, with $p$ and
$\rho$ related by the _radiation equation of state_

$p=\frac{\rho}{3}.$

The fluid variables for the FLRW metric take the form, in Cartesian
coordinates,

$u=\partial_{t},\qquad\rho=\frac{3}{4t^{2}}.$ (1.10)

It is known that the solution (1.10) is unstable to formation of shock waves
for the Euler equations $\nabla_{\mu}T^{\mu\nu}=0$ on the fixed FLRW
background $(\mathcal{M}_{\circ},g_{\circ})$ (see Chapter 9 of [34]).

Note also the linear _Jeans instability_ for the Einstein–Euler system
linearised around FLRW [24] (see also the discussion in [6]).

###  1.4 Previous works

There have been previous works on the Einstein–Vlasov system in cosmological
settings for Einstein equations with a positive cosmological constant [31],
and also for perturbations of the vacuum Milne solution [1, 3, 13]. Note also
related works [19, 32, 33, 15] on the Einstein–Euler system. See also [2, 16,
17]. 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)–(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 [9, 30], and later to general perturbations [35,
26, 14]. The present work is based on the approach of [9]. _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 [28].

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 works [20, 5, 22] on
Landau damping, and also the celebrated works [29, 4, 18] on the torus
$\mathbb{T}^{3}$, 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 [23].

###  1.5 Outline of the paper

Section 2 contains preliminaries on the Einstein–massless Vlasov system in
double null gauge. Section 3 concerns further properties of the FLRW
spacetimes, introduced in Section 1.2. A double null gauge in FLRW is defined
in Section 3.1 and, as a precursor to the proof Theorem 1.4, properties of the
massless Vlasov equation (1.2)–(1.3) on the FLRW background spacetimes are
presented in Section 3.2. In Section 4 a more precise version of Theorem 1.4
is formulated, and in Section 5 its proof is given.

### Acknowledgements

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

##  2 The spherically symmetric Einstein–massless Vlasov system

This section concerns facts about the Einstein–massless Vlasov system in
spherical symmetry. In Section 2.1 the spherical symmetry assumption is
introduced. In Section 2.2 spherically symmetric double null gauges are
introduced, and the residual spherically symmetric double null freedom is
discussed. In Section 2.3 the Einstein–massless Vlasov system in a spherically
symmetric double null gauge is presented. In Section 2.4 functions $t$ and $r$
are introduced. Finally, in Section 2.5 the Cauchy problem is discussed and an
extension principle, which will be used in the proof of Theorem 1.4, is
stated.

###  2.1 Spherical symmetry

A solution $(\mathcal{M},g,f)$ of (1.2)–(1.5) is _spherically symmetric_ if
$SO(3)$ acts by isometry on $(\mathcal{M},g)$ and preserves $f$. More
precisely, it is assumed that there is a smooth isometric action
$\mathcal{O}\colon\mathcal{M}\times SO(3)\to\mathcal{M}$ on $(\mathcal{M},g)$
such that, for each $p\in\mathcal{M}$, the orbit of the action satisfies
either $Orb(p)\simeq S^{2}$, or $Orb(p)=\\{p\\}$ and, moreover, for each of
the generators $\Omega_{1},\Omega_{2},\Omega_{3}$ of the $SO(3)$ action, the
corresponding flows $\Phi^{i}\colon(0,2\pi)\times\mathcal{M}\to\mathcal{M}$
satisfy

$f(x,p)=f(\Phi^{i}_{s}(x),(\Phi^{i}_{s})_{*}p),\qquad\text{for all
}s\in(0,2\pi),\quad i=1,2,3,$ (2.1)

for all $(x,p)\in\mathcal{P}$. It is assumed that the centre

$\Gamma=\\{x\in\mathcal{M}\mid\mathcal{O}(x,\mathfrak{s})=x\text{ for all
}\mathfrak{s}\in SO(3)\\},$

consists of a single timelike geodesic, and that
$\mathcal{M}\smallsetminus\Gamma$ splits diffeomorphically as

$\mathcal{M}\smallsetminus\Gamma\simeq\mathcal{Q}\times S^{2},$

for a smooth 2-manifold $\mathcal{Q}=(\mathcal{M}\smallsetminus\Gamma)/SO(3)$.
The quotient space $\mathcal{M}/SO(3)$, also denoted $\mathcal{Q}$, is viewed
a manifold with boundary, with boundary identified with the centre $\Gamma$.

###  2.2 Double null gauge

A double null gauge consists of functions $u,v\colon\mathcal{Q}\to\mathbb{R}$
such that the level hypersurfaces of $u$ foliate $\mathcal{Q}$ by outgoing
lines which are null with respect to the induced metric on $\mathcal{Q}$, and
the level hypersurfaces of $v$ foliate $\mathcal{Q}$ by ingoing null lines.
The level hypersurfaces of $u$ and $v$ lift to outgoing and incoming null
cones of $\mathcal{M}$ and such a double null gauge can be complemented with
local coordinates $(\theta^{1},\theta^{2})$ on $S^{2}$ to local coordinates
$(u,v,\theta^{1},\theta^{2})$ for $\mathcal{M}$.

For a given double null gauge, the metric $g$ can be written in double null
form

$g=-\Omega^{2}dudv+R^{2}\gamma,$ (2.2)

where $\gamma$ is the unit round metric on $S^{2}$, where $\Omega$ is a
function on $\mathcal{Q}$ and $R\colon\mathcal{Q}\to\mathbb{R}$ is the area
radius function

$R(u,v)=\sqrt{\mathrm{Area}(S_{u,v})/4\pi},$

where $S_{u,v}=\\{(u,v)\\}\times S^{2}\subset\mathcal{M}$. The area radius $R$
extends regularly to $0$ on the centre $\Gamma$. Define

$\lambda=\partial_{v}R,\qquad\nu=\partial_{u}R.$

######  Remark 2.1 (Residual gauge freedom).

There is residual double null gauge freedom present in (2.2). If
$\widetilde{u},\widetilde{v}\colon I\to\mathbb{R}$ are increasing functions,
with $I\subset\mathbb{R}$ a suitable interval, then the change
$u\mapsto\widetilde{u}(u)$, $v\mapsto\widetilde{v}(v)$ preserves the double
null form (2.2). Under such a change, the metric takes the form

$g=-\widetilde{\Omega}^{2}(\widetilde{u},\widetilde{v})d\widetilde{u}d\widetilde{v}+R(\widetilde{u},\widetilde{v})^{2}\gamma=-\widetilde{\Omega}^{2}(\widetilde{u}(u),\widetilde{v}(v))\frac{d\widetilde{u}}{du}(u)\frac{d\widetilde{v}}{dv}(v)dudv+R(\widetilde{u}(u),\widetilde{v}(v))^{2}\gamma.$

It thus follows that

$\displaystyle\Omega^{2}(u,v)$
$\displaystyle=\widetilde{\Omega}^{2}(\widetilde{u}(u),\widetilde{v}(v))\widetilde{u}^{\prime}(u)\widetilde{v}^{\prime}(v),$
(2.3) $\displaystyle R(u,v)$
$\displaystyle=\widetilde{R}(\widetilde{u}(u),\widetilde{v}(v)),$ (2.4)
$\displaystyle\nu(u,v)$
$\displaystyle=\widetilde{\nu}(\widetilde{u}(u),\widetilde{v}(v))\widetilde{u}^{\prime}(u),$
(2.5) $\displaystyle\lambda(u,v)$
$\displaystyle=\widetilde{\lambda}(\widetilde{u}(u),\widetilde{v}(v))\widetilde{v}^{\prime}(v).$
(2.6)

The inverse metric of (2.2) takes the form

$g^{-1}=-2\Omega^{-2}(\partial_{u}\otimes\partial_{v}+\partial_{v}\otimes\partial_{u})+R^{-2}\gamma^{AB}\partial_{\theta^{A}}\otimes\partial_{\theta^{B}},$

where $A$ and $B$ range over $1$ and $2$. The nonvanishing Christoffel symbols
of $g$ are

$\Gamma^{u}_{AB}=\frac{2\lambda
R}{\Omega^{2}}\gamma_{AB},\quad\Gamma^{v}_{AB}=\frac{2\nu
R}{\Omega^{2}}\gamma_{AB},\quad\Gamma^{A}_{Bv}=\frac{\lambda}{R}\delta^{A}_{B},\quad\Gamma^{A}_{Bu}=\frac{\nu}{R}\delta^{A}_{B},\quad\Gamma^{u}_{uu}=\partial_{u}\log\Omega^{2},\quad\Gamma^{v}_{vv}=\partial_{u}\log\Omega^{2},$

and

$\Gamma^{A}_{BC}=\not{\Gamma}^{A}_{BC},$

where $\not{\Gamma}^{A}_{BC}$ are the Christoffel symbols of $(S^{2},\gamma)$.

For a given spherically symmetric double null gauge
$(u,v,\theta^{1},\theta^{2})$ for $(\mathcal{M},g)$, one defines a
corresponding conjugate coordinate system
$(u,v,\theta^{1},\theta^{2},p^{u},p^{v},p^{1},p^{2})$ for $T\mathcal{M}$,
whereby $(u,v,\theta^{1},\theta^{2},p^{u},p^{v},p^{1},p^{2})$ defines the
point

$\big{(}p^{u}\partial_{u}+p^{v}\partial_{v}+p^{A}\partial_{\theta^{A}}\big{)}\big{|}_{(u,v,\theta^{1},\theta^{2})}\in
T\mathcal{M}.$

Such a coordinate system induces a coordinate system
$(u,v,\theta^{1},\theta^{2},p^{v},p^{1},p^{2})$ on the mass shell $P$, with
$p^{u}$ defined by the mass shell relation (1.3), which takes the form

$p^{u}=\frac{R^{2}\gamma_{AB}p^{A}p^{B}}{\Omega^{2}p^{v}}=\frac{L^{2}}{\Omega^{2}R^{2}p^{v}},$
(2.7)

where the angular momentum $L\colon P\to[0,\infty)$ is defined by

$L(x,p)^{2}=R^{4}\gamma_{AB}p^{A}p^{B}.$ (2.8)

The volume form on

$P_{x}=\\{p\in T_{x}\mathcal{M}\mid g_{x}(p,p)=0\\}\subset T_{x}\mathcal{M},$

in this induced coordinate system takes the form

$\frac{R^{2}\sqrt{\det\gamma}}{p^{v}}dp^{v}\wedge dp^{1}\wedge dp^{2}.$

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\colon\mathcal{Q}\times[0,\infty)\times[0,\infty)\to[0,\infty)$ — of $u$,
$v$, $p^{v}$, and $L=(R^{4}\gamma_{AB}p^{A}p^{B})^{\frac{1}{2}}$,

$f(x,p)=f(u,v,p^{v},L).$ (2.9)

The energy momentum tensor $T$ on $\mathcal{M}$ takes the form

$T=T_{uu}(u,v)du^{2}+2T_{uv}(u,v)dudv+T_{vv}(u,v)dv^{2}+T_{AB}(u,v)d\theta^{A}d\theta^{B}.$

In double null gauge, the components take the form

$\displaystyle T_{uu}(u,v)$
$\displaystyle=\frac{\Omega^{4}}{R^{2}}\frac{\pi}{2}\int_{0}^{\infty}\int_{0}^{\infty}f(u,v,p^{v},L)\,p^{v}\,L\,dL\,dp^{v},$
(2.10) $\displaystyle T_{uv}(u,v)$
$\displaystyle=\frac{\Omega^{4}}{R^{2}}\frac{\pi}{2}\int_{0}^{\infty}\int_{0}^{\infty}f(u,v,p^{v},L)\,p^{u}\,L\,dL\,dp^{v},$
(2.11) $\displaystyle T_{vv}(u,v)$
$\displaystyle=\frac{\Omega^{4}}{R^{2}}\frac{\pi}{2}\int_{0}^{\infty}\int_{0}^{\infty}f(u,v,p^{v},L)\,\frac{(p^{u})^{2}}{p^{v}}\,L\,dL\,dp^{v},$
(2.12)

with $p^{u}=p^{u}(u,v,p^{v},L)$ defined by (2.7). It moreover follows from the
mass shell relation (2.7) that

$\gamma^{AB}T_{AB}=\frac{R^{2}}{4\Omega^{2}}T_{uv}.$

###  2.3 The spherically symmetric Einstein–massless Vlasov system in double
null gauge

The Einstein equations (1.4) in spherical symmetry take the form of the
following system of equations for $(\Omega,R,f)$,

$\displaystyle\partial_{u}\partial_{v}(R^{2})$
$\displaystyle=-\frac{\Omega^{2}}{2}+R^{2}T_{uv},$ (2.13)
$\displaystyle\partial_{u}\partial_{v}\log\Omega^{2}$
$\displaystyle=\frac{\Omega^{2}}{2R^{2}}\big{(}1+4\Omega^{-2}\partial_{u}R\partial_{v}R\big{)}-T_{uv}-\frac{\Omega^{2}}{4R^{2}}\gamma^{AB}T_{AB},$
(2.14) $\displaystyle\partial_{u}(\Omega^{-2}\partial_{u}R)$
$\displaystyle=-\frac{1}{2}R\Omega^{-2}T_{uu},$ (2.15)
$\displaystyle\partial_{v}(\Omega^{-2}\partial_{v}R)$
$\displaystyle=-\frac{1}{2}R\Omega^{-2}T_{vv},$ (2.16)

where $\not{g}=R^{2}\gamma$. The mass shell relation (1.3) takes the form

$-\Omega^{2}p^{u}p^{v}+R^{2}\gamma_{AB}p^{A}p^{B}=0.$ (2.17)

In particular, $\Omega^{2}T^{uv}=R^{2}\gamma_{AB}T^{AB}$, or
$4\Omega^{-2}T_{uv}=R^{-2}\gamma^{AB}T_{AB}$, and so equation (2.14) can be
rewritten

$\displaystyle\partial_{u}\partial_{v}\log\Omega^{2}$
$\displaystyle=\frac{\Omega^{2}}{2R^{2}}\big{(}1+4\Omega^{-2}\partial_{u}R\partial_{v}R\big{)}-2T_{uv},$
(2.18)

The Vlasov equation (1.2) in spherical symmetry takes the form

$\displaystyle X(f)\equiv
p^{u}\partial_{u}f+p^{v}\partial_{v}f-\Big{(}R^{2}\gamma_{AB}p^{A}p^{B}\frac{2\nu}{R\Omega^{2}}+(p^{v})^{2}\partial_{v}\log\Omega^{2}\Big{)}\partial_{p^{v}}f-2\frac{\lambda
p^{v}+\nu p^{u}}{R}L\partial_{L}f=0.$ (2.19)

The angular momentum $L\colon P\to[0,\infty)$ 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
$(\mathcal{M},g,f)$ of (2.13)–(2.19) defines a solution
$(\mathcal{Q},\Omega^{2},R,f)$ of the reduced system (2.13)–(2.19) and,
conversely, any solution $(\mathcal{Q},\Omega^{2},R,f)$ of the reduced system
(2.13)–(2.19) defines a solution of (1.2)–(1.5) with
$\mathcal{M}=\mathcal{Q}\times S^{2}$ and $g$ defined by (2.2).

See [27] for more on the spherically symmetric Einstein–massless Vlasov
system.

###  2.4 The functions $t$ and $r$

Define functions $t$ and $r$ by

$t(u,v)=\frac{1}{4}(v+u)^{2},\qquad r(u,v)=v-u.$

When the FLRW spacetime (1.7)–(1.8) is expressed in the double null gauge
introduced in Section 3.1 below, the functions $t$ and $r$ coincide with those
of Section 1.2.

###  2.5 The Cauchy problem

The Einstein–massless Vlasov system (1.2)–(1.5) admits the following local
existence theorem. Initial data consists of a Riemannian $3$ manifold
$(\Sigma,\overline{g})$, along with a symmetric $(0,2)$ tensor $K$ on
$\Sigma$, and a function $f_{1}\colon T\Sigma\to[0,\infty)$, satisfying
constraint equations

$\displaystyle\overline{R}-|K|_{\overline{g}}^{2}+(\overline{\mathrm{tr}}K)^{2}=2\sqrt{\det\overline{g}}\int
f_{1}|p|_{\overline{g}}dp^{1}dp^{2}dp^{3},\qquad\overline{\mathrm{div}}K^{i}-\nabla^{i}\overline{\mathrm{tr}}K=\sqrt{\det\overline{g}}\int
f_{1}p^{i}dp^{1}dp^{2}dp^{3},$

where $\overline{R}$ is the scalar curvature and $\overline{\nabla}$ is the
Levi-Civita connection of $\overline{g}$. A development of initial data
$(\Sigma,\overline{g},K,f_{1})$ is a globally hyperbolic spacetime
$(\mathcal{M},g)$ which admits $\Sigma$ as a Cauchy hypersurface, with induced
first and second fundamental form $\overline{g}$ and $K$ respectively, such
that the restriction of $f$ to the mass shell over $\Sigma$, $P|_{\Sigma}$,
coincides with $f_{1}$, when $T\Sigma$ is appropriately identified with
$P|_{\Sigma}$. See [31] 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 system [7, 8]).

For any smooth initial data set $(\Sigma,\overline{g},K,f_{1})$ for the
Einstein–massless Vlasov system, as above, there exists a unique maximal
development $(\mathcal{M},g,f)$ solving the system (1.2)–(1.5). Moreover, if
$(\Sigma,\overline{g},K,f_{1})$ is spherically symmetric then the maximal
development $(\mathcal{M},g,f)$ is spherically symmetric.

See also [27] 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 $(\mathcal{M},g,f)$ is the maximal development of
a smooth spherically symmetric initial data set
$(\Sigma,\overline{g},K,f_{1})$ for the Einstein–massless Vlasov system
(1.2)–(1.5), and that $\mathcal{Q}$ is the quotient manifold of $\mathcal{M}$
by the action of the $SO(3)$ isometry, as in Section 2.1 above.

######  Theorem 2.3 (Extension principle, around _non-central points_ , for
the spherically symmetric Einstein–massless Vlasov system [11]).

Let $\mathcal{Q}$ be as above. Let $(u_{*},v_{*})$ be such that there exists
$U<u_{*}$ and $V<v_{*}$ such that the characteristic diamond
$\mathcal{D}_{U,V}^{u_{*},v_{*}}=\\{U\leq u<u_{*},V\leq v<v_{*}\\}$ is
contained in $\mathcal{Q}$,

$\mathcal{D}_{U,V}^{u_{*},v_{*}}\subset\mathcal{Q}.$ (2.21)

Assume also that there exists $0<R_{0}<R_{1}$ such that

$\int^{u_{*}}_{U}\int^{v_{*}}_{V}\Omega^{2}(u,v)dudv<\infty,\qquad\text{and}\qquad
R_{0}\leq R(u,v)\leq R_{1},\quad\text{for all
}(u,v)\in\mathcal{D}_{U,V}^{u_{*},v_{*}},$

and that $f(u,v,\cdot,\cdot)$ is compactly supported for all
$(u,v)\in\mathcal{D}_{U,V}^{u_{*},v_{*}}$. Then $(u_{*},v_{*})\in\mathcal{Q}$.

For a proof of Theorem 2.3 see Section 4.3 of [11] (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 [10] in view of the presence of the anti-trapped surfaces in the
spacetimes under consideration (see Remark 3.1 below).

##  3 The FLRW spacetimes

Recall the FLRW spacetimes introduced in Section 1.2. In Section 3.1 a double
null gauge is introduced in each of these spacetimes. Section 3.2 concerns
properties of the massless Vlasov equation on these spacetimes. Though the
main theorem of Section 3.2, 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.

###  3.1 The FLRW metrics in double null gauge

Recall that the FLRW spacetime $(\mathcal{M}_{\circ},g_{\circ})$ takes the
form

$\mathcal{M}_{\circ}=(0,\infty)\times\mathbb{R}^{3},\qquad
g_{\circ}=-dt^{2}+t\left(\frac{4\varrho}{3}\right)^{\frac{1}{2}}\big{(}(dx^{1})^{2}+(dx^{2})^{2}+(dx^{3})^{2}\big{)}.$

For $t\in(0,\infty)$, define

$\Sigma_{t}=\\{t\\}\times\mathbb{R}^{3}\subset\mathcal{M}_{\circ}.$

Define double null coordinates

$u=t^{\frac{1}{2}}-\frac{r}{2},\qquad
v=t^{\frac{1}{2}}+\frac{r}{2},\qquad\text{where
}r=\left(\frac{4\varrho}{3}\right)^{\frac{1}{4}}\sqrt{(x^{1})^{2}+(x^{2})^{2}+(x^{3})^{2}},$

and note that, since $t\geq 0$ and $r\geq 0$ in $\mathcal{M}_{\circ}$,

$v\geq 0,\qquad v\geq u,\qquad v\geq-u,\qquad\text{in }\mathcal{M}_{\circ}.$

The FLRW metric $g_{\circ}$ in the above double null gauge takes the form

$g_{\circ}=-4tdudv+tr^{2}\gamma,\qquad t=\frac{1}{4}(v+u)^{2},\quad r=v-u,$
(3.1)

defined on the quotient manifold

$\mathcal{Q}_{\circ}=\\{(u,v)\in\mathbb{R}^{2}\mid v\geq 0,v\geq
u,v\geq-u\\}.$

The metric $g_{\circ}$ can be written

$g_{\circ}=-\Omega_{\circ}^{2}dudv+R_{\circ}^{2}\gamma,\qquad\text{where}\qquad\Omega^{2}_{\circ}=4t,\qquad
R_{\circ}=t^{\frac{1}{2}}r,$

and moreover

$\lambda_{\circ}=\partial_{v}R_{\circ}=v=\frac{r}{2}+t^{\frac{1}{2}},\qquad\nu_{\circ}=\partial_{u}R_{\circ}=-u=\frac{r}{2}-t^{\frac{1}{2}},\qquad\sqrt{-\det
g_{\circ}}=2t^{2}r^{2}\sqrt{\det\gamma}.$

######  Remark 3.1 (Anti-trapped spheres).

Note that the FLRW spacetime $(\mathcal{M}_{\circ},g_{\circ})$ contains _anti-
trapped_ spheres of symmetry: if $(u,v)\in\mathcal{Q}_{\circ}$ is such that
$u<0$, then $\lambda_{\circ}(u,v)>0$ and $\nu_{\circ}(u,v)>0$.

The non-vanishing Christoffel symbols of (3.1) take the form

$\mathring{\Gamma}^{u}_{AB}=\frac{rv}{2t^{\frac{1}{2}}}\gamma_{AB}=\frac{1}{2}\left(\frac{1}{2t^{\frac{1}{2}}}+\frac{1}{r}\right)r^{2}\gamma_{AB},\qquad\mathring{\Gamma}^{v}_{AB}=-\frac{ru}{2t^{\frac{1}{2}}}\gamma_{AB}=\frac{1}{2}\left(\frac{1}{2t^{\frac{1}{2}}}-\frac{1}{r}\right)r^{2}\gamma_{AB},$
$\mathring{\Gamma}_{uu}^{u}=\mathring{\Gamma}_{vv}^{v}=t^{-\frac{1}{2}},\qquad\mathring{\Gamma}^{A}_{Bv}=\left(\frac{1}{2t^{\frac{1}{2}}}+\frac{1}{r}\right)\delta_{A}^{B},\qquad\mathring{\Gamma}^{A}_{Bu}=\left(\frac{1}{2t^{\frac{1}{2}}}-\frac{1}{r}\right)\delta_{A}^{B},\qquad\mathring{\Gamma}^{A}_{BC}=\not{\Gamma}^{A}_{BC},$

where $\not{\Gamma}^{A}_{BC}$ are the Christoffel symbols of $(S^{2},\gamma)$.

It follows from (1.9) that the null components of the energy momentum tensor
of $f_{\circ}$ satisfy

$T^{\circ}_{uu}=\frac{1}{t},\qquad T^{\circ}_{uv}=\frac{1}{2t},\qquad
T^{\circ}_{vv}=\frac{1}{t}.$ (3.2)

The mass shell relation (2.17) takes the form

$-4p^{u}p^{v}+r^{2}\gamma_{AB}p^{A}p^{B}=0.$ (3.3)

Recall the expression (2.7) for $p^{u}$, and define

$\mathring{p}^{u}=\frac{r^{2}\gamma_{AB}p^{A}p^{B}}{4p^{v}}.$

Where there is no ambiguity, for example in Section 3.2, $\mathring{p}^{u}$
will also be denoted $p^{u}$.

Given a solution $(\mathcal{Q},\Omega^{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)–(1.8) using
values of the coordinates $(u,v)$ coordinates for $\mathcal{Q}$, and the
values of the coordinates $(u,v,p^{v},L)$ for
$\mathcal{Q}\times[0,\infty)\times[0,\infty)$, along with the above double
null gauge for the spatially homogeneous FLRW solutions. For example

$(\Omega^{2}-\Omega^{2}_{\circ})(u,v)=\Omega^{2}(u,v)-(v+u)^{2}.$

Note in particular that, if $f=f_{\circ}$, then

$\Omega^{-4}R^{-2}T_{uu}=\Omega^{-4}_{\circ}R^{-2}_{\circ}T^{\circ}_{uu},\qquad\Omega^{-2}R^{-4}T_{uv}=\Omega^{-2}_{\circ}R^{-4}_{\circ}T^{\circ}_{uv},\qquad
R^{-6}T_{vv}=R^{-6}_{\circ}T^{\circ}_{vv}.$

The asymmetry in $u$ and $v$ arises from the choice of parameterising the mass
shell by $p^{v}$, rather than $p^{u}$.

###  3.2 The massless Vlasov equation on an FLRW background

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 $\mathbb{R}^{3}$ spatial topology).

Let $f$ be a solution of the massless Vlasov equation (1.2) on
$((0,\infty)\times\mathbb{R}^{3},g_{\circ})$, where $g_{\circ}$ is the FLRW
metric (3.1), such that $f_{1}=f|_{\\{t=1\\}}$ is compactly supported. The
components of the energy momentum tensor satisfy

$T_{uu}\lesssim\frac{\|f_{1}\|_{L^{\infty}}}{t^{2}},\qquad
T_{uv}\lesssim\frac{\|f_{1}\|_{L^{\infty}}}{t^{3}},\qquad
T_{vv}\lesssim\frac{\|f_{1}\|_{L^{\infty}}}{t^{4}},$ (3.4)

for $t\geq 1$. Moreover, there exist $U_{0}\leq U_{1}$ such that

$\pi(\mathrm{supp}(f))\subset\\{U_{0}\leq u\leq U_{1}\\},$ (3.5)

where $\pi\colon P\to\mathcal{M}$ denotes the natural projection.

The support of such solutions are depicted in Figure 1.

Figure 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.

The proof of Theorem 3.2 is based on properties of null geodesics in FLRW.
Given a null geodesic $\gamma\colon[1,\infty)\to\mathcal{M}$, write the
tangent vector to $\gamma$ as

$\dot{\gamma}(s)=p^{u}(s)\partial_{u}+p^{v}(s)\partial_{v}+p^{A}(s)\partial_{\theta^{A}}.$

The geodesic equations for $p^{u}(s)$ and $p^{v}(s)$ take the form

$\displaystyle\dot{p}^{u}(s)+t^{-\frac{1}{2}}\big{(}p^{u}(s)\big{)}^{2}+\frac{1}{2}\left(\frac{1}{2t^{\frac{1}{2}}}+\frac{1}{r}\right)r^{2}\gamma_{AB}p^{A}p^{B}$
$\displaystyle=0,$ (3.6)
$\displaystyle\dot{p}^{v}(s)+t^{-\frac{1}{2}}\big{(}p^{v}(s)\big{)}^{2}+\frac{1}{2}\left(\frac{1}{2t^{\frac{1}{2}}}-\frac{1}{r}\right)r^{2}\gamma_{AB}p^{A}p^{B}$
$\displaystyle=0.$ (3.7)

One has equations for $t(s)=t(\gamma(s))$ and $r(s)=r(\gamma(s))$ along the
geodesic

$\dot{t}(s)=t(s)^{\frac{1}{2}}\big{(}p^{v}(s)+p^{u}(s)\big{)},\qquad\dot{r}(s)=p^{v}(s)-p^{u}(s).$
(3.8)

Note in particular that

$\frac{d}{ds}\big{(}t(p^{v}-p^{u})\big{)}=\frac{t}{r}r^{2}\gamma_{AB}p^{A}p^{B}\geq
0.$ (3.9)

######  Proposition 3.3 (Null geodesics in FLRW).

Let $\mathcal{B}\subset P|_{\Sigma_{1}}$ be a compact subset of the mass shell
over $\Sigma_{1}$ and let $\gamma\colon[1,\infty)\to\mathcal{M}$ be a future
directed null geodesic such that $(\gamma(1),\dot{\gamma}(1))\in\mathcal{B}$
and let $v_{0}$ be sufficiently large. There exist $0<c<1<C$, $L_{0}\geq 0$,
and $s_{0}\in[1,\infty)$ such that $v(\gamma(s_{0}))=v_{0}$ and the components
of the tangent vector to $\gamma$ satisfy, for all $s\geq s_{0}$,

$\left(\gamma_{AB}p^{A}(s)p^{B}(s)\right)^{\frac{1}{2}}\leq\frac{L_{0}}{t(s)r(s)^{2}},\qquad
0\leq p^{u}(s)\leq\frac{Cp^{v}(s_{0})}{t(s)r(s)^{2}},\qquad
c\frac{p^{v}(s_{0})}{t(s)}\leq p^{v}(s)\leq C\frac{p^{v}(s_{0})}{t(s)}.$
(3.10)

Moreover, there exists retarded times $U_{0}<0<U_{1}$ such that the $u$
component of $\gamma$ satisfies

$U_{0}\leq u(s)\leq U_{1},\qquad\text{for all }s\geq 1.$ (3.11)

###### Proof.

Provided $v_{0}$ is suitably large, the existence of $s_{0}\in[1,\infty)$ such
that $v(\gamma(s_{0}))=v_{0}$ follows from the compactness of $\mathcal{B}$.

Consider now the bounds (3.10). For simplicity, suppose that $r(\gamma(1))\neq
0$ (otherwise replace $1$ in the following with $1+\epsilon$ for some
$\epsilon>0$). By the conservation of angular momentum (2.20), one has the
conservation law

$L(s)^{2}:=t(s)^{2}r(s)^{4}\gamma_{AB}p^{A}(s)p^{B}(s)=t(1)^{2}r(1)^{4}\gamma_{AB}p^{A}(1)p^{B}(1)\qquad\text{for
all }s\geq 1,$ (3.12)

from which the first of (3.10) trivially follows.

For the remaining inequalities of (3.10), suppose first that $\gamma$ is
radial, so that the conserved quantity $L(s)$ vanishes. It follows from the
mass shell relation (3.3) that either $p^{u}(1)=0$ or $p^{v}(1)=0$. The
geodesic equations (3.6)–(3.7) take the form

$\frac{d}{ds}\left(tp^{u}\right)=0,\qquad\frac{d}{ds}\left(tp^{v}\right)=0,$

and so, as long as the geodesic remains away from the centre where the
coordinate system breaks down, one has $t(s)p^{u}(s)=p^{u}(1)$ and
$t(s)p^{v}(s)=p^{v}(1)$. If $p^{u}(1)=0$ then, by (3.8), $\dot{r}(s)>0$ for
all $s\geq 1$ and so the solution remains away from the centre. The bounds
(3.10) then trivially hold. If $p^{v}(1)=0$ then equations (3.8) imply

$t(s)=\Big{(}\frac{3}{2}(s-1)p^{u}(1)+1\Big{)}^{\frac{2}{3}},\qquad
r(s)=r(1)-\int_{1}^{s}p^{u}(1)\Big{(}\frac{3}{2}(s^{\prime}-1)p^{u}(1)+1\Big{)}^{-\frac{2}{3}}ds^{\prime}.$

It follows that the there is a time $s_{*}\geq 1$ at which the geodesic hits
the centre, i.e. a time at which $r(s_{*})=0$. The geodesic then becomes
outgoing, with

$p^{u}(s)=0\text{ for }s>s_{*},\qquad\text{and}\qquad\lim_{s\downarrow
s_{*}}p^{v}(s)=\lim_{s\uparrow s_{*}}p^{u}(s)=\frac{p^{u}(1)}{t(s_{*})}.$

The bounds (3.10) then again trivially hold, as before.

Suppose now that $\gamma$ is not radial, so that $L\neq 0$, where $L$ is
defined by (3.12). Recall (see (3.9)) that

$\frac{d}{ds}\big{(}t(p^{v}-p^{u})\big{)}=\frac{L^{2}}{tr^{3}}\geq 0.$

If $p^{v}(1)\geq p^{u}(1)$ then it follows that $p^{v}(s)\geq p^{u}(s)$ for
all $s\geq 1$ and, by the latter of (3.8), $r(s)$ is non-decreasing in $s$ and
the geodesic remains away from the centre. If $p^{v}(1)<p^{u}(1)$ then, since
$\int_{0}^{\epsilon}r^{-3}dr=\infty$ for all $\epsilon>0$, it follows that
there exists a time $s_{*}\geq 1$ such that $p^{v}(s_{*})=p^{u}(s_{*})$, and
moreover $r(s)>0$ for all $1\leq s\leq s_{*}$. As above $r(s)$ is then non-
decreasing for $s\geq s_{*}$ and so the geodesic again remains away from the
centre. Using now the mass shell relation (3.3) and the equations (3.6)–(3.7),

$\frac{d}{ds}\left(\frac{r^{2}p^{u}}{p^{v}}\right)=-4r\frac{(p^{u})^{2}}{p^{v}}\leq
0,$

and so (defining $s_{*}=1$ if $p^{v}(1)\geq p^{u}(1)$ and integrating from
$s=s_{*}$),

$r(s)^{2}p^{u}(s)\leq r(s_{*})^{2}p^{v}(s),\qquad\text{for all }s\geq s_{*}.$
(3.13)

Consider the latter of (3.10). It follows from (3.7) and (3.8) that

$\displaystyle\frac{d}{ds}\left(tp^{v}\right)=\frac{t}{2r}r^{2}\gamma_{AB}p^{A}p^{B}=\frac{2}{tr^{3}}\frac{t^{2}r^{4}\gamma_{AB}p^{A}p^{B}}{4}=\frac{2}{t(s)r(s)^{3}}t(s_{0})^{2}r(s_{0})^{2}p^{u}(s_{0})p^{v}(s_{0}),$
(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 $v_{0}$ is sufficiently large (and
hence $s_{0}\geq s_{*}$ is sufficiently large), then $p^{v}(s_{0})\geq
2p^{u}(s_{0})$. Equation (3.9) implies that

$t(s)\big{(}p^{v}(s)-p^{u}(s)\big{)}\geq
t(s_{0})\big{(}p^{v}(s_{0})-p^{u}(s_{0})\big{)}\qquad\text{for all }s\geq
s_{0},$ (3.15)

and so equation (3.14) implies that

$t(s)p^{v}(s)-t(s_{0})p^{v}(s_{0})\leq\int_{s_{0}}^{s}\frac{2}{t(s^{\prime})r(s^{\prime})^{3}}ds^{\prime}t(s_{0})^{2}r(s_{0})^{2}p^{u}(s_{0})p^{v}(s_{0})\\\
\leq\int_{r(s_{0})}^{r(s)}\frac{2}{r^{3}}dr\frac{t(s_{0})r(s_{0})^{2}p^{u}(s_{0})p^{v}(s_{0})}{\big{(}p^{v}(s_{0})-p^{u}(s_{0})\big{)}}\leq
4t(s_{0})p^{v}(s_{0}),$

and the upper bound follows. The estimate for $p^{u}(s)$ in (3.10) finally
follows from returning to (3.13).

Consider now (3.11). The existence of $U_{0}$ follows trivially from the fact
that

$\dot{u}(s)=p^{u}(s)\geq 0.$

The existence of $U_{1}$ follows from the fact that,

$u(s)-u(s_{0})=\int_{s_{0}}^{s}p^{u}(s^{\prime})ds^{\prime}\leq\int_{r(s_{0})}^{r(s)}\frac{C}{tr^{2}}\frac{p^{v}(s_{0})}{p^{v}-p^{u}}dr\leq\frac{2C}{t(s_{0})r(s_{0})},$

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 $\mathcal{B}=\mathrm{supp}(f_{1})$. The
inclusion (3.5) follows immediately from (3.11). Let $v_{0}$ be sufficiently
large so that $R_{\circ}(U_{1},v_{0})>0$. The inclusion (3.5) in particular
implies that,

$cr\leq t^{\frac{1}{2}}\leq Cr,\qquad\text{in }\mathrm{supp}(f)\cap\\{v\geq
v_{0}\\},$

and so (3.10) imply that

$\left(\gamma_{AB}p^{A}p^{B}\right)^{\frac{1}{2}}\leq\frac{L_{0}}{t^{2}},\qquad
0\leq p^{u}\leq\frac{Cp^{v}}{t},\qquad 0\leq p^{v}\leq\frac{C}{t},\qquad
ct^{2}\leq R_{\circ}^{2}\leq Ct^{2}\qquad\text{in
}\mathrm{supp}(f)\cap\\{v\geq v_{0}\\}.$

Hence, in view of the expression (2.10),

$T_{uu}=\frac{\pi}{2}\frac{\Omega_{\circ}^{4}}{R_{\circ}^{2}}\int_{0}^{L_{0}}\int_{0}^{Ct^{-1}}f\,p^{v}\,dp^{v}\,L\,dL\lesssim\frac{\|f_{1}\|_{L^{\infty}}}{t^{2}}.$

Similarly, in view of (2.11) and (2.12),

$T_{uv}=\frac{\pi}{2}\frac{\Omega_{\circ}^{4}}{R_{\circ}^{2}}\int_{0}^{L_{0}}\int_{0}^{Ct^{-1}}f\,p^{u}\,dp^{v}\,L\,dL\lesssim\frac{\|f_{1}\|_{L^{\infty}}}{t^{3}},\qquad
T_{vv}=\frac{\pi}{2}\frac{\Omega_{\circ}^{4}}{R_{\circ}^{2}}\int_{0}^{L_{0}}\int_{0}^{Ct^{-1}}f\,\frac{(p^{u})^{2}}{p^{v}}\,dp^{v}\,L\,dL\lesssim\frac{\|f_{1}\|_{L^{\infty}}}{t^{4}}.$

∎

######  Remark 3.4 (Massless Vlasov equation on Minkowski space).

Consider Minkowski space $(\mathbb{R}^{3+1},m)$, where
$m=-dt^{2}+(dx^{1})^{2}+(dx^{2})^{2}+(dx^{3})^{2}$. In standard double null
coordinates $u=t-r$, $v=t+r$, $r=\sqrt{(x^{1})^{2}+(x^{2})^{2}+(x^{3})^{2}}$,
with respect to which the metric takes the form

$m=-dudv+r^{2}\gamma,$

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

$T_{uu}\lesssim\frac{\|f_{1}\|_{L^{\infty}}}{t^{2}},\qquad
T_{uv}\lesssim\frac{\|f_{1}\|_{L^{\infty}}}{t^{4}},\qquad
T_{vv}\lesssim\frac{\|f_{1}\|_{L^{\infty}}}{t^{6}},$

for $t\geq 1$ [35]. Compare with the rates (3.4) which, with respect to the
double null frame $e_{3}=t^{-\frac{1}{2}}\partial_{u}$,
$e_{4}=t^{-\frac{1}{2}}\partial_{v}$ normalised so that
$g_{\circ}(e_{3},e_{4})=-2$, take the form

$T_{33}\lesssim\frac{\|f_{1}\|_{L^{\infty}}}{t^{3}},\qquad
T_{34}\lesssim\frac{\|f_{1}\|_{L^{\infty}}}{t^{4}},\qquad
T_{44}\lesssim\frac{\|f_{1}\|_{L^{\infty}}}{t^{5}}.$

######  Remark 3.5 (Massless Vlasov equation on FLRW with $\mathbb{T}^{3}$
topology).

The FLRW metric $g_{\circ}$, along with the function $f_{\circ}$, defined by
(1.7)–(1.8), also define a solution of the Einstein–massless Vlasov system on
the manifold $(0,\infty)\times\mathbb{T}^{3}$. 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

$T_{tt}\lesssim\frac{\|f|_{\Sigma_{1}}\|}{t^{2}},\qquad\sum_{i=1}^{3}T_{it}\lesssim\frac{\|f|_{\Sigma_{1}}\|}{t^{\frac{3}{2}}},\qquad\sum_{i,j=1}^{3}T_{ij}\lesssim\frac{\|f|_{\Sigma_{1}}\|}{t},$

for all $t$, for a suitable norm $\|\cdot\|$. Compare again with the rates
(3.4).

##  4 The main theorem

Theorem 1.4 can be more precisely stated as follows. Recall that initial data
for the Einstein–massless Vlasov system consists of
$(\Sigma,\overline{g},K,f_{1})$ (see Section 2.5). Let $\|\cdot\|$ 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
$(\Sigma_{1},g_{\circ}|_{\Sigma_{1}})$,

$\|(\overline{g},K,f_{1})\|=\|\overline{g}\|_{C^{k+1}}+\|K\|_{C^{k}}+\|f_{1}\|_{C^{k}},$

suffices, for example, for $k$ suitably large.

######  Theorem 4.1 (Future stability of FLRW).

Let $\mu:[0,\infty)\to[0,\infty)$ be a smooth function, decaying suitably so
that $\int_{\mathbb{R}^{3}}|p|\mu(|p|^{2})dp$ is finite, such that
$\mu\not\equiv 0$. Let $g_{\circ}$ and $f_{\circ}$ be defined by (1.7), (1.8).
Let $(\Sigma_{1},\overline{g},K,f_{1})$ be a spherically symmetric initial
data set for the Einstein–massless Vlasov system on
$\Sigma_{1}=\\{t=1\\}\subset\mathcal{M}_{\circ}$ such that
$f_{1}-f_{\circ}|_{\Sigma_{1}}$ is compactly supported on $P|_{\Sigma_{1}}$
and $\overline{g}-g_{\circ}|_{\Sigma_{1}}$ and
$K-\frac{1}{2}\partial_{t}g_{\circ}|_{\Sigma_{1}}$ are compactly supported on
$\Sigma_{1}$. There exists $\varepsilon_{0}>0$ such that, if

$\Big{\|}\Big{(}\overline{g}-g_{\circ}|_{\Sigma_{1}},K-\frac{1}{2}\partial_{t}g_{\circ}|_{\Sigma_{1}},f_{1}-f_{\circ}|_{\Sigma_{1}}\Big{)}\Big{\|}\leq\varepsilon_{0},$

then the resulting maximal developement of Theorem 2.2 is future geodesically
complete and there exists $U_{0}<0<U_{1}$ such that, in an appropriately
normalised double null gauge, the associated $(\mathcal{Q},\Omega^{2},R,f)$
satisfies the estimates

$\Big{|}\Omega^{-4}R^{-2}T_{uu}-\Omega^{-4}_{\circ}R^{-2}_{\circ}T^{\circ}_{uu}\Big{|}\lesssim\frac{\varepsilon_{0}}{t^{6}},\quad\Big{|}\Omega^{-2}R^{-4}T_{uv}-\Omega^{-2}_{\circ}R^{-4}_{\circ}T^{\circ}_{uv}\Big{|}\lesssim\frac{\varepsilon_{0}}{t^{8}},\quad\Big{|}R^{-6}T_{vv}-R^{-6}_{\circ}T^{\circ}_{vv}\Big{|}\lesssim\frac{\varepsilon_{0}}{t^{10}},$
(4.1)
$|\Omega^{2}-4t|\lesssim\varepsilon_{0}t,\qquad|R-t^{\frac{1}{2}}r|\lesssim\varepsilon_{0},\qquad|\lambda-v|\lesssim\varepsilon_{0},\qquad|\nu+u|\lesssim\varepsilon_{0},$
(4.2)

in the region $\\{U_{0}<u<U_{1}\\}$, where $t=\frac{1}{4}(v+u)^{2}$ and
$r=v-u$. Moreover,

$\pi(\mathrm{supp}(f-f_{\circ}))\subset\\{U_{0}\leq u\leq U_{1}\\},$

where $\pi\colon P\to\mathcal{M}$ is the natural projection, and the solution
is isometric to the spatially homogeneous FLRW solution (1.7)–(1.8) in the
region $\\{u\geq U_{1}\\}$.

Note that, since the estimates (4.2) imply that $\Omega^{4}R^{2}\sim t^{4}$,
$\Omega^{2}R^{4}\sim t^{5}$ and $R^{6}\sim t^{6}$ in
$\mathrm{supp}(f-f_{\circ})\cap\\{v\geq v_{0}\\}$, 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.

##  5 The proof of the main theorem

This section concerns the proof of Theorem 4.1. In Section 5.1 the proof is
reduced so that the main region of consideration arises from an appropriate
characteristic initial value problem. In Section 5.2 a more convenient
renormalisation of the equations of Section 2.3 are given. Section 5.3
concerns Theorem 5.1, which forms the main content of the proof of Theorem
4.1. In Section 5.4 the proof of Theorem 4.1 is completed.

###  5.1 Cauchy stability and the domain of dependence property

Recall the setting of Theorem 4.1. Let $(\mathcal{M},g)$ denote the unique
maximal development of $(\Sigma_{1},\overline{g},K,f_{1})$, and let
$(\mathcal{Q},\Omega^{2},R,f)$ denote the associated maximal solution of the
reduced system (2.13)–(2.19) (see Section 2.3), in which the residual gauge
freedom, discussed in Remark 2.1, is normalised so that
$\lambda|_{\\{t=1\\}}=\lambda_{\circ}|_{\\{t=1\\}}$ and
$\nu|_{\\{t=1\\}\cup\Gamma}=\nu_{\circ}|_{\\{t=1\\}\cup\Gamma}$.

Consider some $v_{0}>U_{1}$, where $U_{1}>0$ is to be determined later. Given
$\delta_{0}>0$, $L_{0}\geq 0$, define the set

$\mathcal{U}_{\delta_{0},L_{0},v_{0}}=\\{(x,p)\in P\mid v(x)=v_{0},\
2-v_{0}\leq u(x)\leq v_{0}-\delta_{0},\ 0\leq p^{u}\leq 2p^{v},\ L(x,p)\leq
L_{0}\\}\subset P,$

so that, in particular, $r\geq\delta_{0}$ and $t\geq 1$ in
$\mathcal{U}_{\delta_{0},L_{0},v_{0}}$.

By the domain of dependence property, in the region $\\{u\leq U_{0}\\}$, for
some appropriate $U_{0}<0$, the solution coincides with the FLRW solution
(1.7)–(1.8). In particular,

$\Omega^{2}=4t,\qquad\partial_{v}\log\Omega^{2}=t^{-\frac{1}{2}},\qquad
R=t^{\frac{1}{2}}r,\qquad\lambda=v,\qquad\nu=-u,\qquad\text{on
}\\{u=U_{0}\\}.$ (5.1)

Let $t_{0}=t(U_{1},v_{0})=\frac{1}{4}(v_{0}+U_{1})^{2}$. Let $L_{0}>0$ be such
that

$R^{4}\gamma_{AB}p^{A}p^{B}\leq(L_{0})^{2}\qquad\text{in}\qquad\mathrm{supp}(f_{1}-f_{\circ}|_{\Sigma_{1}}),$

and let $\delta_{0}>0$ be such that $r(U_{1},v_{0})=v_{0}-U_{1}>\delta_{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=t_{0}\\}$ and, by
Proposition 3.3, for any future-maximal null geodesic
$\gamma\colon[1,S)\to\mathcal{M}$ such that
$(\gamma(1),\dot{\gamma}(1))\in\mathrm{supp}(f_{1})$, there exists
$s_{0}\in[1,\infty)$ such that

$(\gamma(s_{0}),\dot{\gamma}(s_{0}))\in\mathcal{U}_{\delta_{0},L_{0},v_{0}},$

if $v_{0}$ is sufficiently large and $\varepsilon_{0}>0$ is sufficiently
small. Moreover

$|\Omega^{2}-4t|\lesssim\varepsilon_{0},\quad|\partial_{u}\log\Omega^{2}-t^{-\frac{1}{2}}|\lesssim\varepsilon_{0},\quad|R-t^{\frac{1}{2}}r|\lesssim\varepsilon_{0},\quad|\lambda-v|\lesssim\varepsilon_{0},\quad|\nu+u|\lesssim\varepsilon_{0},\quad\text{on
}\\{v=v_{0}\\}.$ (5.2)

The remainder of the proof will then be concerned with the region $\\{v\geq
v_{0}\\}\cap\\{u\geq U_{0}\\}$.

###  5.2 The renormalised equations

Rather than considering the equations (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,

$\displaystyle\partial_{u}\partial_{v}(R^{2}-R_{\circ}^{2})$
$\displaystyle=-\frac{\Omega^{2}-\Omega^{2}_{\circ}}{2}+(R^{2}-R_{\circ}^{2})T_{uv}+R_{\circ}^{2}(T_{uv}-T^{\circ}_{uv}),$
(5.3)
$\displaystyle\partial_{u}\partial_{v}\log\frac{\Omega^{2}}{\Omega^{2}_{\circ}}$
$\displaystyle=\frac{\Omega^{2}}{2R^{2}}-\frac{\Omega^{2}_{\circ}}{2R^{2}_{\circ}}+\frac{2}{R^{2}}\nu\lambda-\frac{2}{R^{2}_{\circ}}\nu_{\circ}\lambda_{\circ}-2(T_{uv}-T^{\circ}_{uv}),$
(5.4)
$\displaystyle\partial_{u}(\Omega^{-2}\nu-\Omega^{-2}_{\circ}\nu_{\circ})$
$\displaystyle=\frac{1}{2}R_{\circ}\Omega^{-2}_{\circ}T^{\circ}_{uu}-\frac{1}{2}R\Omega^{-2}T_{uu},$
(5.5)
$\displaystyle\partial_{v}(\Omega^{-2}\lambda-\Omega^{-2}_{\circ}\lambda_{\circ})$
$\displaystyle=\frac{1}{2}R_{\circ}\Omega^{-2}_{\circ}T^{\circ}_{vv}-\frac{1}{2}R\Omega^{-2}T_{vv},$
(5.6)

where the solution is compared with the FLRW solution using the double null
gauge of Section 3.1, as discussed at the end of Section 3.1.

For any null geodesic $\gamma\colon[1,S)\to\mathcal{M}$, write the tangent
vector to $\gamma$ as

$\dot{\gamma}(s)=p^{u}(s)\partial_{u}+p^{v}(s)\partial_{v}+p^{A}(s)\partial_{\theta^{A}}.$

The geodesic equations for $p^{u}(s)$ and $p^{v}(s)$ take the form

$\displaystyle\dot{p}^{u}(s)+\partial_{u}\log\Omega^{2}\big{(}p^{u}(s)\big{)}^{2}+\frac{2\lambda
R}{\Omega^{2}}\gamma_{AB}p^{A}p^{B}$ $\displaystyle=0,$ (5.7)
$\displaystyle\dot{p}^{v}(s)+\partial_{v}\log\Omega^{2}\big{(}p^{v}(s)\big{)}^{2}+\frac{2\nu
R}{\Omega^{2}}\gamma_{AB}p^{A}p^{B}$ $\displaystyle=0.$ (5.8)

The following renormalisation of (5.8),

$\displaystyle\frac{d}{ds}(tp^{v})=t(t^{-\frac{1}{2}}-\partial_{v}\log\Omega^{2})\big{(}p^{v}\big{)}^{2}+\Big{(}(R-R_{\circ})-2t^{\frac{1}{2}}(\nu-\nu_{\circ})\Big{)}\frac{t^{\frac{1}{2}}R}{\Omega^{2}}\gamma_{AB}p^{A}p^{B}+\frac{2t^{\frac{3}{2}}R}{\Omega^{2}}\gamma_{AB}p^{A}p^{B},$
(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 $\dot{t}=t^{\frac{1}{2}}(p^{v}+p^{u})$.

###  5.3 The bootstrap theorem

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

Given $T>t_{0}$ and $U_{1}>0$, define the region

$\mathcal{Q}_{U_{1},T}=\\{(u,v)\in\mathcal{Q}\mid U_{0}\leq u\leq U_{1},v\geq
v_{0},\text{and }t_{0}\leq t(u,v)<T\\}.$ (5.10)

######  Theorem 5.1 (Bootstrap theorem).

Suppose that $T>t_{0}$ is such that the maximal development exists up to time
$T$ — in the sense that $(u(t,r),v(t,r))\in\mathcal{Q}$ for all $t_{0}\leq
t<T$, $r\geq 0$ — the difference $f-f_{\circ}$ is supported in the region

$\mathrm{supp}(\pi(f-f_{\circ}))\subset\\{U_{0}\leq u\leq U_{1}\\},$ (5.11)

and, for some fixed $\varepsilon>0$ and $U_{1}>0$, the solution moreover
satisfies, for all $(u,v)\in\mathcal{Q}_{U_{1},T}$,

$\Big{|}\Omega^{-4}R^{-2}T_{uu}-\Omega^{-4}_{\circ}R^{-2}_{\circ}T^{\circ}_{uu}\Big{|}<\frac{\varepsilon}{t^{6}},\quad\Big{|}\Omega^{-2}R^{-4}T_{uv}-\Omega^{-2}_{\circ}R^{-4}_{\circ}T^{\circ}_{uv}\Big{|}<\frac{\varepsilon}{t^{8}},\quad\Big{|}R^{-6}T_{vv}-R^{-6}_{\circ}T^{\circ}_{vv}\Big{|}<\frac{\varepsilon}{t^{10}},$
(5.12) $|\Omega^{2}-4t|<\varepsilon
t,\qquad\Big{|}\partial_{v}\log\Omega^{2}-\frac{1}{t^{\frac{1}{2}}}\Big{|}<\frac{\varepsilon}{t},\qquad\Big{|}\partial_{u}\log\Omega^{2}-\frac{1}{t^{\frac{1}{2}}}\Big{|}<\varepsilon,$
(5.13) $|R-t^{\frac{1}{2}}r|<\varepsilon,\qquad|R^{2}-tr^{2}|<\varepsilon
t,\qquad|\lambda-v|<\varepsilon,\qquad|\nu+u|<\varepsilon.$ (5.14)

Then, if $U_{1}$ is suitably large and $\varepsilon$ is suitably small, there
exists a constant $C>0$ (independent of $\varepsilon$ and $U_{1}$) such that
the inequalities (5.12)–(5.14) hold in $\mathcal{Q}_{U_{1},T}$ with
$\varepsilon$ replaced by $C\varepsilon_{0}$ and (5.11) holds with $U_{1}$
replaced by $U_{1}/2$. Moreover, the solution is isometric to the spatially
homogeneous FLRW solution (1.7)–(1.8) in the region $\\{t_{0}\leq
t(u,v)<T\\}\cap\\{u\geq U_{1}\\}$.

###### Proof.

The proof is divided into several steps. First, the size of the support of
$f-f_{\circ}$ 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 $\\{U_{0}\leq u\leq U_{1}\\}$. Finally, the
solution is shown to be isometric to the spatially homogeneous FLRW solution
(1.7)–(1.8) in the region $\\{u\geq U_{1}\\}$.

Throughout the proof the notation $A\lesssim B$ is used when there exists a
constant $K$, which may depend on $U_{1}-U_{0}$, such that $A\leq KB$.
Constants $0<c\ll 1\ll C$ are always independent of $U_{1}-U_{0}$.

Estimates for the support of $f-f_{\circ}$: Let $\delta_{0}$, $L_{0}$ and
$v_{0}$ be as in Section 5.1, and recall that, for each null geodesic
$\gamma\colon[1,S)\to\mathcal{M}$ emanating from $\mathrm{supp}(f_{1})$, if
$S$ is sufficiently large there exists a time $s_{0}\in[1,S)$ such that
$(\gamma(s_{0}),\dot{\gamma}(s_{0}))\in\mathcal{U}_{\delta_{0},L_{0},v_{0}}$.

For such a null geodesic, write

$\dot{\gamma}(s)=p^{u}(s)\partial_{u}+p^{v}(s)\partial_{v}+p^{A}(s)\partial_{\theta^{A}}.$

By the bootstrap assumption (5.11), the $u$ coordinate of $\gamma$ satisfies

$U_{0}\leq u(s)\leq U_{1},$

for all $s_{0}\leq s<S$. Recall moreover that $v_{0}>U_{1}$ and so
$R(v_{0},U_{1})>0$. It follows in particular that

$t^{\frac{1}{2}}\sim r\sim R^{\frac{1}{2}}\sim v,\qquad\text{along }\gamma.$
(5.15)

There exist constants $0<c\ll 1\ll C$ such that the components of the tangent
vector $\dot{\gamma}$ of any such null geodesic
$\gamma\colon[1,S)\to\mathcal{M}$ emanating from $\mathrm{supp}(f_{1})$
satisfy

$\left(\gamma_{AB}p^{A}(s)p^{B}(s)\right)^{\frac{1}{2}}\leq\frac{2L_{0}}{t(s)r(s)^{2}},\qquad
0\leq p^{u}(s)\leq\frac{Cp^{v}(s_{0})}{t(s)r(s)^{2}},\qquad
c\frac{p^{v}(s_{0})}{t(s)}\leq p^{v}(s)\leq C\frac{p^{v}(s_{0})}{t(s)},$
(5.16)

for all $s\geq s_{0}$. Indeed, let $\gamma\colon[s_{0},S)\to\mathcal{M}$ be a
null geodesic such that
$(\gamma(s_{0}),\dot{\gamma}(s_{0}))\in\mathcal{U}_{\delta_{0},L_{0},v_{0}}$.
The first of (5.16) follows from the conservation of angular momentum (2.20),
which implies that

$R(s)^{4}\gamma_{AB}p^{A}(s)p^{B}(s)=R(s_{0})^{4}\gamma_{AB}p^{A}(s_{0})p^{B}(s_{0})\leq(L_{0})^{2},$
(5.17)

along with the fact that $t^{2}r^{4}\leq 2R^{4}$, which follows from the
bootstrap assumption (5.14).

It follows from the geodesic equation (5.9) that

$\frac{d}{ds}\Big{(}tp^{v}\exp\Big{(}\int_{s_{0}}^{s}h(s^{\prime})p^{v}(s^{\prime})ds^{\prime}\Big{)}\Big{)}=\exp\Big{(}\int_{s_{0}}^{s}h(s^{\prime})p^{v}(s^{\prime})ds^{\prime}\Big{)}\frac{2t^{\frac{3}{2}}R}{\Omega^{2}}\gamma_{AB}p^{A}p^{B},$
(5.18)

where

$h(s)=(t^{-\frac{1}{2}}-\partial_{v}\log\Omega^{2})+\big{(}(R-R_{\circ})-2t^{\frac{1}{2}}(\nu-\nu_{\circ})\big{)}\frac{t^{-\frac{1}{2}}R}{\Omega^{2}}\frac{\gamma_{AB}p^{A}p^{B}}{(p^{v})^{2}}.$

In order to obtain the lower bound of the latter of (5.16), consider the set

$\mathcal{E}=\left\\{s\in[s_{0},S)\mid(tp^{v})(s^{\prime})\geq\frac{1}{2}(tp^{v})(s_{0})\text{
for all }s_{0}\leq s^{\prime}\leq s\right\\}\subset[s_{0},S).$

The set $\mathcal{E}$ is manifestly a closed, connected, non-empty subset of
$[s_{0},S)$. Moreover, if $s\in\mathcal{E}$, then

$\int_{s_{0}}^{s}|h(s^{\prime})|p^{v}(s^{\prime})ds^{\prime}=\int_{v_{0}}^{v(s)}|h|dv\leq
C\int_{v_{0}}^{v(s)}\frac{\varepsilon}{v^{2}}dv\leq C\varepsilon,$ (5.19)

by the fact that

$R(s)^{4}\gamma_{AB}p^{A}(s)p^{B}(s)=R(s_{0})^{4}\gamma_{AB}p^{A}(s_{0})p^{B}(s_{0})=R(s_{0})^{2}\Omega(s_{0})^{2}p^{u}(s_{0})p^{v}(s_{0}),$
(5.20)

along with the bootstrap assumptions (5.13), (5.14) and the equivalence
(5.15). Since

$\frac{d}{ds}\Big{(}tp^{v}\exp\Big{(}\int_{s_{0}}^{s}h(s^{\prime})p^{v}(s^{\prime})ds^{\prime}\Big{)}\Big{)}\geq
0,$

it follows that

$(tp^{v})(s_{0})\leq e^{C\varepsilon}(tp^{v})(s),$

and thus, if $\varepsilon$ is sufficiently small, the set
$\mathcal{E}\subset[s_{0},S)$ is open, and hence equal to $[s_{0},S)$.

For the upper bound of the latter of (5.16), recall equation (5.18) and note
that

$\int_{s_{0}}^{s}\frac{2t^{\frac{3}{2}}R}{\Omega^{2}}\gamma_{AB}p^{A}p^{B}ds^{\prime}=\int_{v_{0}}^{v(s)}\frac{2t^{\frac{5}{2}}}{R^{3}\Omega^{2}}\frac{R^{4}\gamma_{AB}p^{A}p^{B}}{tp^{v}}dv\leq
C\int_{v_{0}}^{v(s)}\frac{p^{v}(s_{0})}{v^{3}}dv\leq Cp^{v}(s_{0}),$

by (5.15), (5.20) and the lower bound for $(tp^{v})(s)$. The upper bound for
$(tp^{v})(s)$ then follows from the estimate (5.19), after integrating
equation (5.18).

The estimate for $p^{u}(s)$ in (5.16) finally follows from the mass shell
relation (2.17), which gives

$\frac{\Omega(s)^{2}}{t(s)}p^{u}(s)=\frac{\big{(}R^{4}\gamma_{AB}p^{A}p^{B}\big{)}(s)}{t(s)p^{v}(s)},$

along with the former and latter of (5.16).

It in particular follows from (5.16) that

$\left(\gamma_{AB}p^{A}p^{B}\right)^{\frac{1}{2}}\leq\frac{CL_{0}}{t^{2}},\qquad
0\leq p^{u}\leq\frac{Cp^{v}}{t},\qquad 0\leq
p^{v}\leq\frac{C}{t},\qquad\text{in }\mathrm{supp}(f-f_{\circ}).$ (5.21)

Consider now the support property (5.11). For any null geodesic
$\gamma\colon[s_{0},S)\to\mathcal{M}$ as above, it follows from (5.16) that,

$u(s)-u(s_{0})=\int_{s_{0}}^{s}p^{u}(s^{\prime})ds^{\prime}\leq\int_{v(s_{0})}^{v(s)}\frac{C}{tv^{2}}\frac{p^{v}(s_{0})}{p^{v}}dv\leq\frac{C}{t(s_{0})v_{0}}.$

Thus, since $t(s_{0}),v_{0}\geq 1$,

$\mathrm{supp}(\pi(f-f_{\circ}))\subset\\{U_{0}\leq u\leq C\\},$ (5.22)

and so (5.11) in fact holds with $U_{1}$ replaced by $U_{1}/2$, provided that
$U_{1}\geq 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_{\circ}$. By (5.22), for any $(u,v,p^{v},L)$,

$|f(u,v,p^{v},L)-f_{\circ}(u,v,p^{v},L)|\leq\varepsilon_{0}\mathbbm{1}_{\\{U_{0}\leq
u\leq U_{1}/2\\}}.$ (5.23)

Define

$L_{\circ}^{2}=R_{\circ}^{4}\gamma_{AB}p^{A}p^{B},$

and note that $R^{-4}L^{2}=R_{\circ}^{-4}L_{\circ}^{2}$. It follows from the
support property (5.21) that

$\displaystyle\Big{|}\Omega^{-4}R^{-2}T_{uu}-\Omega^{-4}_{\circ}R^{-2}_{\circ}T^{\circ}_{uu}\Big{|}\leq\frac{\pi}{2}\int_{0}^{L_{0}}\int_{0}^{Ct^{-1}}|f-f_{\circ}|\,p^{v}dp^{v}R^{-4}L\,dL\lesssim\frac{\varepsilon_{0}}{t^{6}}\mathbbm{1}_{\\{U_{0}\leq
u\leq U_{1}/2\\}}.$

Noting that
$\Omega^{2}R^{-2}p^{u}=\Omega^{2}_{\circ}R^{-2}_{\circ}\mathring{p}^{u}$, it
similarly follows that

$\displaystyle\Big{|}\Omega^{-2}R^{-4}T_{uv}-\Omega^{-2}_{\circ}R^{-4}_{\circ}T^{\circ}_{uv}\Big{|}\leq\frac{\pi}{2}\int_{0}^{L_{0}}\int_{0}^{Ct^{-1}}|f-f_{\circ}|\,\frac{\Omega^{2}p^{u}}{R^{2}}dp^{v}R^{-4}L\,dL\lesssim\frac{\varepsilon_{0}}{t^{8}}\mathbbm{1}_{\\{U_{0}\leq
u\leq U_{1}/2\\}},$

and

$\displaystyle\Big{|}R^{-6}T_{vv}-R^{-6}_{\circ}T^{\circ}_{vv}\Big{|}$
$\displaystyle\leq\frac{\pi}{2}\int_{0}^{L_{0}}\int_{0}^{Ct^{-1}}|f-f_{\circ}|\,\frac{\Omega^{4}(p^{u})^{2}}{R^{4}p^{v}}dp^{v}R^{-4}L\,dL\lesssim\frac{\varepsilon_{0}}{t^{10}}\mathbbm{1}_{\\{U_{0}\leq
u\leq U_{1}/2\\}}.$

Thus

$\displaystyle\Big{|}T_{uu}-\frac{1}{t}\Big{|}$
$\displaystyle\lesssim\frac{\varepsilon_{0}}{t^{2}}\mathbbm{1}_{\\{U_{0}\leq
u\leq U_{1}/2\\}}+t^{-2}|R-R_{\circ}|+t^{-2}|\Omega^{2}-\Omega_{\circ}^{2}|,$
(5.24) $\displaystyle\Big{|}T_{uv}-\frac{1}{2t}\Big{|}$
$\displaystyle\lesssim\frac{\varepsilon_{0}}{t^{3}}\mathbbm{1}_{\\{U_{0}\leq
u\leq U_{1}/2\\}}+t^{-2}|R-R_{\circ}|+t^{-2}|\Omega^{2}-\Omega_{\circ}^{2}|,$
(5.25) $\displaystyle\Big{|}T_{vv}-\frac{1}{t}\Big{|}$
$\displaystyle\lesssim\frac{\varepsilon_{0}}{t^{4}}\mathbbm{1}_{\\{U_{0}\leq
u\leq U_{1}/2\\}}+t^{-2}|R-R_{\circ}|.$ (5.26)

Estimates for the metric quantities in the region $u\leq U_{1}$: For the
metric quantities, consider the region $U_{0}\leq u\leq U_{1}$ and recall that

$t^{\frac{1}{2}}\sim r\sim R^{\frac{1}{2}}\sim v,\qquad\text{for }U_{0}\leq
u\leq U_{1},\quad v\geq v_{0}.$

Let $u$ and $v$ be such that $U_{0}\leq u\leq U_{1}$, $v\geq v_{0}$ and
$t_{0}\leq t(u,v)<T$.

Note that

$\frac{1}{2}\partial_{v}(R^{2}-R_{\circ}^{2})=R_{\circ}(\lambda-\lambda_{\circ}\big{)}+(R-R_{\circ})\lambda.$

Equation (5.3) and the estimate (5.25) then imply that

$\big{|}\partial_{u}\big{(}R_{\circ}(\lambda-\lambda_{\circ}\big{)}+(R-R_{\circ})\lambda\big{|}\lesssim|\Omega^{2}-\Omega_{\circ}^{2}|+|R-R_{\circ}|+\varepsilon_{0}v^{-2}.$

Thus, integrating and using the boundary condition (5.1),

$|\lambda-\lambda_{\circ}|(u,v)\lesssim
v^{-1}|R-R_{\circ}|(u,v)+v^{-2}\int_{U_{0}}^{u}|\Omega^{2}-\Omega_{\circ}^{2}|(u^{\prime},v)+|R-R_{\circ}|(u^{\prime},v)du^{\prime}+\varepsilon_{0}v^{-4}.$
(5.27)

Similarly, since

$\nu-\nu_{\circ}=\Omega^{2}\big{(}\Omega^{-2}\nu-\Omega^{-2}_{\circ}\nu_{\circ}\big{)}-\nu_{\circ}(1-\Omega^{2}\Omega^{-2}_{\circ}),$

it follows from equation (5.5) and the estimate (5.25) and the boundary
condition (5.1) that

$|\nu-\nu_{\circ}|(u,v)\lesssim
v^{-2}|\Omega^{2}-\Omega_{\circ}^{2}|(u,v)+v^{-2}\int_{U_{0}}^{u}|\Omega^{2}-\Omega_{\circ}^{2}|(u^{\prime},v)+|R-R_{\circ}|(u^{\prime},v)du^{\prime}+\varepsilon_{0}v^{-2}.$
(5.28)

Moreover,

$\partial_{u}(R-R_{\circ})=\nu-\nu_{\circ},$

and so (5.28) and the Grönwall inequality imply that

$|R-R_{\circ}|(u,v)\lesssim
v^{-2}\int_{U_{0}}^{u}|\Omega^{2}-\Omega_{\circ}^{2}|(u^{\prime},v)du^{\prime}+\varepsilon_{0}v^{-2}.$
(5.29)

Returning to (5.27) and (5.28), one then has

$|\lambda-\lambda_{\circ}|(u,v)\lesssim
v^{-2}\int_{U_{0}}^{u}|\Omega^{2}-\Omega_{\circ}^{2}|(u^{\prime},v)du^{\prime}+\varepsilon_{0}v^{-3},$
(5.30)

and

$|\nu-\nu_{\circ}|(u,v)\lesssim
v^{-2}|\Omega^{2}-\Omega_{\circ}^{2}|(u,v)+v^{-2}\int_{U_{0}}^{u}|\Omega^{2}-\Omega_{\circ}^{2}|(u^{\prime},v)du^{\prime}+\varepsilon_{0}v^{-2}.$
(5.31)

Now, by equation (5.4) and the estimate (5.25)

$\Big{|}\partial_{u}\partial_{v}\log\frac{\Omega^{2}}{\Omega^{2}_{\circ}}\Big{|}\lesssim
v^{-4}|\Omega^{2}-\Omega_{\circ}^{2}|+v^{-4}|R-R_{\circ}|+v^{-3}|\nu-\nu_{\circ}|+v^{-4}|\lambda-\lambda_{\circ}|+\varepsilon_{0}v^{-4},$
(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

$\Big{|}\log\frac{\Omega^{2}}{\Omega^{2}_{\circ}}\Big{|}\lesssim\varepsilon_{0}+\int_{v_{0}}^{v}\int_{U_{0}}^{u}(v^{\prime})^{-4}|\Omega^{2}-\Omega_{\circ}^{2}|(v^{\prime},u^{\prime})du^{\prime}dv^{\prime}$
(5.33)

Now since, for $\varepsilon$ suitably small,

$|\Omega^{2}-\Omega_{\circ}^{2}|\lesssim\Omega^{2}\Big{|}\log\frac{\Omega^{2}}{\Omega^{2}_{\circ}}\Big{|}\lesssim
v^{2}\Big{|}\log\frac{\Omega^{2}}{\Omega^{2}_{\circ}}\Big{|}\qquad\text{for
}U_{0}\leq u\leq U_{1},$

it follows from the Grönwall inequality, after integrating (5.33) once more in
$u$, that

$|\Omega^{2}-\Omega_{\circ}^{2}|\lesssim\varepsilon_{0}v^{2}.$

Returning to (5.29)–(5.31), it then follows that

$|R-R_{\circ}|\lesssim\varepsilon_{0},\qquad|\lambda-\lambda_{\circ}|\lesssim\varepsilon_{0},\qquad|\nu-\nu_{\circ}|\lesssim\varepsilon_{0}.$

Similarly, returning to (5.32),

$\Big{|}\partial_{v}\log\Omega^{2}-\frac{1}{t^{\frac{1}{2}}}\Big{|}\lesssim\varepsilon_{0}v^{-2},\qquad\Big{|}\partial_{u}\log\Omega^{2}-\frac{1}{t^{\frac{1}{2}}}\Big{|}\lesssim\varepsilon_{0}.$

The region $u\geq U_{1}$: In order to see that the solution is isometric to
the spatially homogeneous FLRW solution (1.7)–(1.8) in the region $u\geq
U_{1}$, recall first (5.23), from which it follows that

$\big{(}R^{-6}T_{vv}\big{)}(U_{1},v)=\big{(}R_{\circ}^{-6}T^{\circ}_{vv}\big{)}(U_{1},v)=t(U_{1},v)^{-4}r(U_{1},v)^{-6}\qquad\text{for
all }v\geq V_{1},$

where $V_{1}$ is such that $R(U_{1},V_{1})=0$. Recall the residual gauge
freedom of Remark 2.1 and consider the following change of gauge. Define
increasing functions $\widetilde{v}(v)$ and $\widetilde{u}(u)$ as solutions of

$\widetilde{v}(V_{1})=U_{1},\qquad\widetilde{v}^{\prime}(v)=\lambda(U_{1},v)\lambda_{\circ}^{-1}(U_{1},\widetilde{v}(v)),\qquad\text{for
}v\geq V_{1},$

and

$\widetilde{u}(u)=u\quad\text{for }u\leq
U_{1},\qquad\widetilde{u}^{\prime}(u)=\frac{1}{\widetilde{v}^{\prime}(V_{1})}\Omega^{2}(u,V_{1})\Omega_{\circ}^{-2}(\widetilde{u}(u),V_{1})\quad\text{for
}u>U_{1},$

respectively. It follows that $\\{u=U_{1}\\}=\\{\widetilde{u}=U_{1}\\}$, and
moreover the relations (2.3) and (2.6) imply that the solution in the
resulting gauge satisfies

$R(U_{1},U_{1})=0,\qquad\Omega^{2}(U_{1},U_{1})=\Omega_{\circ}^{2}(U_{1},U_{1})\qquad\partial_{v}R(U_{1},v)=\partial_{v}R_{\circ}(U_{1},v),\qquad\text{for
all }v\geq U_{1}.$

The Raychaudhuri equation (2.16), restricted to $u=U_{1}$, takes the form

$\partial_{v}(\Omega^{-2}\partial_{v}R)(U_{1},v)=-\frac{1}{2}R^{7}\Omega^{-2}\big{(}R_{\circ}^{-6}T^{\circ}_{vv}\big{)}(U_{1},v).$

It then follows that

$R(U_{1},v)=R_{\circ}(U_{1},v),\qquad\Omega^{2}(U_{1},v)=\Omega_{\circ}^{2}(U_{1},v)\qquad\text{for
all }v\geq U_{1}.$

By uniqueness for the characteristic initial value problem, it follows that
$R=R_{\circ}$ and $\Omega^{2}=\Omega^{2}_{\circ}$ for all $u\geq U_{1}$. ∎

###  5.4 The conclusion of the proof of the main theorem

Recall the set $\mathcal{Q}_{U_{1},T}\subset{Q}$ defined by (5.10). Consider
the set $\mathcal{A}\subset[t_{0},\infty)$ of times $T\in[t_{0},\infty)$ such
that the solution exists up to time $T$ — in the sense that
$(u(t,r),v(t,r))\in\mathcal{Q}$ for all $t_{0}\leq t<T$, $r\geq 0$ — and
moreover satisfies the estimates (5.12)–(5.14) for all
$(u,v)\in\mathcal{Q}_{U_{1},T}$ for some $\varepsilon>0$ sufficiently small so
that the conclusion of Theorem 5.1 holds. The set
$\mathcal{A}\subset[t_{0},\infty)$ is manifestly connected and open, and is
moreover non-empty by local existence and Cauchy stability, provided
$\varepsilon_{0}$ is sufficiently small. By Theorem 2.3, the bounds
(5.13)–(5.14) imply that $(u(T,r),v(T,r))\in\mathcal{Q}$ for all $r\geq 0$
such that $u(T,r)\leq U_{1}$. By Theorem 5.1, the solution is isometric to the
spatially homogeneous FLRW solution (1.7)–(1.8) in the region $\\{t_{0}\leq
t(u,v)<T\\}\cap\\{u\geq U_{1}\\}$ and thus $(u(T,r),v(T,r))\in\mathcal{Q}$ for
all $r\geq 0$. Theorem 5.1 moreover guarantees that the bounds (5.12)–(5.14)
cannot be saturated, and thus the set $\mathcal{A}$ is moreover closed, and
hence equal to $[t_{0},\infty)$.

The future geodesic completeness follows from the bounds on the Christoffel
symbols arising from (5.13)–(5.14) (as, for example in Section 16 of [25]).

## References

  * [1] Lars Andersson and David Fajman.  Nonlinear stability of the Milne model with matter.  Communications in Mathematical Physics, 378(1):261–298, 2020.
  * [2] Lars Andersson and Vincent Moncrief.  Einstein spaces as attractors for the Einstein flow.  Journal of Differential Geometry, 89(1):1–47, 2011.
  * [3] Hamed Barzegar and David Fajman.  Stable cosmologies with collisionless charged matter.  Journal of Hyperbolic Differential Equations, 19(04):587–634, 2022\.
  * [4] Jacob Bedrossian, Nader Masmoudi, and Clément Mouhot.  Landau damping: paraproducts and Gevrey regularity.  Annals of PDE, 2:1–71, 2016.
  * [5] Jacob Bedrossian, Nader Masmoudi, and Clément Mouhot.  Linearized wave-damping structure of Vlasov–Poisson in $\mathbb{R}^{3}$.  SIAM Journal on Mathematical Analysis, 54(4):4379–4406, 2022.
  * [6] Gerhard Börner.  The early universe: facts and fiction.  Springer Science & Business Media, 2013.
  * [7] Yvonne Choquet-Bruhat.  Problème de Cauchy pour le système intégro-différentiel d’Einstein–Liouville.  In Annales de l’institut Fourier, volume 21, pages 181–201, 1971\.
  * [8] Yvonne Choquet-Bruhat and Robert Geroch.  Global aspects of the Cauchy problem in general relativity.  Communications in Mathematical Physics, 14:329–335, 1969.
  * [9] Mihalis Dafermos.  A note on the collapse of small data self-gravitating massless collisionless matter.  Journal of Hyperbolic Differential Equations, 3(04):589–598, 2006\.
  * [10] Mihalis Dafermos and Alan D Rendall.  An extension principle for the Einstein–Vlasov system in spherical symmetry.  In Annales Henri Poincaré, volume 6, pages 1137–1155. Springer, 2005.
  * [11] Mihalis Dafermos and Alan D Rendall.  Strong cosmic censorship for surface-symmetric cosmological spacetimes with collisionless matter.  Communications on Pure and Applied Mathematics, 69(5):815–908, 2016\.
  * [12] Jurgen Ehlers, P Geren, and Rainer K Sachs.  Isotropic solutions of the Einstein–Liouville equations.  Journal of Mathematical Physics, 9(9):1344–1349, 1968.
  * [13] David Fajman.  The nonvacuum Einstein flow on surfaces of negative curvature and nonlinear stability.  Communications in Mathematical Physics, 353(2):905–961, 2017.
  * [14] David Fajman, Jérémie Joudioux, and Jacques Smulevici.  The stability of the Minkowski space for the Einstein–Vlasov system.  Analysis & PDE, 14(2):425–531, 2021.
  * [15] David Fajman, Maximilian Ofner, and Zoe Wyatt.  Slowly expanding stable dust spacetimes.  arXiv preprint arXiv:2107.00457, 2021.
  * [16] David Fajman and Zoe Wyatt.  Attractors of the Einstein–Klein–Gordon system.  Communications in Partial Differential Equations, 46(1):1–30, 2021\.
  * [17] Grigorios Fournodavlos.  Future dynamics of FLRW for the massless-scalar field system with positive cosmological constant.  Journal of Mathematical Physics, 63(3):032502, 2022.
  * [18] Emmanuel Grenier, Toan T Nguyen, and Igor Rodnianski.  Landau damping for analytic and Gevrey data.  arXiv preprint arXiv:2004.05979, 2020.
  * [19] Mahir Hadžić and Jared Speck.  The global future stability of the FLRW solutions to the dust-Einstein system with a positive cosmological constant.  Journal of Hyperbolic Differential Equations, 12(01):87–188, 2015\.
  * [20] Daniel Han-Kwan, Toan T Nguyen, and Frédéric Rousset.  On the linearized Vlasov–Poisson system on the whole space around stable homogeneous equilibria.  Communications in Mathematical Physics, 387:1405–1440, 2021.
  * [21] Stephen W Hawking and George FR Ellis.  The large scale structure of space-time.  Cambridge university press, 2023.
  * [22] Alexandru Ionescu, Benoit Pausader, Xuecheng Wang, and Klaus Widmayer.  Nonlinear Landau damping for the Vlasov-Poisson system in $\mathbb{R}^{3}$: the Poisson equilibrium.  arXiv preprint arXiv:2205.04540, 2022.
  * [23] Michael K-H Kiessling.  The “Jeans swindle”: a true story–mathematically speaking.  Advances in Applied Mathematics, 31(1):132–149, 2003.
  * [24] Evgenii Mikhailovich Lifshitz.  On the gravitational stability of the expanding universe.  J. Phys., 10:116–129, 1946.
  * [25] Hans Lindblad and Igor Rodnianski.  Global existence for the Einstein vacuum equations in wave coordinates.  Communications in mathematical physics, 256:43–110, 2005.
  * [26] Hans Lindblad and Martin Taylor.  Global stability of Minkowski space for the Einstein–Vlasov system in the harmonic gauge.  Archive for rational mechanics and analysis, 235(1):517–633, 2020\.
  * [27] Georgios Moschidis.  The characteristic initial-boundary value problem for the Einstein–massless Vlasov system in spherical symmetry.  arXiv preprint arXiv:1812.04274, 2018.
  * [28] Georgios Moschidis.  A proof of the instability of AdS for the Einstein-massless Vlasov system.  Inventiones mathematicae, pages 1–206, 2022.
  * [29] Clément Mouhot and Cédric Villani.  On Landau damping.  Acta Mathematica, 207:29–201, 2011.
  * [30] Gerhard Rein and Alan D Rendall.  Global existence of solutions of the spherically symmetric Vlasov–Einstein system with small initial data.  Communications in mathematical physics, 150:561–583, 1992.
  * [31] Hans Ringström.  On the topology and future stability of the universe.  OUP Oxford, 2013.
  * [32] Igor Rodnianski and Jared Speck.  The nonlinear future stability of the FLRW family of solutions to the irrotational Euler–Einstein system with a positive cosmological constant.  Journal of the European Mathematical Society, 15(6):2369–2462, 2013\.
  * [33] Jared Speck.  The nonlinear future stability of the FLRW family of solutions to the Euler–Einstein system with a positive cosmological constant.  Selecta Mathematica, 18(3):633–715, 2012.
  * [34] Jared Speck.  The stabilizing effect of spacetime expansion on relativistic fluids with sharp results for the radiation equation of state.  Archive for Rational Mechanics and Analysis, 210:535–579, 2013\.
  * [35] Martin Taylor.  The global nonlinear stability of Minkowski space for the massless Einstein–Vlasov system.  Annals of PDE, 3(1):9, 2017.


# Graph attentive feature aggregation
for text-independent speaker verification

###### Abstract

The objective of this paper is to combine multiple frame-level features into a
single utterance-level representation considering pairwise relationships. For
this purpose, we propose a graph attentive feature aggregation module by
interpreting each frame-level feature as a node of a graph. The inter-
relationships between all possible pairs of features, typically exploited
indirectly, can be directly modeled using a graph. The module comprises a
graph attention layer and a graph pooling layer followed by a readout
operation. The graph attention layer first models the non-Euclidean data
manifold between different nodes. Then, the graph pooling layer discards less
informative nodes considering the significance of the nodes. Finally, the
readout operation combines the remaining nodes into a single representation.
We employ two recent systems with different input features and architectures,
SE-ResNet and RawNet2, and demonstrate that the proposed feature aggregation
module consistently shows a relative improvement over $10\%$, compared to the
baseline.

Index Terms—  speaker recognition, speaker verification, attention, graph
attention networks, deep learning

##  1 Introduction

Speaker verification (SV) can be used for various authentication scenarios
where the identity of a given speech is compared to that of the claimed
speaker. In SV, speaker representations are derived by first extracting the
frame-level features and then aggregating them. The network extracting the
frame-level features is referred to as trunk network (e.g., convolutional
neural network (CNN) and x-vector [jung2020improving, chung2020in,
lin2020wav2spk, snyder2018x, tang2019deep, wu2020vector]). After extracting
the frame-level features, various techniques including gated recurrent network
(GRU), learnable dictionary encoding (LDE) are used on top of the trunk
network to aggregate frame-level features into a single utterance-level
feature [li2017deep, jung2019short, zhang2019seq2seq, jung2019rawnet,
jung2020improved, snyder2017deep, zhu2018self, okabe2018attentive,
kye2021cross, cai2018exploring, xie2019utterance, jung2019spatial]. The
condensed utterance-level representation ultimately represents the entire
character of an utterance, unlike speech recognition where each frame-level
representation is of worth. Hence, feature aggregation plays an essential role
and a number of studies in the recent literature have focused on developing
feature aggregation methods.

A group of studies employed recurrent layers for feature aggregation. These
studies modeled the frame-level representations in sequential order
[li2017deep, jung2019short, zhang2019seq2seq, jung2019rawnet,
jung2020improved]. Another group of research integrated frame-level features
using LDE [cai2018exploring, xie2019utterance, jung2019spatial]. Similar to
the recurrent layers, the LDE layer aggregates features over time whereas it
further utilizes statistics. Sequential modeling, a foundation of the above
two approaches, has been used in SV based on its success in language
processing and speech recognition, where the sequential order could be prior
knowledge. However, recent studies found that sequential information may not
be the key, especially in text-independent SV (TI-SV) [okabe2018attentive,
desplanques2020ecapa].

Other groups of research exploited attention-based approaches-
[snyder2017deep, zhu2018self, okabe2018attentive, kye2021cross], which have
recently become the mainstream in TI-SV. Attention-based methods compute
attention coefficients with projection, followed by a softmax function. This
approach can select informative features regardless of their sequence. The
importance of each feature is decided by a dot product between all frame-level
features and the projection weight vector. Overall features are aggregated
with a relative comparison of attention coefficients although softmax non-
linearity exclusively increases the difference between values.

The field of graph neural networks is recently getting attention by
incorporating the advantages of graph structures and deep neural networks
[goyal2018graph, wu2020comprehensive]. Interpreting high-dimensional
representations as nodes of a graph, graph neural networks can model the non-
Euclidean data manifold within nodes, including inter-node relationships.
Particularly, recent architectures such as graph convolutional network
[kipf2016semi] and graph attention network (GAT) [velivckovic2017graph]
combined with various graph pooling layers [gao2019graph, lee2019self,
yuan2020structpool] have demonstrated competitive performance even in image
and audio domains [jung2021graph, tak2021graph, tak2021end].

Inspired by the the recent success of graph neural networks, we argue that the
feature aggregation could be further improved by explicitly comparing and
modeling the inter-relationship between frame-level features. The feature
aggregation is a task to express characteristics relative to other features
and condense them into a single representation. Hence, the inter-relationship
between features extends the ability of its expressiveness and could be the
vital information itself. To this end, we propose a novel feature aggregation
module to model inter-relationship by interpreting each frame-level feature as
a node of a graph. The proposed module leverages the correlation of all
possible pairs of nodes to obtain an utterance-level representation. In other
words, entire frame-level features are aggregated considering their inter-
relationships in addition to their intrinsic characteristics by utilizing a
graph. Specifically, the aggregation involves a GAT, a graph pooling layer,
and a readout operation [velivckovic2017graph, gao2019graph]. The GAT assigns
different weights to each connection of nodes (edges) and models their
relationships. Then, the graph pooling layer discards less informative nodes.
Finally, a single utterance-level representation is derived using the readout
operation.

Fig. 1:  Proposed graph attentive feature aggregation module. After extracting
frame-level features using a trunk network, a graph is constructed using each
sequence element as a node. Each node can have a dimensionality equal to
either the number of channels (raw waveform, 1D-CNN) or the number of channels
$\times$ the number of frequency bins (mel-spectrogram, 2D-CNN). Graph
attention network first aggregates nodes features, and then the graph pooling
reduces the number of nodes. The readout function integrates nodes into a
single representation (best viewed in color).

The effectiveness of the proposed feature aggregation framework is validated
using two strong baselines, SE-ResNet [heo2020clova] and RawNet2
[jung2020improved], comprising different input features and architectures.
Experimental results on the VoxCeleb datasets [nagrani2017voxceleb,
chung2018voxceleb2] demonstrate more than 10% consistent improvement over the
corresponding baseline. Using the proposed graph aggregation module, we could
compose more lightweight models compared to self-attentive pooling and GRU-
based aggregations.

##  2 graph attentive feature aggregation

In this section, we introduce the graph attentive feature aggregation module.
The proposed module is located after the extraction of frame-level features
using a trunk network (RawNet2 or SE-ResNet). The proposed approach consists
of three components: i) graph attention layer, ii) graph pooling layer, and
iii) readout. We first describe how we formulate a graph from a feature map
and perform GAT to integrate the overall features. Then, we explain how the
graph pooling layer discards less informative nodes and the readout combines
the results into a single representative feature. The overall scheme is
illustrated in Figure 1 using RawNet2 (the scheme is identical for SE-ResNet
except the dimensionality of each node).

###  2.1 Graph attention layer

We first formulate a graph using frame-level features extracted from a trunk
network.Specifically, we interpret each sequence element as a node of a graph;
for RawNet2, each node has a dimensionality equal to the number of filters and
for SE-ResNet, the dimensionality equals to the number of filters $\times$ the
number of frequency bins. Note that we design a fully-connected (complete)
graph where all possible edges exist between pairs of nodes. Nodes can be
directly compared considering all combinations/relationships because all nodes
have bi-directional edges. Hence, entire frame-level features are aggregated
considering their inter-relationships as well as intrinsic characteristics of
their own by utilizing the graph attention layer. The GAT process can be
described as follows:

$\displaystyle\mathbf{n^{\prime}}=\mathbf{x}W,$ (1) $\displaystyle
e_{ij}=LeakyReLU(\sigma(\mathbf{n^{\prime}}_{i},\mathbf{n^{\prime}}_{j})),$
(2)
$\displaystyle\mathbf{a}_{ij}=\mathsf{softmax}(e_{ij})=\frac{\mathsf{exp}(e_{ij})}{\Sigma_{t\in
N}\mathsf{exp}(e_{it})},$ (3) $\displaystyle\mathbf{n}_{i}=\sum_{j\in
N}\mathbf{a}_{ij}\mathbf{n^{\prime}}_{j}.$ (4)

Let $\mathbb{G}$ be a complete graph comprising a set of nodes $\mathbf{x}$,
$\mathbf{x}\in\mathbb{R}^{N\times F}$, where $N$ is the number of nodes and
$F$ refers to the dimensionality of each node. In the GAT layer, $\mathbf{x}$
is first projected into F’ dimensional space using a matrix multiplication
with $W$, $W\in\mathbb{R}^{F\times F^{\prime}}$ resulting in
$\mathbf{n^{\prime}}$, $\mathbf{n^{\prime}}\in\mathbb{R}^{N\times F^{\prime}}$
(Equation (1)). Then, we calculate attention scores where an attention score
$\mathbf{a}_{ij}$ between $i$-th node $\mathbf{n^{\prime}}_{i}$ and $j$-th
node $\mathbf{n^{\prime}}_{j}$ ($\mathbf{n^{\prime}}_{i}$,
$\mathbf{n^{\prime}}_{j}$ $\in\mathbb{R}^{F^{\prime}}$) is calculated through
Equation (2) and (3). We derive $e_{ij}$, by first concatenating two nodes
$n^{\prime}_{i}$ and $n^{\prime}_{j}$ and then projecting to a scalar value
through a dot product with $\sigma$, a learnable parameter,
$\sigma\in\mathbb{R}^{2F^{\prime}\times 1}$ followed by a Leaky ReLU non-
linearity [maas2013rectifier]. The softmax function is applied to all edge
scores associated with $\mathbf{n^{\prime}}_{i}$ to exclusively emphasize the
scores. Using the calculated attention scores, the GAT performs self-attention
on node $\mathbf{n^{\prime}}_{i}$ as in Equation (4). The GAT output,
$\mathbf{n}$, is then fed to a graph pooling layer.

Table 1:  SE-ResNet architecture applying the proposed graph aggregation method. ReLU and batch normalization layers are omitted. Squeeze and excitation [hu2018squeeze] is performed after all convolutional blocks. The numbers denoted in brackets refer to [size$\times$ size, filters] and those in the output refer to (filters, frequency, time). Output of the last fully-connected layer is used as the speaker embedding.  layer name | Filters | Output
---|---|---
conv1 |  3$\times$3, 32, stride 1 |  (32$\times$40$\times$T)
conv2 |  $\left[\begin{tabular}[]{c}3\times 3, 32\\\ 3\times 3, 32\\\ \end{tabular}\right]\times 3$, stride 1 |  (32$\times$40$\times$T)
conv3 |  $\left[\begin{tabular}[]{c}3\times 3, 64\\\ 3\times 3, 64\\\ \end{tabular}\right]\times 4$, stride 2 |  (64$\times$20$\times$T/2)
conv4 |  $\left[\begin{tabular}[]{c}3\times 3, 128\\\ 3\times 3, 128\\\ \end{tabular}\right]\times 6$, stride 2 |  (128$\times$10$\times$T/4)
conv5 |  $\left[\begin{tabular}[]{c}3\times 3, 128\\\ 3\times 3, 128\\\ \end{tabular}\right]\times 3$, stride 2 |  (128 $\times$ 5 $\times$T/8)
aggregation | GAT | (640)
fc | FCN, 256 | (256)
Table 2: Application of proposed graph module to SE-ResNet and RawNet2. SAP in SE-ResNet and GRU in RawNet2 are both replaced by GAT.  Feature extractor | Aggregation | # Params | EER (%)
---|---|---|---
SE-ResNet | SAP | 6.0M | 1.98
SE-ResNet | GAT | 5.4M | 1.86
RawNet2 | GRU | 13.2M | 2.48
RawNet2 | GAT | 9.9M | 2.23
Table 3:  Comparison of applying different pooling ratios in graph pooling layer. A higher pooling ratio means more nodes remain. Pooling method | EER (%)
---|---
SE-ResNet | RawNet2
w/o gPool | 1.86 | 2.23
Top 11% gPool | 2.00 | 2.39
Top 33% gPool | 1.91 | 2.17
Top 80% gPool | 1.75 | 2.15
Table 4: Comparison of various readout methods. “Combine by concat” indicates that the results of summation, standard deviation (std), minimum, and maximum are concatenated with reduced dimensions in the quarter. Operation | EER (%)
---|---
Mean | 2.38
Sum | 2.23
Max | 2.3
Combine by concat | 2.38

###  2.2 Graph pooling layer and readout

We utilize a graph pooling layer followed by the readout operation to obtain
an utterance-level feature. The graph pooling layer reduces the original graph
into a sub-graph by removing less informative nodes as depicted in Figure 1.
For the graph pooling, we exploit gPool introduced in the Graph U-Net
architecture [gao2019graph]. There are various other graph pooling methods,
but many of them assume the hierarchical architecture of the graph. As it is
difficult to define a hierarchical relationship between frame-level features
in audio data, we adopt gPool, which does not require the graph to have a
hierarchical structure. gPool consists of the following processes: projection,
top-$k$ node selection, and a gating mechanism. Formally, gPool can be denoted
as follows:

$\displaystyle\mathrm{y}=\mathbf{n}^{l}\mathrm{p}^{l}/\lVert{\mathrm{p}}^{l}\lVert,$
(5) $\displaystyle\mathrm{index}=\mathrm{rank}(\mathrm{y},K),$
$\displaystyle\mathrm{\tilde{y}}=\mathrm{sigmoid}(\mathrm{y}(\mathrm{index})),$
$\displaystyle\tilde{\mathbf{n}}^{l}=\mathbf{n}^{l}(\mathrm{index},:)$
$\displaystyle\mathbf{n}^{l+1}=\tilde{\mathbf{n}}^{l}\odot\mathrm{\tilde{y}}$

A set of nodes in the graph pooling layer $l$ is represented as
${\mathbf{n}}^{l}$. Note that the adjacency matrix is omitted because we
utilize a complete graph which is different with the original paper of
[gao2019graph]. In the projection stage, gPool obtains score vector
$\mathrm{y}$ by estimating the scalar projection values of each node, using a
trainable projection vector $\mathrm{p}\in\mathbb{R}^{F^{\prime}}$. Based on
the calculated score vector $\mathrm{y}$, top-$k$ nodes are selected to get
the indices of the highest $k$ nodes. In the gate stage, element-wise product
$\odot$, is conducted between the pooled feature map $\mathbf{\tilde{n}}^{l}$
and the selected node scores $\mathrm{\tilde{y}}$.

After the graph pooling, we conduct the readout operation, which combines the
processed nodes into a single node. We explore various types of readout
operation where the best result is obtained using the summation of nodes:

$\displaystyle U=\sum_{k=1}^{K}\mathbf{n}_{k}^{l+1}.$ (6)

Here, $U$ is the aggregated utterance-level feature and $K$ nodes $n_{k}$ are
integrated to $\mathbf{n}$. The results of various readout mechanisms can be
found in Section 4.

Table 5: Results of comparison with recent state-of-the-art systems. The two baselines of this study are each compared with various architectures that use the same input feature. | Input Feature | Front-end | Aggregation | Loss | EER (%)
---|---|---|---|---|---
Jung et al. [jung2020improving] | Mel-spectrogram | Fast RexNet-34 | LDE | NP + Softmax | 1.98
Kye et al. [kye2021cross] | Fast RexNet-34 | CAP | NP + Softmax | 1.88
Chung et al. [chung2020in] | Thin ResNet-34 | SAP | AP | 2.21
Liu et al. [liu2021effective] | Dense-Residual | ABP | Softmax | 2.54
India et al. [india2021double] | adapted VGG | Double MHA | AM-Softmax | 3.19
Ours | SE-ResNet | SAP | AP + AAM-Softmax | 1.98
Ours(Proposed) | SE-ResNet | GAT | AP + AAM-Softmax | 1.75
Lin et al. [lin2020wav2spk] | Raw waveform | wav2spk | Gating + SP | AM-Softmax | 3.00
Zhu et al. [zhu2021vector] | Y-vector | SP | AM-Softmax | 2.60
Jung et al. [jung2020improved] | RawNet2 | GRU | Softmax | 2.48
Ours(Proposed) | RawNet2 | GAT | AM-Softmax | 2.15
Table 6: Results of ablation studies. The first row shows the result of applying GAT and GRU sequentially. The lower two rows show the result of applying two GAT layers with different pooling architectures. System | EER (%)
---|---
GAT - GRU | 4.16
2GAT - Global architecture | 2.66
2GAT - Hierarchical architecture | 2.55

##  3 Experiments

###  3.1 Datasets

All experiments are performed using the VoxCeleb1&2 datasets
[nagrani2017voxceleb, chung2018voxceleb2]. We train the model using the
development subset of VoxCeleb2 that includes the utterances from 5,994
speakers. Then, the evaluation is performed using the original trial that uses
VoxCeleb1’s evaluation subset.

###  3.2 Implementation details

Both baselines using a mel-spectrogram and a raw waveform are implemented
based on the PyTorch framework [paszke2019pytorch]. Graphs are implemented
using Deep Graph Library [wang2019deep]. We adjust the number of attention
heads and use 16 heads in RawNet2 and 32 heads in SE-ResNet based on empirical
results.

SE-ResNet. Our mel-spectrogram baseline is most similar to the architecture of
[heo2020clova], but overall details are adjusted. We use log mel-spectrograms
as inputs extracted with 40-dimensional mel-filterbanks and 1,024 point FFT.

Then, log mel-spectrogram is generated using a hamming window of width 25ms
and step 10ms. The mel-spectrogram baseline is optimized using Adam optimizer
which uses a learning rate of 0.001, and the learning rate is decreased by 5%
in each epoch. We use two types of loss functions, additive angular margin
softmax (AAM-softmax) [deng2019arcface] and angular prototypical loss (AP)
[wang2017deep] to consider inter-class as well as intra-class covariance. We
employ a margin of 0.3 and a scale of 30 for AAM-softmax. For AP, we use a
mini-batch size of 200, where each batch contains two utterances of the
speakers. The system is trained for 100 epochs. Our overall architecture
applying the proposed method is presented in Table 1.

RawNet2. RawNet2 [jung2020improved] is the end-to-end system that are fed by
raw waveforms directly without preprocessing techniques. For mini-batch
construction, utterances are either cropped or duplicated (concatenated) into
59,049 samples ($\approx$ 3.69s) in the training phase, following
[jung2020improved]. In the evaluation phase, no adjustments are made to the
length. We modify several details from RawNet2 [jung2020improved] in the
process of adjusting the GAT as follows: i) exclude sinc-conv layer, as it
slows down training time in spite of showing similar performance, ii) replace
softmax with additive margin softmax (AM-softmax) [wang2018cosface], iii)
reduce the dimensions by half for the last fully-connected layer before
speaker embedding (1024 dimensions to 512 dimensions).

##  4 Results

Application of GAT. In Table 2, we first address the effectiveness of the
graph attention layer on our two baselines: SE-ResNet and RawNet2. Here, the
graph pooling is temporarily excluded and summation readout operation is
applied. Adoption of the graph attention layer reduced the EER relatively by
6% and 10%, respectively. In addition, we confirm that the number of
parameters also decreases (by almost 25% in the case of RawNet2 (13.2M to
9.9M)).

Graph pooling. In Table 3, we explore the effect of applying the graph pooling
layer between the graph attention layer and the readout operation with diverse
pooling ratios. The denoted percentage of pooling refers to the remaining
ratios of nodes; higher ratios mean that less nodes have been removed. Top 80%
pooling removing 20% of irrelevant nodes shows the best results where the
relative improvement compared to the baseline becomes 11.6% and 13.3%. Through
a comparison of the three different ratios of pooling, we find that discarding
too many nodes worsens the performance (11% pool) but removing a few of the
most irrelevant nodes improves the performance (80% pool). We thus conclude
that the discriminative information is well combined into the majority of
nodes (80% in our case) via the GAT layer.

Comparison of readout operations. In Table 4, we further compare several
readout operations using the RawNet2 baseline. We find that the proposed
graph-based aggregation shows consistent effectiveness, regardless of the
adopted readout approach because all node aggregation methods performed better
than GRU. Among various methods, summation achieved the best result.

Comparison with state-of-the-art systems. In Table 5, we compare the proposed
system with other recent state-of-the-art systems in the literature that adopt
various feature aggregation methods. None of the systems apply data
augmentation techniques. In both input features, the proposed models using the
proposed graph aggregation module demonstrated the best performance, showing
the effectiveness of the proposed approach. Given that we did not use data
augmentation methods, there is room to improve the system using various data
augmentation.

Ablations. Table 6 addresses two ablation experiments using the RawNet2
baseline. First, we explore whether adopting a GRU after the GAT is beneficial
and show the result in the first row. To maintain the overall complexity of
the model, we omitted the last residual block and then placed the GAT and the
GRU in sequence. However, this worsens the performance, as presented in the
first row. In our analysis, this shows that applying GRU directly to the GAT’s
output is not effective, because the GAT’s output is not sequential.

Second, we explore using two GAT layers with different pooling architectures
proposed in [lee2019self] and denote the result in second and third rows. Both
architectures adopt two GAT layers in sequence. The global architecture
concatenates two GAT layers’ output and feed it to the graph pooling layer
followed by readout, whereas hierarchical readout performs graph pooling and
readout after each GAT layers and then add them element-wisely. Through
experiments, we found that both modified pooling architectures did not bring
further improvements.

##  5 Conclusions

In this paper, we proposed a graph feature aggregation method for SV.
Utilizing a GAT, graph pooling layer, and readout operation, we directly
modeled the inter-relationship between entire frame-level features, which is
partially or indirectly utilized in the existing methods. As this is the first
work employing graph neural networks for feature aggregation, we also explored
various configurations to optimize the system. Consistent improvements over
the baselines with different aggregation modules demonstrate the effectiveness
of the proposed approach.
Acknowledgement. We would like to thank Jee-weon Jung at Naver Coporation for
his help with the conceptualization and editing.