Future attractors of Bianchi types II and V cosmologies with massless Vlasov matter

In this paper, we consider a Lorentzian manifold $\bar{M}=I{\times}M$, with an open interval $I\subset \mathbb{R}$ and M being a three-dimensional Riemannian manifold that is also a simply connected three-dimensional Lie group. The Lorentzian and Riemannian metrics are denoted by $\bar{g}$ and g with their associated covariant derivatives $\bar{\nabla }$ and ∇, respectively. Greek indices stand for the spacetime coordinates on $\bar{M}$ whereas Latin indices denote the spatial coordinates on M. We occasionally use the notation |⋅|_g for the norm of a vector field or a symmetric tensor with respect to a given metric g. Finally, C denotes any positive constant which is uniform in the sense that it does not depend on the solution. However, the value of C may change from line to line.

In this section, we briefly introduce the Bianchi models and refer the reader to [36] for a detailed discussion. The Bianchi models are a subclass of SH spacetimes whose isometry group has three-dimensional subgroup M which acts simply transitively on the spacelike orbits that foliate the Lorentzian manifold $\bar{M}$ (cf section 2.1) and are called surfaces of homogeneity. Therefore, in Bianchi models we have $\bar{M}=I{\times}M$. Let M_t be the Lie group M_t := {t} × M which admits a left-invariant frame {E_i } with its dual {Wⁱ }. Moreover, Bianchi models admit a Lie algebra of Killing vector fields K₁, K₂, and K₃ which are tangent to the orbits of the group and satisfy the commutation relation

$$\begin{equation*}\left[{K}_{i},{K}_{j}\right]={C}_{ij}^{\hspace{4pt}k}{K}_{k},\end{equation*}$$

where ${C}_{ij}^{\hspace{4pt}k}$ are the structure constants. Let E₀ be a unit vector field normal to the group orbits. In general, we have $\left[{E}_{\mu },{E}_{\nu }\right]={\gamma }_{\mu \nu }^{\alpha }{E}_{\alpha }$ where ${\gamma }_{\mu \nu }^{\alpha }$ are the commutation functions. Then, one has a natural choice for the time coordinate such that it commutes with the left-invariant frame {E_i }. This has the consequence that the commutation functions are constants and can be identified with the structure constants by a constant linear transformation, i.e.

$$\begin{equation*}\left[{E}_{i},{E}_{j}\right]={C}_{ij}^{\hspace{4pt}k}{E}_{k}.\end{equation*}$$

In this way {E_i } and hence {Wⁱ } are time-independent and the metric in the left-invariant frame reads

$$\begin{equation*}\bar{g}=-\mathrm{d}{t}^{2}+{g}_{ij}\left(t\right){W}^{i}{W}^{j}.\end{equation*}$$

Bianchi models can be classified by classifying the structure constants. The structure constants can be decomposed uniquely as

$$\begin{equation*}{C}_{ij}^{\hspace{4pt}k}={{\epsilon}}_{ij\ell }{n}^{k\ell }+{a}_{i}{\delta }_{j}^{\hspace{4pt}k}-{a}_{j}{\delta }_{i}^{k},\end{equation*}$$

where n^ij is a constant symmetric matrix and a_i are constant. There are two Bianchi classes based on this decomposition; Bianchi class A where a_i = 0 and Bianchi class B where a_i ≠ 0 for at least one i ∈ {1, 2, 3}. Moreover, in Bianchi class A trace of the structure constant vanishes, i.e. ${C}_{ij}^{\hspace{4pt}j}=0$. Finally, by a unimodular transformation we can diagonalize the matrix n^ij as ${n}^{ij}=\mathrm{diag}\left({\hat{n}}_{1},{\hat{n}}_{2},{\hat{n}}_{3}\right)$ (cf e.g. [32]). Therefore, one can classify Bianchi class A only by ${\hat{n}}_{i}$ whereas for Bianchi class B one needs also a_i .

Remark 2.1. The one-forms Wⁱ which are dual to the frame {E_i } satisfy the Maurer–Cartan equation

$$\begin{equation*}\mathrm{d}{W}^{i}=-\frac{1}{2}{C}_{jk}^{\hspace{4pt}i}{W}^{j}\wedge {W}^{k},\end{equation*}$$

which allows us to determine the one-forms in terms of local coordinates for each Bianchi type (cf table 8.2 in [34]). In the Bianchi types studied in the present paper the one-forms have the form

$$\begin{align*}\hfill \left\{\mathrm{d}x-z\enspace \mathrm{d}y,\mathrm{d}y,\mathrm{d}z\right\}& \quad \text{for}\enspace \text{Bianchi}\;\text{type}\;\text{II},\hfill \\ \hfill \left\{\mathrm{d}x,{\mathrm{e}}^{x}\enspace \mathrm{d}y,{\mathrm{e}}^{x}\enspace \mathrm{d}z\right\}& \quad \text{for}\enspace \text{Bianchi}\;\text{type}\;\text{V}.\hfill \end{align*}$$

Remark 2.2. Because a cosmological spacetime is supposed to be the one which admits a compact Cauchy hypersurface, we consider the locally spatially homogeneous spacetimes which constitute a broader class of spacetimes that admit a compact Cauchy hypersurface. The reason is that if we restrict ourselves to the Bianchi models that admit a compact Cauchy hypersurface only Bianchi types I and IX should be considered. A spacetime is locally spatially homogeneous if each point has a neighbourhood that is isometric to an open subset of a spatially homogeneous spacetime. For more technical discussion we refer the reader to [26].

In the frame introduced in previous section, we define (cf [29])

$$\begin{equation*}\rho {:=}{T}^{00},\qquad {j}_{i}{:=}{T}_{i}^{0},\qquad {S}_{ij}{:=}{T}_{ij},\end{equation*}$$

which are the energy density, matter current density, and the spatial part of the energy–momentum tensor, respectively, in terms of the energy–momentum tensor T_αβ which satisfies the Einstein equation G_αβ = 8πT_αβ where G_αβ is the Einstein tensor. We write the Einstein equations in (3 + 1)-decomposition in the Hamiltonian and momentum constraint equations

$$\begin{equation*}\begin{aligned}\hfill R-{k}^{ij}{k}_{ij}+{k}^{2}& =16\pi \rho ,\hfill \\ \hfill {\nabla }^{j}{k}_{ij}& =8\pi {j}_{i},\hfill \end{aligned}\end{equation*}$$

and the evolution equations

$$\begin{equation}{\partial }_{t}{g}_{ij}=-2{k}_{ij},\end{equation} \tag{ 2.1 }$$

$$\begin{equation}{\partial }_{t}{k}_{ij}={R}_{ij}+k{k}_{ij}-2{k}_{i\ell }{k}_{j}^{\ell }-8\pi \left({S}_{ij}-\frac{1}{2}S{g}_{ij}\right)-4\pi \rho {g}_{ij}.\end{equation} \tag{ 2.2 }$$

Here, R_ij is the Ricci tensor associated to the metric g_ij , R = g^ij R_ij is the Ricci scalar, k_ij and k = g^ij k_ij are the second fundamental form of the hypersurface of spatial homogeneity M_t and its trace, respectively, and finally S = g^ij S_ij . Another useful evolution equation is achieved by taking the trace of (2.1b) and using the Hamiltonian constraint

$$\begin{equation}{\partial }_{t}k={k}_{ij}{k}^{ij}+4\pi \left(S+\rho \right).\end{equation} \tag{ 2.2 }$$

Further, we introduce variables that are usually used in SH cosmological models. First, we decompose the second fundamental form as

$$\begin{equation}{k}_{ij}={\sigma }_{ij}-H{g}_{ij}\end{equation} \tag{ 2.3 }$$

to introduce σ_ij which is the trace-free part of k_ij and is referred to as shear tensor in cosmological context, and the Hubble parameter by

$$\begin{equation*}H{:=}-\frac{1}{3}k.\end{equation*}$$

Remark 2.3. Since the future geodesic completeness is known for these models [26], one can choose the time origin without loss of generality as ${t}_{0}=-3{\left[\beta k\left({t}_{0}\right)\right]}^{-1}$, where β > 0 (cf (4.3) below) depends on the specific fixed point under study. This requires to assume that k < 0 for all time. This is indeed the case for the Bianchi models in the present work [25]. Thus, H > 0 for all time and ${t}_{0}={\left[\beta H\left({t}_{0}\right)\right]}^{-1}$.

We next define the Hubble normalized variable

$$\begin{equation*}{{\Sigma}}_{j}^{i}{:=}{H}^{-1}{\sigma }_{j}^{i},\end{equation*}$$

with the evolution equation

$$\begin{equation}{\partial }_{t}{{\Sigma}}_{j}^{i}=-H\left[\left(3-{\partial }_{t}\left({H}^{-1}\right)\right)\left({{\Sigma}}_{j}^{i}-{\delta }_{j}^{i}\right)-\frac{1}{{H}^{2}}{R}_{j}^{i}+\frac{8\pi }{{H}^{2}}{S}_{j}^{i}-\frac{4\pi }{{H}^{2}}\left(S-\rho \right){\delta }_{j}^{i}\right].\end{equation} \tag{ 2.4 }$$

It turns out that the computations will be simplified if one works with

$$\begin{equation*}\begin{aligned}\hfill {{\Sigma}}_{+}& {:=}-\frac{1}{2}\left({{\Sigma}}_{2}^{2}+{{\Sigma}}_{3}^{3}\right),\hfill \\ \hfill {{\Sigma}}_{-}& {:=}-\frac{1}{2\sqrt{3}}\left({{\Sigma}}_{2}^{2}-{{\Sigma}}_{3}^{3}\right),\hfill \end{aligned}\end{equation*}$$

because ${{\Sigma}}_{j}^{i}$ is trace-free. Note that ${{\Sigma}}_{1}^{1}=2{{\Sigma}}_{+}$, ${{\Sigma}}_{2}^{2}=-{{\Sigma}}_{+}-\sqrt{3}{{\Sigma}}_{-}$, and ${{\Sigma}}_{3}^{3}=-{{\Sigma}}_{+}+\sqrt{3}{{\Sigma}}_{-}$. Moreover, in the diagonal case we have $\vert {\Sigma}{\vert }_{g}^{2}=6\left({{\Sigma}}_{+}^{2}+{{\Sigma}}_{-}^{2}\right)$. We also introduce the simplifying notations

$$\begin{equation*}{\Omega}{:=}\frac{8\pi \rho }{3{H}^{2}},\end{equation*}$$

and

$$\begin{equation}q{:=}-1-\frac{\dot {H}}{{H}^{2}},\end{equation} \tag{ 2.5 }$$

which are called density parameter and deceleration parameter, respectively, in cosmological context. Finally, we define

$$\begin{equation*}{\mathfrak{R}}_{j}^{i}{:=}\frac{1}{6{H}^{2}}{R}_{j}^{i},\quad \text{with}\quad \mathfrak{R}{:=}\frac{1}{6{H}^{2}}R.\end{equation*}$$

We give a quick introduction to the Vlasov matter in the context of the general relativity and we refer to [3, 28, 29] for more details.

In this paper, we assume that all particles are massless and modelled by a distribution function which is defined on the mass shell

$$\begin{equation*}\mathcal{P}{:=}\left\{\left(x,p\right):\vert p{\vert }_{\bar{g}}^{2}=0,{p}^{0}{< }0\right\}\subset T\bar{M},\end{equation*}$$

with p = p^μ ∂_μ and p^μ being the canonical coordinates on the tangent bundle of $\bar{M}$. Then, one can associate an energy–momentum tensor to a distribution function $\bar{f}:\mathcal{P}\to \left[0,\infty \right)$ by

$$\begin{equation*}{T}^{\alpha \beta }\left[\bar{f}\right]{:=}{\int }_{{\mathcal{P}}_{x}}\bar{f}\hspace{-1pt}\left(x,p\right){p}^{\alpha }{p}^{\beta }\enspace \mathrm{d}{\mu }_{{\mathcal{P}}_{x}},\end{equation*}$$

where $\mathrm{d}{\mu }_{{\mathcal{P}}_{x}}$ is the Riemannian measure induced on ${\mathcal{P}}_{x}$ by the Lorentzian metric $\bar{g}$ at a given point x, and is given by

$$\begin{equation*}\mathrm{d}{\mu }_{{\mathcal{P}}_{x}}{:=}\frac{\sqrt{\vert \mathrm{det}\enspace \bar{g}\vert }}{-{p}_{0}}\mathrm{d}{p}^{1}\wedge \mathrm{d}{p}^{2}\wedge \mathrm{d}{p}^{3}=\frac{1}{-{p}_{0}}\mathrm{d}{\mu }_{p},\end{equation*}$$

where $\mathrm{d}{\mu }_{p}{:=}\sqrt{\mathrm{det}\enspace g}\enspace \mathrm{d}{p}^{1}\wedge \mathrm{d}{p}^{2}\wedge \mathrm{d}{p}^{3}$. Note that if the metric g is diagonal then −p₀ = p⁰ = |p|_g .

Instead of dealing with $\bar{f}$ we consider $f{:=}\bar{f}{\circ}{\mathrm{p}\mathrm{r}}^{-1}$ which is referred to as distribution function for the remainder of the paper. Here, the projection map is defined by $\mathrm{p}\mathrm{r}:\mathcal{P}\to T\bar{M}$ with (x, p⁰, pⁱ ) ↦ (x, pⁱ ).

Then, the transport equation reads

$$\begin{equation*}{p}^{\mu }{\partial }_{\mu }f-{\bar{{\Gamma}}}_{\mu \nu }^{i}{p}^{\mu }{p}^{\nu }{\partial }_{{p}^{i}}f=0,\end{equation*}$$

where the components of the connection are given by ${\bar{\nabla }}_{{E}_{\alpha }}{E}_{\beta }={\bar{{\Gamma}}}_{\beta \alpha }^{\lambda }{E}_{\lambda }$ and read

$$\begin{equation}{\bar{{\Gamma}}}_{\beta \alpha }^{\lambda }=\frac{1}{2}{\bar{g}}^{\lambda \mu }\left[{E}_{\beta }\left({\bar{g}}_{\mu \alpha }\right)+{E}_{\alpha }\left({\bar{g}}_{\beta \mu }\right)-{E}_{\mu }\left({\bar{g}}_{\alpha \beta }\right)+{\gamma }_{\alpha \beta }^{\sigma }{\bar{g}}_{\mu \sigma }+{\gamma }_{\mu \alpha }^{\sigma }{\bar{g}}_{\beta \sigma }-{\gamma }_{\beta \mu }^{\sigma }{\bar{g}}_{\alpha \sigma }\right].\end{equation} \tag{ 2.6 }$$

Moreover, we assume that f(x, p^a ) has compact support in momentum space at time t₀.

In the case of Bianchi type II symmetry we restrict ourselves to the solutions of the EVS which admit the reflection symmetry, i.e. on the initial data we impose the following conditions (cf [27])

$$\begin{equation*}\begin{aligned}\hfill {f}_{0}\left({p}_{1},{p}_{2},{p}_{3}\right)& ={f}_{0}\left({p}_{1},-{p}_{2},-{p}_{3}\right)={f}_{0}\left(-{p}_{1},-{p}_{2},{p}_{3}\right),\hfill \\ \hfill g\left({t}_{0}\right)& =\mathrm{diag}\left({g}_{11}\left({t}_{0}\right),{g}_{22}\left({t}_{0}\right),{g}_{33}\left({t}_{0}\right)\right),\hfill \\ \hfill k\left({t}_{0}\right)& =\mathrm{diag}\left({k}_{11}\left({t}_{0}\right),{k}_{22}\left({t}_{0}\right),{k}_{33}\left({t}_{0}\right)\right).\hfill \end{aligned}\end{equation*}$$

Under these conditions we have j_i (t₀) = 0 and S_ij (t₀) is diagonal. The reflection symmetry is preserved by the evolution equations for all time, i.e. g_ij , k_ij , and S_ij will remain diagonal and j_i = 0 for all time. In general, it is true that initial data with symmetry yield solutions of the Einstein–Vlasov equations with symmetry (cf [29]).

Let κ_i be the eigenvalues of k_ij with respect to g_ij , i.e. they solve the equation

$$\begin{equation}\enspace \mathrm{det}\left({k}_{j}^{i}-\kappa {\delta }_{j}^{i}\right)=0.\end{equation} \tag{ 2.7 }$$

Then, the generalized Kasner exponents are given by

$$\begin{equation*}{\mathfrak{p}}_{i}{:=}\frac{{\kappa }_{i}}{{\sum }_{\ell }\hspace{2pt}{\kappa }_{\ell }}=-\frac{{\kappa }_{i}}{3H}.\end{equation*}$$

This is a useful concept from which one can extract information, e.g. about the spacetime anisotropy. However, in contrast to the Kasner exponents, the generalized Kasner exponents are, in general, functions of time t. Moreover, they satisfy the first Kasner condition automatically, i.e. ${\sum }_{i}\hspace{2pt}{\mathfrak{p}}_{i}=1$, but not the second one necessarily, i.e. the condition ${\sum }_{i}\hspace{2pt}{\left({\mathfrak{p}}_{i}\right)}^{2}=1$ could be violated. In the case of the Kasner spacetime the generalized Kasner exponents do coincide with the Kasner exponents.

Assuming that the distribution function f, which determines the energy–momentum tensor of the Vlasov matter, is compatible with the Bianchi symmetry, one obtains the transport equation in the following form

$$\begin{equation}{\partial }_{t}f\left(t,{p}^{a}\right)+\left(2{k}_{j}^{i}{p}^{j}-\vert p{\vert }_{g}^{-1}{\omega }_{jk}^{i}{p}^{j}{p}^{k}\right){\partial }_{{p}^{i}}f\left(t,{p}^{a}\right)=0,\end{equation} \tag{ 2.8 }$$

where ${\omega }_{jk}^{i}$ are the Ricci rotation coefficients defined in terms of the structure constants ${C}_{jk}^{\hspace{4pt}i}$ of the Lie algebra by

$$\begin{equation}{\omega }_{jk}^{i}{:=}\frac{1}{2}{g}^{im}\left({C}_{mjk}+{C}_{mkj}-{C}_{jkm}\right)\quad \text{with}\quad {C}_{ijk}{:=}{g}_{km}{C}_{ij}^{\hspace{4pt}m}.\end{equation} \tag{ 2.9 }$$

This follows directly from (2.6). In general, we have ${\omega }_{ij}^{i}=0$ and ${\omega }_{ji}^{i}=-{C}_{ji}^{\hspace{4pt}i}$. Moreover, in Bianchi class A we have ${\omega }_{ji}^{i}=0$.

The relevant components of the energy momentum tensor entering the Einstein equations then are (cf [29])

$$\begin{equation*}\begin{aligned}\hfill \rho & ={\int }_{{\mathbb{R}}^{3}{\backslash}\left\{0\right\}}f\hspace{-1pt}\left(t,{p}^{a}\right)\vert p{\vert }_{g}\enspace \mathrm{d}{\mu }_{p},\hfill \\ \hfill {j}_{i}& ={\int }_{{\mathbb{R}}^{3}{\backslash}\left\{0\right\}}f\hspace{-1pt}\left(t,{p}^{a}\right){p}_{i}\enspace \mathrm{d}{\mu }_{p},\hfill \\ \hfill {S}_{ij}& ={\int }_{{\mathbb{R}}^{3}{\backslash}\left\{0\right\}}f\hspace{-1pt}\left(t,{p}^{a}\right)\frac{{p}_{i}{p}_{j}}{\vert p{\vert }_{g}}\enspace \mathrm{d}{\mu }_{p},\hfill \end{aligned}\end{equation*}$$

where the integrals exclude the element p^a = 0 to assure regularity of the integrand. We drop the domain for simplicity in the following. Note that in the massless case S = ρ. Therefore, the constraint and evolution equations take the form

$$\begin{equation}{\Omega}=1-\frac{1}{6}\vert {\Sigma}{\vert }_{g}^{2}+\mathfrak{R},\end{equation} \tag{ 2.11 }$$

$$\begin{equation}{\partial }_{t}{{\Sigma}}_{j}^{i}=-H\left[\left(2-q\right)\left({{\Sigma}}_{j}^{i}-{\delta }_{j}^{i}\right)-6{\mathfrak{R}}_{j}^{i}+8\pi {H}^{-2}{S}_{j}^{i}\right],\end{equation} \tag{ 2.12 }$$

$$\begin{equation}{\partial }_{t}H=-{H}^{2}\left(1+q\right),\end{equation} \tag{ 2.13 }$$

where

$$\begin{equation*}q=1+\mathfrak{R}+\frac{1}{6}\vert {\Sigma}{\vert }_{g}^{2}.\end{equation*}$$

Equation (2.10b) can be written in terms of Σ₊ and Σ₋ for diagonal components

$$\begin{align*}\hfill {\partial }_{t}{{\Sigma}}_{+}& =-H\left[{{\Sigma}}_{+}\left(2-q\right)+{\Omega}-2\mathfrak{R}+3\left({\mathfrak{R}}_{2}^{2}+{\mathfrak{R}}_{3}^{3}\right)\right]+\frac{4\pi }{H}\left({S}_{2}^{2}+{S}_{3}^{3}\right),\hfill \\ \hfill {\partial }_{t}{{\Sigma}}_{-}& =-H\left[{{\Sigma}}_{-}\left(2-q\right)-\sqrt{3}\left({\mathfrak{R}}_{3}^{3}-{\mathfrak{R}}_{2}^{2}\right)\right]+\frac{4\pi }{\sqrt{3}H}\left({S}_{2}^{2}-{S}_{3}^{3}\right).\hfill \end{align*}$$

Extra evolution equations shall be coupled to the above system for analyzing the future stability that depend on the specific Bianchi model, i.e. the evolution equations of the following variable that are used in diagonal case in Bianchi type II (cf [11, 32]) (no summation on repeated indices in the following equations is assumed)

$$\begin{equation*}{N}_{i}{:=}\frac{{n}_{i}}{H},\quad \text{with}\quad {n}_{i}{:=}{\hat{n}}_{i}\frac{{g}_{ii}}{\sqrt{\mathrm{det}\enspace g}},\end{equation*}$$

based on which one can express the Ricci tensor as

$$\begin{equation*}{R}_{i}^{i}=\frac{1}{2}\left[{n}_{i}^{2}-{\left({n}_{j}-{n}_{k}\right)}^{2}\right],\end{equation*}$$

where (ijk) denotes a cyclic permutation of (123). In the Bianchi type V we consider the evolution equation for $\mathfrak{R}$ instead.

In the following we summarize the concrete form of the variables defined until now for Bianchi types II and V.

2.7.1. Bianchi type II

In this case we have

$$\begin{equation*}{\hat{n}}_{1}=1,\qquad {\hat{n}}_{2}=0,\qquad {\hat{n}}_{3}=0,\end{equation*}$$

and the only non-vanishing components of the structure constants are

$$\begin{equation}{C}_{23}^{\hspace{4pt}1}=1=-{C}_{32}^{\hspace{4pt}1}.\end{equation} \tag{ 2.11 }$$

The components of the Ricci tensor in diagonal case are

$$\begin{equation*}{R}_{1}^{1}=-{R}_{2}^{2}=-{R}_{3}^{3}=-R=\frac{1}{2}{n}_{1}^{2}.\end{equation*}$$

Note that $\mathfrak{R}=-{N}_{1}^{2}/12$. Then, the evolution equations for (Σ₊, Σ₋, N₁) in Bianchi type II with reflection symmetry read

$$\begin{equation}{\partial }_{t}{{\Sigma}}_{+}=-H\left[1+{{\Sigma}}_{+}\left(2-q\right)-\frac{5}{12}{N}_{1}^{2}-\left({{\Sigma}}_{+}^{2}+{{\Sigma}}_{-}^{2}\right)\right]+\frac{4\pi }{H}\left({S}_{2}^{2}+{S}_{3}^{3}\right),\end{equation} \tag{ 2.15 }$$

$$\begin{equation}{\partial }_{t}{{\Sigma}}_{-}=-H{{\Sigma}}_{-}\left(2-q\right)+\frac{4\pi }{\sqrt{3}H}\left({S}_{2}^{2}-{S}_{3}^{3}\right),\end{equation} \tag{ 2.16 }$$

$$\begin{equation}{\partial }_{t}{N}_{1}=-H{N}_{1}\left(4{{\Sigma}}_{+}-q\right).\end{equation} \tag{ 2.17 }$$

2.7.2. Bianchi type V

The Bianchi type V is of class B; this means that a₁ ≠ 0. On the other hand, ${\hat{n}}_{i}=0$ for i = 1, 2, 3. The non-vanishing components of the structure constants are then

$$\begin{equation*}{C}_{1j}^{\hspace{4pt}i}={a}_{1}{\delta }_{j}^{i}=-{C}_{j1}^{\hspace{4pt}i},\quad i,j\in \left\{2,3\right\}.\end{equation*}$$

Without loss of generality we set a₁ = 1. Recall that in this case we do not assume the reflection symmetry. Moreover, we have (cf equation (1.95) in [36])

$$\begin{equation*}{R}_{j}^{i}=\frac{1}{3}R{\delta }_{j}^{i}\quad \text{with}\quad R=-6{a}_{1}{a}^{1}=-6{a}_{1}^{2}{g}^{11}=-6{g}^{11},\end{equation*}$$

because the trace-free part of R_ij vanishes. Then,

$$\begin{equation*}{\mathfrak{R}}_{j}^{i}=\frac{1}{3}\mathfrak{R}{\delta }_{j}^{i}\quad \text{with}\quad \mathfrak{R}=-\frac{1}{{H}^{2}}{g}^{11},\end{equation*}$$

and the evolution equations read (cf (2.10b))

$$\begin{equation}{\partial }_{t}{{\Sigma}}_{j}^{i}=-H\left(2-q\right){{\Sigma}}_{j}^{i}-\frac{8\pi }{H}{S}_{j}^{i};\quad i\ne j,\end{equation} \tag{ 2.18 }$$

$$\begin{equation}{\partial }_{t}{{\Sigma}}_{+}=-H\left[{{\Sigma}}_{+}\left(2-q\right)\right]+\frac{4\pi }{H}\left({S}_{2}^{2}+{S}_{3}^{3}\right)+H{\Omega},\end{equation} \tag{ 2.19 }$$

$$\begin{equation}{\partial }_{t}{{\Sigma}}_{-}=-H{{\Sigma}}_{-}\left(2-q\right)+\frac{4\pi }{\sqrt{3}H}\left({S}_{2}^{2}-{S}_{3}^{3}\right),\end{equation} \tag{ 2.20 }$$

$$\begin{equation}{\partial }_{t}\mathfrak{R}=2H\mathfrak{R}\left(2{{\Sigma}}_{+}+q\right)-\frac{2}{H}\left({{\Sigma}}_{2}^{1}{g}^{12}+{{\Sigma}}_{3}^{1}{g}^{13}\right).\end{equation} \tag{ 2.21 }$$

Note that in this case we pulled out Ω in the evolution equation for Σ₊ and did not exploit the Hamiltonian constraint because the energy density estimate is good enough to control the matter terms.

Remark 2.4. To all systems of evolution equations mentioned above, one should couple the evolution equation for the Hubble parameter H. Moreover, they are constrained by the Hamiltonian constraint (2.10a) and they are coupled to the Vlasov equation through the components of the energy–momentum tensor.

In [11] Calogero and Heinzle showed that the Bianchi type II models having the LRS with massless collisionless particles have a future fixed point with the metric ¹

$$\begin{equation*}{\overline{g}}_{\text{GCS}}=-\mathrm{d}{t}^{2}+{\mathfrak{g}}_{ij}\left(t\right){W}^{i}\otimes {W}^{j},\quad \text{with}\quad \mathfrak{g}\left(t\right)=\mathrm{diag}\left(\frac{8}{9}{t}^{\frac{2}{3}},{t}^{\frac{4}{3}},{t}^{\frac{4}{3}}\right).\end{equation*}$$

The dynamical variables at this fixed point read

$$\begin{equation}H=\frac{5}{9}{t}^{-1},\qquad {{\Sigma}}_{+}=\frac{1}{5},\qquad {{\Sigma}}_{-}=0,\qquad {N}_{1}=\frac{6\sqrt{2}}{5},\qquad q=\frac{4}{5},\end{equation} \tag{ 2.14 }$$

and the generalized Kasner exponents are

$$\begin{equation*}{\mathfrak{p}}_{\text{GCS}}=\left(\frac{1}{5},\frac{2}{5},\frac{2}{5}\right).\end{equation*}$$

The Milne solution is a negatively curved FLRW model with vanishing energy density and pressure whose scale factor is t. The Milne model can be considered as a part of Minkowski spacetime described in comoving coordinates adapted to the wordline of a particle. Its metric in the left-invariant frame is given by

$$\begin{equation*}{\bar{g}}_{M}=-\mathrm{d}{t}^{2}+{\mathbf{g}}_{ij}\left(t\right){W}^{i}\otimes {W}^{j},\quad \text{with}\quad {\mathbf{g}}_{ij}\left(t\right)={t}^{2}{\delta }_{ij}.\end{equation*}$$

The dynamical variables at this fixed point are

$$\begin{equation}H={t}^{-1},\qquad {{\Sigma}}_{j}^{i}=0,\qquad \mathfrak{R}=-1,\qquad q=0,\end{equation} \tag{ 2.15 }$$

and the generalized Kasner exponents, as in other FLRW models, are

$$\begin{equation*}{\mathfrak{p}}_{M}=\left(\frac{1}{3},\frac{1}{3},\frac{1}{3}\right).\end{equation*}$$

We make use of bootstrap argument to prove the main theorems. We refer the reader to [29] for a discussion and to [35] for a rigorous proof of the bootstrap principle. In the following we say that the bootstrap assumptions hold if there are smooth solutions to the systems (2.12) and (2.13) such that on the time interval [t₀, T) for some T > t₀ the following conditions hold:

• (a)  
  In Bianchi type II

  $$\begin{equation*}\begin{aligned}\hfill \left\vert {{\Sigma}}_{+}-\frac{1}{5}\right\vert & {\leqslant}C{\left(1+t\right)}^{-\alpha },\hfill \\ \hfill \left\vert {N}_{1}-\frac{6\sqrt{2}}{5}\right\vert & {\leqslant}C{\left(1+t\right)}^{-\alpha },\hfill \\ \hfill \left\vert {{\Sigma}}_{-}\right\vert & {\leqslant}C{\left(1+t\right)}^{-\alpha },\hfill \end{aligned}\end{equation*}$$

  for 1/5 < α < 4/15, and
• (b)  
  In Bianchi type V

  $$\begin{equation*}\begin{aligned}\hfill \left\vert {{\Sigma}}_{j}^{i}\right\vert & {\leqslant}C{\left(1+t\right)}^{-\alpha },\quad \text{for}\enspace i,j=1,2,3\hfill \\ \hfill \left\vert \mathfrak{R}+1\right\vert & {\leqslant}C{\left(1+t\right)}^{-\alpha },\hfill \end{aligned}\end{equation*}$$

  for 2 < α < 7/2.

We then wish to show that T = ∞. This is done by improving the bootstrap assumptions in the sense that the decay rates would be higher than the bootstrap assumptions after analyzing the system of evolution equations. Then, a continuation criterion finishes the proof.

Note that in what follows α has the values given in this section for different Bianchi types.

Remark 4.1. In the same way one shows the global existence which is already done by Rendall in [25]. Here, we rather examine the stability by the bootstrap argument.

Remark 4.2. The bootstrap assumptions reflect the small-data assumptions. That is, the initial data for the massless EVS are ɛ-close to the GCS and Milne solutions in Bianchi types II and V cases, respectively.

We first linearize the systems of evolution equations from section 2.7 around their future attractors. For this purpose we introduce

$$\begin{equation*}\mathcal{O}\left(\vert \hat{X}\vert \right){:=}\mathcal{O}\left(\vert \hat{{\Sigma}}{\vert }_{g}\right)+\mathcal{O}\left(\vert {\hat{N}}_{1}\vert \right)+\mathcal{O}\left(\vert \hat{\mathfrak{R}}\vert \right),\end{equation*}$$

where the hatted variables are the shifted ones which are defined in the next sections. In Bianchi type V the second term does not exist. We use also $\mathcal{O}\left(\vert \hat{X}{\vert }^{2}\right)$ to denote the sum of all mixed second order terms.

4.2.1. Bianchi type II

We linearize the system (2.12) around the GCS solution by shifting the system to the origin using the variables (cf (2.14))

$$\begin{equation*}{\hat{{\Sigma}}}_{+}{:=}{{\Sigma}}_{+}-\frac{1}{5},\qquad {\hat{{\Sigma}}}_{-}{:=}{{\Sigma}}_{-},\qquad {\hat{N}}_{1}{:=}{N}_{1}-\frac{6\sqrt{2}}{5}.\end{equation*}$$

Hence,

$$\begin{equation*}q=\frac{4}{5}+\frac{2}{5}{\hat{{\Sigma}}}_{+}-\frac{\sqrt{2}}{5}{\hat{N}}_{1}+{\hat{{\Sigma}}}_{+}^{2}+{\hat{{\Sigma}}}_{-}^{2}-\frac{1}{12}{\hat{N}}_{1}^{2},\end{equation*}$$

and

$$\begin{equation}{\partial }_{t}{\hat{{\Sigma}}}_{+}=-H\left(\frac{18}{25}{\hat{{\Sigma}}}_{+}-\frac{24\sqrt{2}}{25}{\hat{N}}_{1}\right)+\frac{4\pi }{H}\left({S}_{2}^{2}+{S}_{3}^{3}\right)+H\cdot \mathcal{O}\left(\vert \hat{X}{\vert }^{2}\right),\end{equation} \tag{ 4.1 }$$

$$\begin{equation}{\partial }_{t}{\hat{{\Sigma}}}_{-}=-\frac{6}{5}H{\hat{{\Sigma}}}_{-}+\frac{4\pi }{\sqrt{3}H}\left({S}_{2}^{2}-{S}_{3}^{3}\right)+H\cdot \mathcal{O}\left(\vert \hat{X}{\vert }^{2}\right),\end{equation} \tag{ 4.2 }$$

$$\begin{equation}{\partial }_{t}{\hat{N}}_{1}=-H\left(\frac{12}{25}{\hat{N}}_{1}+\frac{108\sqrt{2}}{25}{\hat{{\Sigma}}}_{+}\right)+H\cdot \mathcal{O}\left(\vert \hat{X}{\vert }^{2}\right).\end{equation} \tag{ 4.3 }$$

4.2.2. Bianchi type V

We linearize the system (2.13) around the Milne solution by shifting the system to the origin using the variables (cf (2.15))

$$\begin{equation*}{\hat{{\Sigma}}}_{j}^{i}{:=}{{\Sigma}}_{j}^{i},\qquad \hat{\mathfrak{R}}{:=}\mathfrak{R}+1.\end{equation*}$$

Thus,

$$\begin{equation*}q=\hat{\mathfrak{R}}+\frac{1}{6}\vert \hat{{\Sigma}}{\vert }_{g}^{2},\end{equation*}$$

and

$$\begin{equation}{\partial }_{t}{\hat{{\Sigma}}}_{j}^{i}=-2H{\hat{{\Sigma}}}_{j}^{i}-\frac{8\pi }{H}{S}_{j}^{i}+H\mathcal{O}\left(\vert \hat{X}{\vert }^{2}\right);\quad i\ne j,\end{equation} \tag{ 4.4 }$$

$$\begin{equation}{\partial }_{t}{\hat{{\Sigma}}}_{+}=-2H{\hat{{\Sigma}}}_{+}+\frac{4\pi }{H}\left({S}_{2}^{2}+{S}_{3}^{3}\right)+H{\Omega}+H\mathcal{O}\left(\vert \hat{X}{\vert }^{2}\right),\end{equation} \tag{ 4.5 }$$

$$\begin{equation}{\partial }_{t}{\hat{{\Sigma}}}_{-}=-2H{\hat{{\Sigma}}}_{-}+\frac{4\pi }{\sqrt{3}H}\left({S}_{2}^{2}-{S}_{3}^{3}\right)+H\mathcal{O}\left(\vert \hat{X}{\vert }^{2}\right),\end{equation} \tag{ 4.6 }$$

$$\begin{equation}{\partial }_{t}\hat{\mathfrak{R}}=-2H\left(\hat{\mathfrak{R}}+2{\hat{{\Sigma}}}_{+}\right)-\frac{2}{H}\left({{\Sigma}}_{2}^{1}{g}^{12}+{{\Sigma}}_{3}^{1}{g}^{13}\right)+H\mathcal{O}\left(\vert \hat{X}{\vert }^{2}\right).\end{equation} \tag{ 4.7 }$$

Before we give an estimate for the Hubble parameter we need to specify the constant β and hence the time origin t₀ discussed in remark 2.3. We define

$$\begin{equation}\beta {:=}\begin{cases}\frac{9}{5};\quad \text{in}\;\text{Bianchi}\;\text{type}\;\text{II},\quad \hfill \\ 1;\quad \text{in}\;\text{Bianchi}\;\text{type}\;\text{V}.\quad \hfill \end{cases}\end{equation} \tag{ 4.3 }$$

Then, the Hubble parameter has the following estimate.

Lemma 4.3. Assume that the bootstrap assumptions hold. Then,

$$\begin{equation}H={\left(\beta t\right)}^{-1}\left[1+\mathcal{O}\left({t}^{-\alpha }\right)\right],\end{equation} \tag{ 4.4 }$$

with β defined in (4.3).

Proof. The evolution equation for the Hubble parameter (2.10c) can be rewritten as

$$\begin{equation*}{\partial }_{t}\left({H}^{-1}\right)=1+q=\beta +\mathcal{O}\left(\vert \hat{X}\vert \right),\end{equation*}$$

where q is written in terms of the shifted variables given in the previous section. Hence,

$$\begin{equation*}\beta -C{t}^{-\alpha }{\leqslant}{\partial }_{t}\left({H}^{-1}\right){\leqslant}\beta +C{t}^{-\alpha }.\end{equation*}$$

Integrating this inequality on [t₀, t) and using the time origin mentioned in remark 2.3 with β defined in (4.3) completes the proof. □

In this subsection we summarize the estimate of the metric components in different cases. We use the evolution equation for g^ij , i.e.

$$\begin{equation}{\partial }_{t}{g}^{ij}=2H\left({{\Sigma}}_{\ell }^{i}-{\delta }_{\ell }^{i}\right){g}^{\ell j}.\end{equation} \tag{ 4.5 }$$

Lemma 4.4. Assume that bootstrap assumptions hold. Then, for Bianchi type II we have

$$\begin{equation}{g}^{11}={t}^{-\frac{2}{3}}\left[{\mathcal{G}}^{11}+\mathcal{O}\left({t}^{-\alpha }\right)\right],\end{equation} \tag{ 4.6 }$$

$$\begin{equation}{g}^{ii}={t}^{-\frac{4}{3}}\left[{\mathcal{G}}^{ii}+\mathcal{O}\left({t}^{-\alpha }\right)\right];\quad i=2,3,\end{equation} \tag{ 4.7 }$$

whereas in Bianchi type V we have

$$\begin{equation}{g}^{ij}={t}^{-2}\left[{\mathcal{G}}^{ij}+\mathcal{O}\left({t}^{-\alpha }\right)\right].\end{equation} \tag{ 4.8 }$$

Here, ${\mathcal{G}}^{ij}$ are some positive constants.

Proof. Let us start with Bianchi types II. From (4.5) and lemma 4.3 we get

$$\begin{align*}\hfill {\partial }_{t}{g}^{11}& =2H\left(2{{\Sigma}}_{+}-1\right){g}^{11}\hfill \\ \hfill & =2H\left(2{\hat{{\Sigma}}}_{+}+\frac{2}{5}-1\right){g}^{11}\hfill \\ \hfill & {\leqslant}-\frac{6}{5}{\beta }^{-1}{t}^{-1}{g}^{11}+C{t}^{-1-\alpha }{g}^{11},\hfill \end{align*}$$

where we imposed the bootstrap assumptions. Recall that in Bianchi type II β = 9/5. Hence,

$$\begin{equation}{\partial }_{t}\left({t}^{\frac{2}{3}}{g}^{11}\right){\leqslant}C{t}^{-1-\alpha }\left({t}^{\frac{2}{3}}{g}^{11}\right).\end{equation} \tag{ 4.9 }$$

Integrating this inequality on [t₀, t) yields the claim of the lemma.

Similarly, for gⁱⁱ with i = 2, 3 we find

$$\begin{align*}\hfill {\partial }_{t}{g}^{ii}& =2H\left(-{{\Sigma}}_{+}{\pm}\sqrt{3}{{\Sigma}}_{-}-1\right){g}^{ii}\hfill \\ \hfill & =-2H\left(\frac{6}{5}+{\hat{{\Sigma}}}_{+}{\pm}\sqrt{3}{\hat{{\Sigma}}}_{-}\right){g}^{ii}\hfill \\ \hfill & {\leqslant}-\frac{4}{3}{t}^{-1}{g}^{ii}+C{t}^{-1-\alpha }{g}^{ii}.\hfill \end{align*}$$

Integrating this inequality on [t₀, t) gives the estimate for g²² and g³³.

For Bianchi type V we first look at the diagonal components

$$\begin{align*}\hfill {\partial }_{t}{g}^{ii}& =2H\left({{\Sigma}}_{i}^{\hspace{4pt}i}-1\right){g}^{ii}\hfill \\ \hfill & =2H\left({\hat{{\Sigma}}}_{i}^{\hspace{4pt}i}-1\right){g}^{ii}\hfill \\ \hfill & {\leqslant}-2{\beta }^{-1}{t}^{-1}{g}^{ii}+C{t}^{-1-\alpha }{g}^{ii},\hfill \end{align*}$$

for i = 1, 2, 3. Recall that β = 1 in Bianchi type V. Integrating the inequality on [t₀, t) gives ${g}^{ii}={t}^{-2}\left[{\mathcal{G}}^{ii}+\mathcal{O}\left({t}^{-\alpha }\right)\right]$. Then, for the off-diagonal components (i ≠ j) we find

$$\begin{align*}\hfill {\partial }_{t}{g}^{ij}& =2H\left({{\Sigma}}_{i}^{\hspace{4pt}i}-1\right){g}^{ij}+2H\sum\limits _{i\ne \ell }{{\Sigma}}_{\ell }^{i}{g}^{\ell j}\hfill \\ \hfill & {\leqslant}-2H{g}^{ij}+2H\left\vert {\hat{{\Sigma}}}_{\ell }^{i}{g}^{\ell j}\right\vert \hfill \\ \hfill & {\leqslant}-2{t}^{-1}{g}^{ij}+C{t}^{-1-\alpha }{g}^{\text{i}\;\text{j}}+C{t}^{-3-\alpha },\hfill \end{align*}$$

where we used $\vert {g}^{\ell j}\vert {\leqslant}{\left({g}^{\ell \ell }{g}^{jj}\right)}^{1/2}$ since g^ij is a Riemannian metric. Finally, integrating the last inequality on [t₀, t) and using the Grönwall's inequality completes the proof of the lemma. □

We need to estimate the components of the energy–momentum that appear in the systems of the evolution equations. This can be done in several ways. However, it is crucial to achieve the best decay rate such that one can show the stability. Therefore, we begin this subsection by analyzing the behaviour of the momenta of the distribution function f along the characteristic curve Pⁱ (t) of the transport equation (2.8) which reads

$$\begin{equation*}\frac{\mathrm{d}{P}^{i}}{\mathrm{d}t}=2{k}_{j}^{i}{p}^{j}-\vert P{\vert }_{g}^{-1}{\omega }_{jk}^{i}{P}^{j}{P}^{k},\end{equation*}$$

where ${P}^{i}\left(\bar{t}\right)={\bar{p}}^{i}$ for a given $\bar{t}$. Moreover, the characteristic curve P_i (t) takes a simpler form

$$\begin{equation}\frac{\mathrm{d}{P}_{i}}{\mathrm{d}t}=-\vert P{\vert }_{g}^{-1}{\omega }_{mn}^{\ell }{P}_{j}{P}_{k}{g}^{jm}{g}^{kn}{g}_{i\ell },\end{equation} \tag{ 4.10 }$$

where ${P}_{i}\left(\bar{t}\right)={\bar{p}}_{i}$ for a given $\bar{t}$. In the following Pⁱ and P_i indicate that pⁱ and p_i are parameterized by t, respectively. We prove a result about the characteristic curve P_i which will be used later to estimate the energy density.

Lemma 4.5. Assume that bootstrap assumptions hold. Then, for the Bianchi types II and V and for large time t we have

$$\begin{equation}\left\vert {P}_{i}\right\vert ={\mathcal{P}}_{i}+\mathcal{O}\left({t}^{-\alpha }\right);\quad \text{for}\enspace i=1,2,3,\end{equation} \tag{ 4.11 }$$

where ${\mathcal{P}}_{i}$ are positive constants.

Proof. Define the positive function ${F}_{i}{:=}{g}^{ii}{P}_{i}^{2}$ for i = 1, 2, 3 in Bianchi type II. Then, by (4.5) we find

$$\begin{equation}{\partial }_{t}{F}_{i}=2H\left({{\Sigma}}_{i}^{\hspace{4pt}i}-1\right){F}_{i},\end{equation} \tag{ 4.12 }$$

where we used the characteristic equation and the structure constants (2.11) to get

$$\begin{equation*}\frac{\mathrm{d}{P}_{i}}{\mathrm{d}t}=0;\quad \text{for}\enspace i=1,2,3.\end{equation*}$$

Thus, using lemma 4.3, applying the bootstrap assumptions, integrating (4.12) on [t₀, t) and finally comparing the result with that of lemma 4.4 yield the claim of the proof for the Bianchi type II. In a similar vein, we define P² := g^ij P_i P_j for Bianchi type V. Then,

$$\begin{equation*}\frac{\mathrm{d}}{\mathrm{d}t}{P}^{2}={\partial }_{t}{g}^{ij}{P}_{i}{P}_{j}=-2H{P}^{2}+2H{{\Sigma}}^{ij}{P}_{i}{P}_{j},\end{equation*}$$

where we again used the characteristic equation (4.10). Repeating the same argument as before completes the proof. □

Now, we are ready to estimate the components of the energy–momentum tensor.

Proposition 4.6. Assume that bootstrap assumptions hold. Then,

$$\begin{equation}H{S}_{j}^{i}={t}^{-\lambda }\left[{\mathcal{S}}_{j}^{i}+\mathcal{O}\left({t}^{-\alpha }\right)\right],\end{equation} \tag{ 4.13 }$$

where

$$\begin{equation}\lambda =\left\{\begin{array}{c}\frac{4}{3};\quad \text{in}\;\text{Bianchi}\;\text{type}\;\text{II}\;\text{for}\enspace i=j\in \left\{2,3\right\},\quad \\ 3;\quad \text{in}\;\text{Bianchi}\;\text{type}\;\text{V}\;\text{for}\;\enspace i,j\in \left\{1,2,3\right\}.\quad \end{array}\right.\end{equation} \tag{ 4.14 }$$

Here, ${\mathcal{S}}_{j}^{i}$ are positive constants.

Proof. We start with the time derivative of S^ij . Note that pⁱ are time-independent. Then,

$$\begin{equation*}{\partial }_{t}{S}^{ij}=\int {\partial }_{t}f\hspace{-1pt}\left(t,{p}^{a}\right)\frac{{p}^{i}{p}^{j}}{\vert p{\vert }_{g}}\enspace \mathrm{d}{\mu }_{p}+\int f\hspace{-1pt}\left(t,{p}^{a}\right){p}^{i}{p}^{j}{\partial }_{t}\left(\vert p{\vert }_{g}^{-1}\right)\enspace \mathrm{d}{\mu }_{p}+\int f\hspace{-1pt}\left(t,{p}^{a}\right)\frac{{p}^{i}{p}^{j}}{\vert p{\vert }_{g}}{\partial }_{t}\left(\mathrm{d}{\mu }_{p}\right).\end{equation*}$$

Using the Vlasov equation (2.8) for the first term, integrating it by parts, and finally using

$$\begin{align*}{\partial }_{{p}^{i}}\left(\vert p{\vert }_{g}^{-1}\right)& =-\vert p{\vert }_{g}^{-3}{p}_{i},& {\partial }_{t}\left(\vert p{\vert }_{g}^{-1}\right)& ={k}_{ij}{p}^{i}{p}^{j}\vert p{\vert }_{g}^{-3},& {\omega }_{\hspace{2pt}jk}^{i}{p}^{j}{p}^{k}{p}_{i}& =0,\\ {\partial }_{{p}^{i}}\left(\vert p{\vert }_{g}^{-2}\right)& =-2\vert p{\vert }_{g}^{-4}{p}_{i},& {\partial }_{t}\left(\mathrm{d}{\mu }_{p}\right)& =-k\enspace \mathrm{d}{\mu }_{p},& & \end{align*}$$

where we used ${\partial }_{t}\sqrt{g}=-k\sqrt{g}$ in the last equality, we arrive at

$$\begin{equation*}{\partial }_{t}{S}_{j}^{i}={\partial }_{t}\left({g}_{j\ell }{S}^{i\ell }\right)=-4H{S}_{j}^{i}+2H{{\Sigma}}_{\ell }^{i}{S}_{j}^{\ell }+{Y}_{j}^{i}-H{Z}_{j}^{i},\end{equation*}$$

where

$$\begin{align*}\hfill {Y}^{ij}& {:=}2{C}_{ab}^{\hspace{4pt}c}\int f\hspace{-1pt}\left(t,{p}^{k}\right){g}^{a\left(i\right.}{p}^{\left.j\right)}{p}^{b}{p}_{c}\vert p{\vert }_{g}^{-2}\enspace \mathrm{d}{\mu }_{p},\hfill \\ \hfill {Z}_{j}^{i}& {:=}{{\Sigma}}_{m\ell }\int f\hspace{-1pt}\left(t,{p}^{k}\right){p}^{\ell }{p}^{m}{p}^{i}{p}_{j}\vert p{\vert }_{g}^{-3}\enspace \mathrm{d}{\mu }_{p}.\hfill \end{align*}$$

It is readily seen that in reflection symmetric case the following holds

$$\begin{equation*}{Y}_{i}^{i}=0;\quad i=1,2,3.\end{equation*}$$

On the other hand, we have

$$\begin{align*}\hfill {{\Sigma}}_{m\ell }\hspace{2pt}{p}^{m}{p}^{\ell }& ={{\Sigma}}_{1}^{1}{p}^{1}{p}_{1}+{{\Sigma}}_{2}^{2}{p}^{2}{p}_{2}+{{\Sigma}}_{3}^{3}{p}^{3}{p}_{3}\hfill \\ \hfill & =2{{\Sigma}}_{+}{p}^{1}{p}_{1}-{{\Sigma}}_{+}\left({p}^{2}{p}_{2}+{p}^{3}{p}_{3}\right)-\sqrt{3}{{\Sigma}}_{-}\left({p}^{2}{p}_{2}-{p}^{3}{p}_{3}\right)\hfill \\ \hfill & =2{{\Sigma}}_{+}\left[\vert p{\vert }_{g}^{2}-\frac{3}{2}\left({p}^{2}{p}_{2}+{p}^{3}{p}_{3}\right)\right]-\sqrt{3}{{\Sigma}}_{-}\left(\vert p{\vert }_{g}^{2}-{p}^{1}{p}_{1}-2{p}^{3}{p}_{3}\right).\hfill \end{align*}$$

Then, for i = 2, 3 we have

$$\begin{equation}\begin{array}{cc}{\partial }_{t}{S}_{\hspace{2pt}i}^{i}& =-4H{S}_{\hspace{2pt}i}^{i}-2H\left({{\Sigma}}_{+}{\pm}\sqrt{3}{{\Sigma}}_{-}\right){S}_{\hspace{2pt}i}^{i}-H\left(2{{\Sigma}}_{+}-\sqrt{3}{{\Sigma}}_{-}\right){S}_{\hspace{2pt}i}^{i}\\ & +3H{{\Sigma}}_{+}\int f\hspace{-1pt}\left(t,{p}^{a}\right)\frac{{p}^{2}{p}_{2}+{p}^{3}{p}_{3}}{\vert p{\vert }_{g}^{3}}{p}^{i}{p}_{i}\enspace \mathrm{d}{\mu }_{p}\\ & -\sqrt{3}H{{\Sigma}}_{-}\int f\hspace{-1pt}\left(t,{p}^{a}\right)\frac{{p}^{1}{p}_{1}+2{p}^{3}{p}_{3}}{\vert p{\vert }_{g}^{3}}{p}^{i}{p}_{i}\enspace \mathrm{d}{\mu }_{p}.\end{array}\end{equation} \tag{ 4.15 }$$

Note that ${S}_{\hspace{2pt}i}^{i}$ is non-negative. Using ${p}^{2}{p}_{2}+{p}^{3}{p}_{3}{\leqslant}\vert p{\vert }_{g}^{2}$, we estimate the fourth term

$$\begin{align*}3H{{\Sigma}}_{+}\int f\hspace{-1pt}\left(t,{p}^{a}\right)\frac{{p}^{2}{p}_{2}+{p}^{3}{p}_{3}}{\vert p{\vert }_{g}^{3}}{p}^{i}{p}_{i}\enspace \mathrm{d}{\mu }_{p}& {\leqslant}3\vert {\hat{{\Sigma}}}_{+}\vert {S}_{\hspace{2pt}i}^{i}\\ & \quad +\frac{3}{5}\int f\hspace{-1pt}\left(t,{p}^{a}\right)\frac{{p}^{2}{p}_{2}+{p}^{3}{p}_{3}}{\vert p{\vert }_{g}^{3}}{p}^{i}{p}_{i}\enspace \mathrm{d}{\mu }_{p},\end{align*}$$

The fifth term of (4.15) can be estimated similarly. Hence,

$$\begin{align*}\hfill {\partial }_{t}{S}_{\hspace{2pt}i}^{i}& {\leqslant}-H\left[\frac{21}{5}+{\hat{{\Sigma}}}_{+}+\sqrt{3}\left(1{\pm}2\right){\hat{{\Sigma}}}_{-}\right]{S}_{\hspace{2pt}i}^{i}+\left(3\vert {\hat{{\Sigma}}}_{+}\vert +\sqrt{3}\vert {\hat{{\Sigma}}}_{-}\vert \right)H{S}_{\hspace{2pt}i}^{i}\hfill \\ \hfill & {\leqslant}-\frac{7}{3}{t}^{-1}{S}_{\hspace{2pt}i}^{i}+C{t}^{-1-\alpha }{S}_{\hspace{2pt}i}^{i},\hfill \end{align*}$$

where we used (4.4) and imposed the bootstrap assumptions. Consequently,

$$\begin{equation*}{\partial }_{t}\left({t}^{\frac{7}{3}}{S}_{\hspace{2pt}i}^{i}\right){\leqslant}C{t}^{-1-\alpha }\left({t}^{\frac{7}{3}}{S}_{\hspace{2pt}i}^{i}\right).\end{equation*}$$

Integrating this inequality on [t₀, t) yields the claim of the proposition.

For the Bianchi type V we need another strategy because of the non-diagonality. Since f(t₀, p^a ) has compact support and P_i (t) is uniformly bounded for large time t (cf lemma 4.5), there exists a constant C such that

$$\begin{equation*}f\hspace{-1pt}\left(t,{p}_{i}\right)=0;\quad \text{if}\enspace \vert {p}_{i}\vert {\geqslant}C,\end{equation*}$$

for large time t. Now, choose an (spatial) orthonormal frame $\left\{{\tilde {e}}_{i}\right\}$. We denote the quantities in the orthonormal frame by a tilde, i.e. $\tilde {p}{:=}\left({\tilde {p}}_{1},{\tilde {p}}_{2},{\tilde {p}}_{3}\right)$ and $\mathrm{d}{\mu }_{\tilde {p}}{:=}\mathrm{d}{\tilde {p}}_{1}\enspace \mathrm{d}{\tilde {p}}_{2}\enspace \mathrm{d}{\tilde {p}}_{3}$. Then, in the orthonormal frame for Bianchi type V we have (cf [17])

$$\begin{equation*}f\hspace{-1pt}\left(t,\tilde {p}\right)=0;\quad \text{if}\enspace \vert {\tilde {p}}_{i}\vert {\geqslant}C{t}^{-1},\end{equation*}$$

where we used (4.8). Hence, we find

$$\begin{equation*}\rho ={\int }_{0{< }\vert {\tilde {p}}_{i}\vert {\leqslant}C{t}^{-1}}f\hspace{-1pt}\left(t,\tilde {p}\right)\vert \tilde {p}{\vert }_{\delta }\enspace \mathrm{d}{\mu }_{\tilde {p}}.\end{equation*}$$

On the other hand, $f\left(t,\tilde {p}\right)$ remains constant along the characteristics, i.e.

$$\begin{equation*}f\hspace{-1pt}\left(t,\tilde {p}\right){\leqslant}{\Vert}{f}_{0}{\Vert}{:=}\mathrm{sup}\left\{f\hspace{-1pt}\left({t}_{0},\tilde {p}\right):\tilde {p}\in \mathrm{supp}\enspace f\hspace{-1pt}\left({t}_{0},\cdot \right)\right\}.\end{equation*}$$

We thus have

$$\begin{equation*}\rho {\leqslant}C{\Vert}{f}_{0}{\Vert}{t}^{-1}{\int }_{0{< }\vert {\tilde {p}}_{i}\vert {\leqslant}C{t}^{-1}}\mathrm{d}{\mu }_{\tilde {p}}{\leqslant}C{t}^{-4}.\end{equation*}$$

Finally, by lemma 4.3 and the fact that $\vert {S}_{j}^{i}\vert \le \rho$ we finish the proof. □

Remark 4.7. If we drop the reflection symmetry assumption on the distribution function the terms ${Y}_{j}^{i}$ do not vanish for Bianchi type II. Also the components of ${Z}_{j}^{i}$ need more careful analysis in that case. This causes the main difficulty in non-diagonal case. If one could manage to obtain good estimates for these quantities one is able to generalize the result of the present work to the non-diagonal case. However, in Bianchi type VI₀, even in diagonal case, the method used above for Bianchi type II do not provide a good decay rate (the decay rate would be t⁻¹). In fact, the characteristics method gives a better decay rate (the decay rate would be t^−7/4), but it is not sufficient to conclude the stability of Bianchi type VI₀.

Remark 4.8. One could also take advantage of the continuity equation ${\bar{\nabla }}_{\mu }{T}^{\mu 0}=0$. This equation, however, do not lead to a good decay rate for the energy density both in non-diagonal Bianchi types II and VI₀. Nevertheless, one gets the same result for the Bianchi type V. This is due to the geometrical property of the Bianchi type V as discussed in section 1. Note also that this result coincides with the known fact about radiation-dominated FLRW universes with the energy density ρ ∼ a⁴(t) where a(t) is the scale factor.

Before we prove the main part of the theorems 3.1 and 3.2, we first prove the second part of the them. Let ${\hat{\sigma }}_{j}^{i}$ denote the shifted values of ${\sigma }_{j}^{i}$. Then, from (2.3) and (2.7) it follows

$$\begin{equation*}\enspace \mathrm{det}\left[{\hat{\sigma }}_{j}^{i}-\left(\kappa +H-BH\right){\delta }_{j}^{i}\right]=0,\end{equation*}$$

where B is the shifted value which depends on the type of Bianchi symmetry and the components of ${\sigma }_{j}^{i}$. For instance, in Bianchi type II and for ${\sigma }_{1}^{1}$ we have ${B}_{1}=\frac{2}{5}$. Therefore, ${\hat{\kappa }}_{i}{:=}{\kappa }_{i}+\left(1-{B}_{i}\right)H$ are the eigenvalues of ${\hat{\sigma }}_{j}^{i}$ and ${\sum }_{i}{\left({\hat{\kappa }}_{i}\right)}^{2}=\vert \hat{\sigma }{\vert }_{g}^{2}={H}^{2}\vert \hat{{\Sigma}}{\vert }_{g}^{2}$ holds. We finally find

$$\begin{equation*}{\mathfrak{p}}_{i}=-\frac{{\kappa }_{i}}{3H}=\frac{1}{3}\left(1-{B}_{i}\right)+\mathcal{O}\left(\vert \hat{{\Sigma}}{\vert }_{g}\right).\end{equation*}$$

This proves the second parts of the theorems.

In this final subsection we define energy functions for the linearized systems of evolution equations. The energy method was first applied on the Bianchi type I with massless Vlasov matter in [6]. It has two significant advantages. First it makes the calculations simple and second it may be relevant in generalizing the method to the inhomogeneous case.

Note that the energy functions defined below are non-negative and only vanish at the fixed points.

4.7.1. Bianchi type II

For the linearized system (4.1) we define the energy by

$$\begin{equation*}{\mathbf{E}}_{\text{II}}{:=}{c}_{1}{\hat{{\Sigma}}}_{+}^{2}+{\hat{{\Sigma}}}_{-}^{2}+{c}_{2}{\hat{N}}_{1}^{2},\end{equation*}$$

where c_i are some positive constants for i = 1, 2. Then,

$$\begin{align*}\hfill {H}^{-1}{\partial }_{t}{E}_{\text{II}}& =-\frac{36}{25}{c}_{1}{\hat{{\Sigma}}}_{+}^{2}-\frac{12}{5}{\hat{{\Sigma}}}_{-}^{2}-\frac{24}{25}{c}_{2}{\hat{N}}_{1}^{2}+2\sqrt{2}\left(\frac{24}{25}{c}_{1}-\frac{108}{25}{c}_{2}\right){\hat{{\Sigma}}}_{+}{\hat{N}}_{1}\hfill \\ \hfill & \quad +\frac{4\pi }{{H}^{2}}\left({S}_{2}^{2}+{S}_{3}^{3}\right){\hat{{\Sigma}}}_{+}+\frac{4\pi }{\sqrt{3}{H}^{2}}\left({S}_{2}^{2}-{S}_{3}^{3}\right){\hat{{\Sigma}}}_{-}+\mathcal{O}\left(\vert \hat{X}{\vert }^{3}\right).\hfill \end{align*}$$

Now, we introduce a decay inducing positive constant μ such that

$$\begin{equation*}{H}^{-1}{\partial }_{t}{\mathbf{E}}_{\text{II}}=-\mu {\mathbf{E}}_{\text{II}}+{Q}_{\text{II}}+\frac{4\pi }{{H}^{2}}\left({S}_{2}^{2}+{S}_{3}^{3}\right){\hat{{\Sigma}}}_{+}+\frac{4\pi }{\sqrt{3}{H}^{2}}\left({S}_{2}^{2}-{S}_{\enspace 3}^{3}\right){\hat{{\Sigma}}}_{-}+\mathcal{O}\left(\vert \hat{X}{\vert }^{3}\right),\end{equation*}$$

with the quadratic form

$$\begin{align*}\hfill {Q}_{\text{II}}& {:=}\left(\mu -\frac{36}{25}\right){c}_{1}{\hat{{\Sigma}}}_{+}^{2}+\left(\mu -\frac{12}{5}\right){\hat{{\Sigma}}}_{-}^{2}+\left(\mu -\frac{24}{25}\right){c}_{2}{\hat{N}}_{1}^{2}\hfill \\ \hfill & \quad +2\sqrt{2}\left(\frac{24}{25}{c}_{1}-\frac{108}{25}{c}_{2}\right){\hat{{\Sigma}}}_{+}{\hat{N}}_{1}\hfill \end{align*}$$

We wish to choose the constants μ, c₁, and c₂ such that the quadratic form Q_II is negative semidefinite and μ attains its maximum possible value. One then finds that this can be done by choosing

$$\begin{equation*}\mu =\frac{24}{25},\qquad {c}_{2}=\frac{2}{9}{c}_{1},\quad \text{for}\enspace {c}_{1}{ >}0.\end{equation*}$$

Thus, imposing the bootstrap assumptions, using lemma 4.3 and the results of proposition 4.6 we arrive at

$$\begin{align*}\hfill {\partial }_{t}{\mathbf{E}}_{\text{II}}& {\leqslant}-\frac{24}{25}H{\mathbf{E}}_{\text{II}}+\frac{4\pi }{H}\left({S}_{2}^{2}+{S}_{3}^{3}\right){\hat{{\Sigma}}}_{+}+\frac{4\pi }{\sqrt{3}H}\left({S}_{2}^{2}-{S}_{3}^{3}\right){\hat{{\Sigma}}}_{-}+H\mathcal{O}\left(\vert \hat{X}{\vert }^{3}\right)\hfill \\ \hfill & {\leqslant}-\frac{8}{15}{t}^{-1}{\mathbf{E}}_{\text{II}}+C{t}^{-4/3-\alpha }+C{t}^{-1-3\alpha }.\hfill \end{align*}$$

Now, recall that 1/5 < α < 4/15. We integrate the following inequality on [t₀, t)

$$\begin{equation*}{\partial }_{t}\left({t}^{8/15}{\mathbf{E}}_{\text{II}}\right){\leqslant}C{t}^{-4/5-\alpha }+C{t}^{-7/15-3\alpha }\end{equation*}$$

and find

$$\begin{equation*}{t}^{8/15}{\mathbf{E}}_{\text{II}}{\leqslant}{t}_{0}^{8/15}{\mathbf{E}}_{\text{II}}\left({t}_{0}\right)+C{t}_{0}^{1/5-\alpha }+C{t}_{0}^{8/15-3\alpha },\end{equation*}$$

which implies

$$\begin{equation*}{\mathbf{E}}_{\text{II}}{\leqslant}C{t}^{-8/15}.\end{equation*}$$

This in turn means

$$\begin{equation*}\left\vert {\hat{{\Sigma}}}_{+}\right\vert =\mathcal{O}\left({t}^{-4/15}\right),\qquad \left\vert {\hat{{\Sigma}}}_{-}\right\vert =\mathcal{O}\left({t}^{-4/15}\right),\qquad \left\vert {\hat{N}}_{1}\right\vert =\mathcal{O}\left({t}^{-4/15}\right),\end{equation*}$$

which improves the bootstrap assumptions and closes the bootstrap argument, and therefore finishes the proof of theorem 3.1.

4.7.2. Bianchi type V

For the linearized system (4.2) we define

$$\begin{equation*}{\mathbf{E}}_{\text{V}}{:=}\sum\limits _{i\ne j}{\left({\hat{{\Sigma}}}_{j}^{i}\right)}^{2}+{c}_{1}{\hat{{\Sigma}}}_{+}^{2}+{\hat{{\Sigma}}}_{-}^{2}+{c}_{2}{\hat{\mathfrak{R}}}^{2},\end{equation*}$$

for some positive constants c₁ and c₂. Then, similar to the Bianchi type II case we find

$$\begin{equation*}{H}^{-1}{\partial }_{t}{\mathbf{E}}_{\text{V}}=-\mu {\mathbf{E}}_{\text{V}}+{Q}_{\text{V}}+\frac{{C}^{\prime }}{{H}^{2}}{S}_{j}^{i}{{\Sigma}}_{j}^{i}-\frac{4}{{H}^{2}}\hat{\mathfrak{R}}\left({{\Sigma}}_{2}^{1}{g}^{12}+{{\Sigma}}_{3}^{1}{g}^{13}\right)+\mathcal{O}\left(\vert \hat{X}{\vert }^{3}\right),\end{equation*}$$

where μ is a positive constant, C' is some constant and

$$\begin{equation*}{Q}_{\text{V}}{:=}\left(\mu -4\right)\sum\limits _{i\ne j}{\left({\hat{{\Sigma}}}_{j}^{i}\right)}^{2}+\left(\mu -4\right){c}_{1}{\hat{{\Sigma}}}_{+}^{2}+\left(\mu -4\right){\hat{{\Sigma}}}_{-}^{2}+\left(\mu -4\right){c}_{2}\hat{\mathfrak{R}}-8{\hat{{\Sigma}}}_{+}\hat{\mathfrak{R}}.\end{equation*}$$

The quadratic form Q_V is negative semidefinite by choosing

$$\begin{equation*}\mu {\leqslant}4-\frac{4}{\sqrt{{c}_{2}}},\quad {c}_{2}{ >}0,\quad {c}_{1}{c}_{2}{ >}1.\end{equation*}$$

One can choose c₂ sufficiently large such that μ = 4 − 2 with 0 < ≪ 1. Hence, by the smallness assumptions and lemmas 4.3 and 4.4 we obtain

$$\begin{equation*}{\partial }_{t}{\mathbf{E}}_{\text{V}}{\leqslant}-\mu {t}^{-1}{\mathbf{E}}_{\text{V}}+C{t}^{-1-2\alpha }+C{t}^{-3-\alpha }.\end{equation*}$$

Integrating

$$\begin{equation*}{\partial }_{t}\left({t}^{\mu }{\mathbf{E}}_{\text{V}}\right){\leqslant}C{t}^{\mu -1-2\alpha }+C{t}^{\mu -3-\alpha }\end{equation*}$$

on [t₀, t) and keeping in mind that 2 < α < 7/2, one finds

$$\begin{equation*}{t}^{\mu }{\mathbf{E}}_{\text{V}}{\leqslant}{t}_{0}^{\mu }{\mathbf{E}}_{\text{V}}\left({t}_{0}\right)+C{t}_{0}^{\mu -2\alpha }+C{t}_{0}^{\mu -2-\alpha },\end{equation*}$$

which implies

$$\begin{equation*}{\mathbf{E}}_{\text{V}}{\leqslant}C{t}^{-\mu },\end{equation*}$$

and consequently

$$\begin{equation*}\left\vert {\hat{{\Sigma}}}_{j}^{i}\right\vert =\mathcal{O}\left({t}^{-2+{\epsilon}}\right),\qquad \left\vert \hat{\mathfrak{R}}\right\vert =\mathcal{O}\left({t}^{-2+{\epsilon}}\right).\end{equation*}$$

This improves the bootstrap assumptions, closes the bootstrap argument, and hence completes the proof of the theorem 3.2.

No new data were created or analysed in this study.