# The nonlinear stability of n+1 dimensional FLRW spacetimes

Puskar Mondal 111e-mail<EMAIL_ADDRESS>Puskar Mondal† †
Center of Mathematical Sciences and Applications,
Department of Mathematics,
Harvard University, Cambridge, MA02138, USA<EMAIL_ADDRESS>
###### Abstract

We prove nonlinear Lyapunov stability of a family of ‘$n+1$’-dimensional
cosmological models of general relativity locally isometric to the Friedman
Lemaître Robertson Walker (FLRW) spacetimes including a positive cosmological
constant. In particular, we show that the perturbed solutions to the Einstein-
Euler field equations around a class of spatially compact FLRW metrics (for
which the spatial slices are compact negative Einstein spaces in general and
hyperbolic for the physically relevant $n=3$ case) arising from regular Cauchy
data remain uniformly bounded and decay to a family of metrics with constant
negative spatial scalar curvature. To accomplish this result, we employ an
energy method for the coupled Einstein-Euler field equations in constant mean
extrinsic curvature spatial harmonic gauge (CMCSH). In order to handle Euler’s
equations, we construct energy from a current that is similar to the one
derived by Christodoulou [1] (and which coincides with Christodoulou’s current
on the Minkowski space) and show that this energy controls the desired norm of
the fluid degrees of freedom. The use of a fluid energy current together with
the CMCSH gauge condition casts the Einstein-Euler field equations into a
coupled elliptic-hyperbolic system. Utilizing the estimates derived from the
elliptic equations, we first show that the gravity-fluid energy functional
remains uniformly bounded in the expanding direction. Using this uniform
boundedness property, we later obtain sharp decay estimates if a positive
cosmological constant $\Lambda$ is included, which suggests that the
accelerated expansion of the physical universe that is induced by the positive
cosmological constant is sufficient to control the non-linearities in the case
of small data. A few physical consequences of this stability result are
discussed.

##  1 Introduction and previous studies

Since its conception, the FLRW cosmological solution has been used extensively
to model the dynamics of our universe. With the advent of accelerated
expansion of the physical universe, a cosmological constant (the simplest form
of the dark energy) is often included in the model and such a model is
designated as the $\Lambda$CDM model [2, 3]. This model is developed based on
the assumption of global homogeneity and isotropy of the physical universe,
the so-called cosmological principle. However, the astronomical observations
that motivate the cosmological principle are restricted to a fraction
(possibly quite small) of the universe and as such do not provide a hint
towards global properties such as the topology of the universe. Therefore, the
assumption of global homogeneity and isotropy seems to be rather restrictive
and the corresponding topology turns out to be rather limited:
$\mathbb{S}^{3}$ (or its two-fold quotient
$\mathbb{RP}^{3}=\mathbb{S}^{3}/\mathbb{Z}_{2}$), $\mathbb{H}^{3}$ or
$\mathbb{E}^{3}$. Each of these spaces is equipped with its standard metric of
constant positive, zero, or negative sectional curvature, respectively. If one
removes the global assumption and rather imposes a local condition, then there
exist numerous closed manifolds that admit locally homogeneous and isotropic
metrics and that can be constructed as quotients of $\mathbb{E}^{3}$,
$\mathbb{S}^{3}$ and $\mathbb{H}^{3}$ by the discrete, proper and torsion-free
subgroups of their respective isometry groups. However, in the current
article, we will only focus on the universe models with spatial slices that
are closed hyperbolic manifolds. More generally, if one considers an $n+1$
($n\geq 3$) dimensional problem, then one may consider spatial manifolds that
are negative Einstein (for the $n=3$ case these negative Einstein spaces
correspond to hyperbolic manifolds due to Mostow rigidity [4]). Another
motivation to consider the models with negative spatial sectional curvature is
that the recent astronomical observations [5, 6] seem to indicate a
possibility that the spatial universe may have slightly negative curvature.
These compact variants of the most general FLRW solutions may be written as
follows

${}^{n+1}\hat{g}=-dt\otimes dt+a(t)^{2}\gamma_{ij}dx^{i}\otimes dx^{j},$ (1)

where $R[\gamma]_{ij}=-\frac{1}{n}\gamma_{ij},~{}t\in(0,\infty),~{}a(t)$ is
the scale factor satisfying $\partial_{t}a(t)>0$.

One important question that arises in the context of FLRW cosmology is whether
these models are predicted by general relativity? To conclude that they are
would seem to hinge on a proof that the purely theoretical variants of FLRW
models (spatially compact models foliated by compact hyperbolic spaces) are
dynamically stable. In other words, if one were to perturb these solutions in
a suitable function space setting, then the natural question would be whether
the perturbed spacetimes exist for all time as solutions of Einstein’s field
equations. Moreover one would like to understand whether all small data
solutions (in suitable function space settings after subtracting the
background solutions given by (1)) of the Einstein’s field equations remain
bounded for all time in terms of a suitable norm of the initial data and
whether they exhibit asymptotic decay. There are other important physical
motivations such as large-scale structure formation in the universe through
cosmological perturbations [7, 8, 9]. So far, the state-of-the-art work on the
structure formation is limited to linearized perturbation theory, numerical
relativity on flat space background ($\mathbb{E}^{3}$), or $N-$ body
simulations where a Newtonian approximation is applied [10, 11, 13, 12]. Each
of these studies reveals a plethora of interesting physics. Moving one step
further, one would want to study the problem of including the complete non-
linearities of Einstein’s field equations, and deduce whether the large-scale
structure formation is hinted at by the Einstein dynamics. A completely
satisfactory study would be one that considers a large data problem, which is
far from being currently solved. In addition, recently several inconsistencies
of the FLRW models have been proposed by a few recent studies [14, 15, 16]
which require a fully non-linear treatment of the field equations. We will not
address such issues here since they require other information in conjunction
with the nonlinear stability result and therefore remain as a matter for
future study. Overall, it is crucial to study the fully non-linear stability
of these cosmological models in order to draw any physical conclusions
whatsoever.

The story of cosmological perturbation theory is not new. A substantial study
has been carried out on the topic by several authors [17, 18, 19, 20]) since
the earliest investigation by Lifsitz [22]. However, these studies are formal
mode stability analyses and as such are unable to yield sufficient conditions
for a ‘true’ linear stability where one is interested in the asymptotic
behavior of all solutions to the linearized Einstein field equations, not
simply the fixed modes. Even though such true linear stability has been
rigorously established [23], it lacks the nonlinear effects which can lead to
the potential formation of naked singularities or shocks (or both). The most
interesting dynamical behaviors are expected to occur in the fully non-linear
regime. In a sense, the feature that makes Einstein’s field equations (in the
presence or absence of sources) tremendously rich is the non-linearities
present. In the regime of fully nonlinear stability of Einstein’s equations
(may that be vacuum or otherwise), there has been considerable progress in the
last three decades. Global non-linear stability proof of de-Sitter spacetime
by [24] marked the initiation of this progress. This was followed by the proof
of the stability of Minkowski space by [25]. Later [26] provided a simpler
proof utilizing the so-called ‘weak null condition’ of Einstein’s equations in
the spacetime harmonic gauge. In the context of cosmological spacetimes,
Andersson and Moncrief [27] established the stability of the $3+1$ dimensional
Milne universe in constant mean extrinsic curvature spatial harmonic gauge
(CMCSH) and later generalized this result to arbitrary dimensions ($n+1$,
$n\geq 3$) [28]. Following these fundamental studies, numerous studies have
been performed regarding small data global existence and stability problems
associated with Einstein’s equations including various sources [30, 29, 31,
32, 34, 33]. Two of the studies that are most relevant to us are the ones by
Rodnianski and Speck [35] and Speck [36]. Rodnianski and Speck [35] studied
the small data perturbations of the $\mathbb{R}\times\mathbb{T}^{3}$ type FLRW
model with a positive cosmological constant in spacetime harmonic gauge, where
an irrotational perfect fluid model was used as the source term. This
essentially reduced the problem to Einstein’s equations coupled to a nonlinear
scalar field. Later Speck [36] extended this analysis on the same topology and
removed the irrotational condition. Both of these studies utilize the rapid
expansion induced by the positive cosmological constant to avoid pathology
such as shock formation associated with the perfect fluid source. The issue of
shock formation is tremendously important in the non-linear setting (this is
another important difference between linear and non-linear stability analysis
of FLRW spacetimes). [37] proved the formation of shocks in relativistic
perfect fluids evolving on a Minkowski background for arbitrarily small
initial data. This may be attributed to the physical fact that the Minkwoski
space does not posses the expansion property and therefore the purely
geometric dispersion alone is not sufficient to suppress the shock formation.
On a compact spatial topology, the energy concentration (and subsequent
formation of singularities) is more favorable since there are no non-compact
directions through which energy may dissipate (decay). Therefore, the
expansion becomes the only factor that may prevent energy concentration by
non-linearities and the problem essentially turns out to be a competition
between the dissipation caused by expansion and energy concentration by non-
linearities. The global existence or finite-time blow-up then depends on which
of these two opposing forces dominate. For accelerated expansion, it turns out
that the expansion dominates the non-linearities for small data perturbations
yielding global existence.

Despite the large amount of studies that are present in the literature, there
are no non-linear stability results for the FLRW spacetimes to our knowledge
that contain spatially compact slices with constant negative sectional
curvature. In addition, the proof of the stability of
$\mathbb{R}\times\mathbb{T}^{3}$ model by [36] uses the spacetime harmonic
gauge while we work in the constant mean extrinsic curvature spatial harmonic
gauge (CMCSH). Both of our proofs use a form of energy current (an $n+1$
vector field on the spacetime) that was derived by Christodoulou [37]. The
reason behind using the energy current is explained below. We construct the
following current vector field bilinear in perturbations and their higher
order derivatives

$\displaystyle\mathcal{C}^{\mu}\partial_{\mu}:=\left\\{\frac{(\gamma_{ad}-1)}{\gamma_{ad}\rho
N}[1+g(v,v)]^{1/2}\dot{\rho}^{2}+2(\gamma_{ad}-1)\frac{g_{kl}v^{k}\dot{v}^{l}}{N[1+g(v,v)]^{1/2}}\dot{\rho}\right.$
$\displaystyle\left.+\frac{\gamma_{ad}\rho[1+g(v,v)]^{1/2}}{N}\left(-\frac{g(v,\dot{v})^{2}}{1+g(v,v)}+g_{kl}\dot{v}^{k}\dot{v}^{l}\right)\right\\}\partial_{t}$
$\displaystyle+\left(\frac{(\gamma_{ad}-1)v^{i}\dot{\rho}^{2}}{\gamma_{ad}\rho}+2(\gamma_{ad}-1)\dot{v}^{i}\dot{\rho}+\gamma_{ad}\rho
v^{i}(-\frac{g(v,\dot{v})^{2}}{1+g(v,v)}+g_{kl}\dot{v}^{k}\dot{v}^{l})\right)\partial_{i}$
(2)

where $v^{i}\partial_{i}$ is the $n-$velocity of the fluid (i.e., the
projection of the $n+1$-velocity on to a constant $t$ hypersurface), $\rho$ is
the mass energy density, $N$ is the lapse function, $g$ is the induced
Riemannian metric on a $t=$constant hypersurface, and $\gamma_{ad}$ is the
adiabatic index lying in the range $\gamma_{ad}\in(1,\frac{n+1}{n})$. Here
$\dot{v}:=\nabla_{\vec{\alpha}}v$ and
$\dot{\rho}=\nabla_{\vec{\alpha}}(\rho-\bar{\rho})$ are the multi spatial
derivatives of the fields $\rho-\bar{\rho}$ ($\bar{\rho}$ is the background
mass energy density) and $v$. $\nabla_{\vec{\alpha}}$ is roughly defined as
$\nabla_{\vec{\alpha}}=\nabla[\gamma]_{I_{1}}\nabla[\gamma]_{I_{2}}\nabla[\gamma]_{I_{3}}....\nabla[\gamma]_{I_{\alpha}}$
where $\nabla[\gamma]$ is the covariant derivative with respect to a
background Riemannian metric $\gamma$. we show, if
$\gamma_{ad}\in(1,\frac{n+1}{n})$, then
$-\int_{\Sigma}\mathcal{C}^{\mu}\mathbf{n}_{\mu}$ is positive definite and
controls the required norm of the matter density (or perturbed matter density
to be precise) and the fluid’s $n$-velocity field for sufficiently small data,
$\mathbf{n}$ is the unit future pointing time-like vector field orthogonal to
a constant $t$ hypersurface $\Sigma$. While evaluating the time derivative of
the energy through a divergence identity, one needs to compute
$\nabla_{\mu}\mathcal{C}^{\mu}$. In order for the energy argument to close,
$\nabla_{\mu}\mathcal{C}^{\mu}$ should not contain any derivative terms of the
fields $\dot{\rho}$ and $\dot{u}$ (or even if there are derivative terms they
should be of total divergence form so that after integration they disappear).
This cancellation of top order derivatives of the physical fields is obtained
if one uses this special energy current (1). This is the main purpose of
constructing this special energy current since the relativistic Euler’s
equations are not of the diagonal form and therefore an ordinary definition of
energy fails to close the energy argument. Once the energy argument closes,
the remaining lower order terms are estimated using available Sobolev and
multiplication inequalities.

Another technicality that arises in our study is that the trivial solutions
(or background solutions (1)) form a finite dimensional sub-manifold (center
manifold) of the infinite dimensional configuration space if the dimension is
strictly larger than $3+1$. This occurs simply because the trivial solutions
are described by $R[\gamma]_{ij}=-\frac{1}{n}\gamma_{ij}$ which form the
Einstein moduli space after a gauge fixing (CMCSH) is performed. The stability
analysis is slightly tricky in this situation since the perturbations that are
parallel to the moduli space do not have detectable dynamics in the sense of
Einstein’s evolution equations since each perturbed metric $\gamma^{*}$ will
satisfy the criteria $R[\gamma^{*}]_{ij}=-\frac{1}{n}\gamma^{*}_{ij}$ and
therefore will be a trivial solution of the Einstein’s equations. In order to
handle this additional technicality, we invoke the shadow gauge condition
introduced by [28]. Implementation of shadow gauge enables us to explicitly
handle the perturbations that are parallel and perpendicular to the moduli
space.

Our analysis indicates that if the expansion factor $a(t)$ appearing in the
background solutions (1) obeys suitable integrability conditions and the
adiabatic index $\gamma_{ad}$ lies in the range $(1,\frac{n+1}{n})$ ($(1,4/3)$
in the physically relevant ‘3+1’ case; $\gamma_{ad}=1$ and $\gamma_{ad}=4/3$
correspond to a pressure-less fluid or ‘dust’ and a ‘radiation’ fluid,
respectively), then the energy functional remains uniformly bounded by its
initial value and moreover the perturbations decay in the expanding direction.
This is in a sense equivalent to the fact that if the perturbations to the
matter part are restricted within a smaller cone contained in the ‘sound’ cone
(in the tangent space; to be defined later), that is, they do not propagate at
a speed higher than $\sqrt{\frac{1}{n}}$ (in the unit of light speed), then
the stability (Lyapunov or asymptotic to be clarified later) holds. On the
other hand, there exists a seemingly contradictory result in the purely non-
relativistic case of a self-gravitating fluid system. For a self-gravitating
non-relativistic 3-dimensional fluid body (Euler-Poisson system),
Chandrasekhar [39] used a virial identity argument to show that the static
isolated compact solution of the Euler-Poisson system is stable for
$\gamma_{ad}>\frac{4}{3}$. However, if one observes carefully, our result and
Chandrasekhar’s result are not contradictory since the notions of stability
are different in the two contexts. In the context of Chandrasekhar’s argument,
stability is defined by the negativity of total energy i.e., the domination of
gravitational energy (stable $\iff$ gravitationally bound system i.e., the
fluid body does not disperse). Now, of course, a system with strong gravity is
stable in the sense of Chandrasekhar (since total energy is dominated by
gravitational energy and therefore negative) which is nonsensical in our fully
relativistic context (ultra-high gravity could focus and blow up forming a
singularity, indicating instability in our context).

The stability criteria are satisfied by the universe model in which one
includes a positive cosmological constant $\Lambda$ (inclusion of $\Lambda$
forces the integrability conditions on $a(t)$ to hold) and in such case, one
obtains a uniform asymptotic decay (in a suitable sense) which matches with
that of the linear theory. However, turning off the cosmological constant
results in losing the uniform boundedness property simply because one does not
have the desired integrability properties for the scale factor. This may be
attributed to the fact that a positive cosmological constant induced
accelerated expansion wins over the gravitational effect at the level of small
data. However, we do not claim that turning off the cosmological constant
leads to instability but simply we are unable to reach a definite conclusion
with the currently available method. It is worth pointing out that [38] proved
the non-linear asymptotic stability of the Milne universe (devoid of
background matter density) including a dust source (slowly expanding
spacetime). In such a case, the individual dust particles travel along their
respective geodesics and the sound speed vanishes. This intuitively suggests
that one may avoid shock formation in such a scenario with the help of
ordinary polynomial expansion. However, we do not know whether such a result
hold for FLRW spacetimes where the background mass-energy density does not
vanish. We will investigate the borderline cases of $\gamma_{ad}=0$ and
$\gamma_{ad}=\frac{n+1}{n}$ in the future using special techniques.

###  1.1 Notations

The ‘$n+1$’ dimensional $C^{\infty}$ spacetime manifold is denoted by
$\hat{M}$. We are interested in spacetime manifolds $\hat{M}$ of globally
hyperbolic type i.e., the topology is the product topology
$\mathbb{R}\times\Sigma$, where $\Sigma$ denotes an $n$-dimensional closed
spatial manifold diffeomorphic to a Cauchy hypersurface (i.e., every
inextendible causal curve intersects $\Sigma$ exactly once). We denote the
space of Riemannian metrics on $\Sigma$ by $\mathcal{M}_{\Sigma}$. We
explicitly work in the $L^{2}$ (with respect to a background metric) Sobolev
space $H^{s}$ for $s>\frac{n}{2}+2$. The general $L^{p}$ Sobolev space of
order $s$ is denoted by $W^{s,p}$. We will consider tensors as sections of a
suitable vector bundle $\mathcal{B}$ over $(\Sigma,g)$, where
$g\in\mathcal{M}_{\Sigma}$ is the dynamical metric that solves the Einstein’s
equations written in $n+1$ form. The $L^{2}$ inner product on a fibre of such
a rank-$k$ bundle is defined as

$\displaystyle\langle
u|v\rangle_{L^{2}}:=\int_{\Sigma}u_{i_{1}i_{2}i_{3}....i_{k}}v_{j_{1}j_{2}j_{3}....j_{k}}\gamma^{i_{1}j_{1}}\gamma^{i_{1}j_{2}}\gamma^{i_{3}j_{3}}....\gamma^{i_{k}j_{k}}\mu_{g},$
(3)

where $\gamma\in\mathcal{M}_{\Sigma}$ is a $C^{\infty}$ background metric
lying in a small enough neighbourhood of $g$ in a suitable space and
$\mu_{g}=\sqrt{\det(g_{ij})}dx^{1}\wedge dx^{2}\wedge
dx^{3}\wedge.........\wedge dx^{n}$ is the volume form associated with the
metric $g$. The standard norms are defined naturally as follows

$\displaystyle||u||_{L^{2}}:=\langle u|u\rangle^{1/2}_{L^{2}}$ (4)
$\displaystyle||u||_{L^{\infty}}:=\sup_{\Sigma}(u_{i_{1}i_{2}i_{3}....i_{k}}u_{j_{1}j_{2}j_{3}....j_{k}}\gamma^{i_{1}j_{1}}\gamma^{i_{1}j_{2}}\gamma^{i_{3}j_{3}}....\gamma^{i_{k}j_{k}})^{1/2}$
(5)

and so on. We adapt the following definition of a twisted rough Laplacian
acting on sections of $\mathcal{B}^{k}$

$\displaystyle\Delta^{\gamma}_{g}u:=-\frac{1}{\mu_{g}}\nabla[\gamma]_{i}(g^{ij}\mu_{g}\nabla[\gamma]_{j}u).$
(6)

This Laplacian is self-adjoint with respect to the following inner product on
derivatives

$\displaystyle\langle\nabla[\gamma]u|\nabla[\gamma]v\rangle_{L^{2}}:=\int_{\Sigma}g^{ij}\nabla[\gamma]_{i}u_{i_{1}i_{2}....i_{k}}\nabla[\gamma]_{j}v_{j_{1}j_{2}....j_{k}}\gamma^{i_{1}j_{1}}\gamma^{i_{1}j_{2}}....\gamma^{i_{k}j_{k}}\mu_{g},$
(7)

where $\mu_{g}$ is the volume form $\sqrt{\det(g)}dx^{1}\wedge dx^{2}\wedge
dx^{3}\wedge....\wedge dx^{n}$. The covariant derivative with respect to the
dynamical metric $g$ will be denoted simply by $\nabla$ while that with
respect to $\gamma$ will be denoted by $\nabla[\gamma]$. We denote the
connection coefficients of metric $g$ and $\gamma$ by $\Gamma[g]$ and
$\Gamma[\gamma]$, respectively. The Riemann and Ricci curvatures are also
denoted similarly. We define the ordinary Laplacian of $g$ in the following
way such that it has non-negative spectrum i.e.,

$\displaystyle\Delta_{g}\equiv-g^{ij}\nabla_{i}\nabla_{j}.$ (8)

We also define a Lichnerowicz type Laplacian $\mathcal{L}_{g,\gamma}$ acting
on sections of a symmetric rank-2 vector bundle, which will play a crucial
role in our analysis

$\displaystyle\mathcal{L}_{g,\gamma}u_{ij}=\Delta^{\gamma}_{g}u_{ij}-2R[\gamma]_{i}~{}^{k}~{}_{j}~{}^{l}u_{kl}.$
(9)

The spectra of this operator is assumed to be non-negative [28].

Lastly, for functions $f(t)\in\mathbb{R}_{\geq 0}$ and
$g(t)\in\mathbb{R}_{\geq 0}$ $f(t)\lesssim g(t)$ means $f(t)\leq Cg(t)$ and
$f(t)\gtrsim g(t)$ means $f(t)\geq Cg(t)$, and $f(t)\approx g(t)$ implies
$C_{1}g(t)\leq f(t)\leq C_{2}g(t)$ for $0<C_{1},C_{2},C<\infty$. The involved
constants may depend only on a fixed background geometry. The spaces of
symmetric covariant 2-tensors and vector fields on $\Sigma$ are denoted by
$S^{0}_{2}(\Sigma)$ and $\mathfrak{X}(\Sigma)$, respectively. By $A^{tr}$ and
$A^{TT}$, we will denote symmetric trace-less and transverse-trace-less
2-tensors, respectively.

###  1.2 Motivation and Overview

In this section we describe our motivation for this study, state a rough
version of the main theorem and discuss the overall structure of the article.
As we have mentioned in the introduction before, the obvious motivation is to
understand whether the purely theoretical FLRW solutions are predicted by
general relativity. If the overall characteristics of this model still persist
after adding perturbations (in a suitable sense), then we may hope that this
purely theoretical model is physically realized. This question is directly
related to several aspects of modern cosmology and therefore very crucial to
answer.

In addition to the above motivation, we want to study the stability properties
of this particular variant of FLRW models (1) because several recent studies
[40, 41, 42] have suggested that the observed local homogeneity and isotropy
(almost) of the physical universe may be derived dynamically from the
Einstein’s equations (including physically reasonable matter sources). This
idea is supported by the temporal behavior of a certain weak Lyapunov
functional (which takes the form of a re-scaled volume functional of the
spatial manifold; the spatial manifold is assumed to be of negative Yamabe
type for a technical reason, see [42] for more detail) that controls the
lowest norm of the gravity and matter fields. Such a Lyapunov functional
monotonically decays in the expanding direction and achieves its infimum on
solutions that are isometric to our background solutions (1). In that limit,
the Lyapunov functional becomes a topological invariant of the manifold (the
so-called $\sigma$-constant that is closely related to the Yamabe invariant).
Applying a well-established theorem from geometry [43], one then observes that
the value of the Lyapunov functional (and hence the volume of the physical
universe) is dominated by the hyperbolic components that are present in the
spatial manifold. Since the compact hyperbolic manifolds are locally
homogeneous and isotropic, this result tends to indicate that there is a
dynamical mechanism at work within the framework of Einstein flow that drives
the physical universe towards an asymptotic state which is characterized by a
locally homogeneous and isotropic spatial metric.

However, since the Lyapunov functional only controls a very rough norm
($H^{1}\times L^{2}$ of (metric, second fundamental form)) of the data, one is
unable to draw a definite conclusion without studying the long time behavior
of the classical solutions (or solutions that lie in the higher Sobolev
spaces) of the Einstein field equations. This is due to the fact that
gravitational energy can concentrate and blow up before the weak Lyapunov
functional ever achieves its infimum. In the presence of matter sources, there
are additional complications since matter-energy could blow up as well. The
prototypical example of such matter singularity is present in the current
context. The matter source in the context of FLRW cosmology is the perfect
fluid and it is known that the relativistic perfect fluid forms shock
singularities in finite time regardless of the size of the initial data on a
Minkowski background ([37]). These are the associated issues and therefore in
order to understand the extent to which the weak Lyapunov functional decays to
its infimum (or maybe obstructed from doing so), it is necessary to study the
global properties of the solutions. As it happens, in the current context we
would only want to study a ‘small data’ problem, which is expected to provide
a level of insight into several physics questions associated with the
asymptotic behavior of the physical universe. Without going into the technical
details, we state the main theorem and discuss a few physical implications.
Let us denote the gravitational degrees of freedom by $(g,k)$, where $g$ and
$k$ are the first and second fundamental forms of the constant time
hypersurfaces foliating the spacetime. The fluid’s mass-energy density and the
spatial velocity fields are denoted by $\rho$ and $v$, respectively. A rough
version of the main stability theorem is as follows
Rough sketch of the main theorem: Let $(g_{0},k^{tr}_{0},\rho_{0},v_{0})\in
H^{s}\times H^{s-1}\times H^{s-1}\times H^{s-1},~{}s>\frac{n}{2}+2$ be the
initial data for the re-scaled Einstein-Euler-$\Lambda$ system with the
cosmological constant $\Lambda>0$. Let us also consider that the Cauchy data
for the re-scaled background solutions are
$(g_{B.G},k^{tr}_{B.G},\rho_{B.G},v_{B.G})$ which satisfy
$R[g_{B.G}]_{ij}=-\frac{1}{n}(g_{B.G})_{ij},k^{tr}_{B.G}=0,\rho_{B.G}=C_{\rho},v^{i}_{B.G}=0$,
$C_{\rho}$ is a constant. Now assume that the following open smallness
condition holds
$||g_{0}-g_{B.G}||_{H^{s}}+||k^{tr}_{0}-k^{tr}_{B.G}||_{H^{s-1}}+||\rho_{0}-\rho_{B.G}||_{H^{s-1}}+||v_{0}-v_{B.G}||_{H^{s-1}}<\epsilon$,
$\epsilon>0$ is chosen sufficiently small. If
$t\mapsto(g(t),k^{tr}(t),\rho(t),v(t))$ is the maximal development of the
Cauchy problem for the re-scaled Einstein-Euler-$\Lambda$ system in constant
mean extrinsic curvature spatial harmonic gauge (CMCSH) (50-55) with initial
data $(g_{0},k^{tr}_{0},\rho_{0},v_{0})$, then the following holds in the
limit of infinite time

$\displaystyle\lim_{t\to\infty}(g(t),k^{tr}(t),\rho(t),v(t))=(\gamma^{{\dagger}},0,\rho^{{}^{\prime}},0),$
(10)

where $\gamma^{{\dagger}}$ satisfies $R[\gamma^{{\dagger}}]=-1$ and
$\rho^{{}^{\prime}}$ depends on the initial re-scaled density $\rho_{0}$ and
$\partial_{i}\rho^{{}^{\prime}}\neq 0$ in general.

Let us now explain the main theorem putting aside the technical details. The
statement essentially deals with the asymptotic stability of the solutions
that are characterized by constant negative scalar curvature ($-1$ to be
precise in this occasion) and sufficiently close to the FLRW solutions. This
corresponds to a Lyapunov stability of the FLRW solutions. In other words, the
solutions of the linearized Einstein-Euler-$\Lambda$ equations about the
background FLRW solutions (1) decay to nearby solutions which have constant
negative scalar curvature while the background FLRW spacetimes (1) are
described by the property that the spatial metric $g_{B.G}$ is negative
Einstein i.e., $R[g_{B.G}]_{ij}=-\frac{1}{n}(g_{B.G})_{ij}$. Now, the space of
negative Einstein metrics is a subspace of the space of metrics with constant
negative scalar curvature. Simple use of the triangle inequality shows that
the asymptotic state should be the space of constant negative scalar curvature
metrics sufficiently close to and containing the background FLRW solutions.
Since the latter is a subspace of the former, one may fine-tune the initial
conditions to obtain the asymptotic solution to be exactly of FLRW type (i.e.,
the solution has a spatial metric that is negative Einstein). But such a
procedure proves to be too restrictive and physically undesirable.

If for the moment we focus on the physical $3+1$ case, then ideally we want
the asymptotic solution to have constant negative sectional curvature (i.e.,
be hyperbolic or negative Einstein by Mostow rigidity) not just constant
negative scalar curvature and have a spatially uniform background energy
density. In other words, if one perturbs an FLRW solution, it does not
generically come back to an FLRW solution in infinite time. Now notice that
the local spatial homogeneity and isotropy criteria required constancy of the
sectional curvature and as such in 3 spatial dimensions, a negative Einstein
metric automatically has constant sectional curvature. However, the space of
metrics with constant negative scalar curvature is a much larger space and
does not necessarily exhibit the local spatial homogeneity and isotropy
criteria (only a subspace that is described by hyperbolic metrics does). This
result is not completely satisfactory but on the other hand, physically
expected as these inhomogeneous and anisotropic characteristics may be
attributed to the structure formation (see [22, 44] about the role of
cosmological perturbation theory in structure formation). This is indeed a
remarkable fact that indicates there is a dynamical mechanism at work within
the Einstein-Euler-$\Lambda$ flow that naturally drives the physical universe
to an anisotropic and inhomogeneous state (small deviation from homogeneity
and isotropy in a suitable sense) leading to cosmological structure formation.
On the other hand, a Lyapunov functional that controls a very rough norm of
the data may be unable to detect any structure that is smoothed out by the
expansion of the universe. However, at the moment, we do not claim that such a
mechanism of structure formation is yet so compelling since we can not prove
that these background solutions are not asymptotically stable and approach a
homogeneous and isotropic state.

We note that [45] reported this behavior of the FLRW spacetimes as well
considering ‘dust’ as the matter source through a rather heuristic argument.
They constructed a finite-dimensional dynamical system by spatially averaging
out the inhomogeneities and showed that FLRW cosmologies are unstable in some
relevant cases: averaged models are driven away from the background through
structure formation and accelerated expansion. The phenomenon of accelerated
expansion driving the solutions away from the background is previously noted
in [46] which dealt with a particular case of vacuum gravity (background
solutions were described by the compact hyperbolic manifolds) including a
positive cosmological constant. The solutions decay asymptotically to the
nearby ones described by constant negative scalar curvature. On the other
hand, in the same problem if the cosmological constant is turned off (i.e.,
the vacuum energy is absent), then the perturbed solutions do come back to the
background solution asymptotically [27, 28]. Similar phenomena of asymptotic
stability are observed in the study of [38], where the background solutions
once again are simply the Milne universe that does not contain a non-zero
energy density. In order to reach a definite conclusion about forming large-
scale structures by non-linear perturbations of a background with non-
vanishing energy density would require further analysis including radiation
sources and possibly more observational inputs.

Now we discuss the overview of the article. Section 2 provides the necessary
detail about the formulation of the gauge fixed Einstein-Euler field equations
including a positive cosmological constant $\Lambda$, the re-scaling of the
field equations, and the background solutions. Here we state a local well-
posedness theorem (local existence, uniqueness, and Cauchy stability of the
solution as well as the conservation of gauge conditions and constraints along
the solution curve) for the coupled Einstein-Euler system in CMCSH gauge. This
yields an ‘in time’ continuation criteria for the solutions and plays a
crucial role in the final step of the global existence proof. We do not
describe in detail the properties of the background solutions rather only the
key features that are essential to our analysis. Since we are studying the
general ‘$n+1$’ dimensional case, the trivial solutions need to be handled
carefully since they form a finite-dimensional manifold for $n>3$ (Einstein
moduli space). Therefore, we briefly describe the center manifold of the re-
scaled dynamics when the Einstein moduli space is non-trivial and the
kinematics of the perturbations (this is crucial) under the assumption of a
shadow gauge condition (introduced by [28]).

Section $3$ is one of the most important sections of the article. In this
section, we construct the energy current for the fluid and from which we
construct a suitable energy. We prove that under the physical assumptions,
this energy is positive definite. Simultaneously, we construct a wave equation
type energy functional for the gravitational degrees of freedom and show that
the total energy of the system controls the suitable norms of the fluid and
gravity degrees of freedom for sufficiently small data (coercive property). In
this section, we explicitly compute the time evolution of the fluid energy
functional as well. In order to do so, we first obtain a set of commuted Euler
equations for the fluid degrees of freedom. Using these commuted equations, we
show that the spacetime divergence of the energy current
$\mathcal{C}^{\mu}\partial_{\mu}$ does not contain any contribution from the
principal terms. This is an essential feature of a hyperbolic system. In a
sense, the relativistic Euler equations are examples of hyperbolic equations
and as such theoretically derivable from a Lagrangian (see [47, 49, 48, 50]).
The gravity part of the total energy is relatively easy to handle since the
application of CMCSH gauge cancels the problematic terms obstructing the
ellipticity of the map $g\mapsto Ricci[g]$ therefore also making the Einstein
evolution equations hyperbolic (essentially we assumed this property while
stating a local existence theorem).

While computing the time derivative of the total energy, we will encounter
terms involving the lapse function and the shift vector field which therefore
need to be estimated by other means. However, the applications of the CMCSH
gauge condition yields elliptic equations for the lapse function and the shift
vector field that are sourced by the matter degrees of freedom. Therefore, we
obtain the necessary elliptic estimates and use those to control the time
derivative of the energy. One subtlety is that apart from estimating the lapse
and the shift vector fields, we also need to estimate the time derivative of
the lapse function. We complete these estimates in section $4.2$. Using the
necessary estimates for the gauge variables, in section $4.3$ and $4.4$, we
obtain a differential inequality for the energy, which upon integration and
using the smallness of the data yields the uniform boundedness property of the
energy (essentially closing a bootstrap argument). Utilizing the local
existence theorem, we immediately then obtain a global existence theorem.
Using the uniform boundedness property, we prove decay separately. We close
this section by stating the fully rigorous version of the main theorem and
sketching a proof of the geodesic completeness of the perturbed spacetimes.

##  2 Gauge fixed field equations and re-scaling

###  2.1 Complete Einstein-Euler system and gauge fixing

The ADM formalism splits the spacetime described by an ‘$n+1$’ dimensional
Lorentzian manifold $\tilde{M}$ into $\mathbb{R}\times M$. Here each level set
$\\{t\\}\times M$ of the time function $t$ is an orientable n-manifold
diffeomorphic to a Cauchy hypersurface (assuming the spacetime admits a Cauchy
hypersurface) and equipped with a Riemannian metric. Such a split may be
executed by introducing a lapse function $N$ and shift vector field $X$
belonging to suitable function spaces on $\tilde{M}$ and defined such that

$\displaystyle\partial_{t}$ $\displaystyle=$ $\displaystyle N\hat{n}+X,$ (11)

where $t$ and $\hat{n}$ are time and a hypersurface orthogonal future directed
timelike unit vector i.e., $\hat{g}(\hat{n},\hat{n})=-1$, respectively. The
above splitting writes the spacetime metric $\hat{g}$ in local coordinates
$\\{x^{\alpha}\\}_{\alpha=0}^{n}=\\{t,x^{1},x^{2},....,x^{n}\\}$ as

$\displaystyle\hat{g}$ $\displaystyle=$ $\displaystyle-N^{2}dt\otimes
dt+g_{ij}(dx^{i}+X^{i}dt)\otimes(dx^{j}+X^{j}dt)$ (12)

and the stress-energy tensor as

$\displaystyle\mathbf{T}$ $\displaystyle=$ $\displaystyle
E\mathbf{n}\otimes\mathbf{n}+2\mathbf{J}\odot\mathbf{n}+\mathbf{S},$ (13)

where $\mathbf{J}\in\Omega^{1}(M)$, $\mathbf{S}\in S^{0}_{2}(M)$, and $A\odot
B=\frac{1}{2}(A\otimes B+B\otimes A$). Here, $\Omega^{1}(M)$ and
$S^{0}_{2}(M)$ are the space of one form fields and the space of symmetric
covariant 2-tensor fields, respectively. Here
$E:=\mathbf{T}(\mathbf{n},\mathbf{n})$ is the energy density observed by a
time-like observer with $n+1$-velocity $\mathbf{n}$,
$\mathbf{J}_{i}=-\mathbf{T}(\partial_{i},\mathbf{n})$ is the momentum density,
$\mathbf{J}_{i}=-\mathbf{T}(\mathbf{n},\partial_{i})$ is the energy flux
density, and $\mathbf{S}_{ij}=\mathbf{T}(\partial_{i},\partial_{j})$ is the
momentum flux density (with respect to the chosen constant $t$ hypersurface
$M$). The choice of a spatial slice in the spacetime leads to consideration of
the second fundamental form $k_{ij}$ which describes how the slice is curved
in the spacetime. The trace of the second fundamental form ($tr_{g}k=\tau$) is
the mean extrinsic curvature of the slice, which will play an important role
in the analysis. Under such decomposition, the Einstein equations

$\displaystyle R_{\mu\nu}-\frac{1}{2}Rg_{\mu\nu}+\Lambda g_{\mu\nu}$
$\displaystyle=$ $\displaystyle T_{\mu\nu}$ (14)

take the form ($8\pi G=c=1$)

$\displaystyle\partial_{t}g_{ij}$ $\displaystyle=$
$\displaystyle-2Nk_{ij}+L_{X}g_{ij},$ (15) $\displaystyle\partial_{t}k_{ij}$
$\displaystyle=$ $\displaystyle-\nabla_{i}\nabla_{j}N+N\\{R_{ij}+\tau
k_{ij}-2k_{ik}k^{k}_{j}-\frac{1}{n-1}(2\Lambda-S+E)g_{ij}$
$\displaystyle-\mathbf{S}_{ij}\\}+L_{X}k_{ij}$

along with the constraints (Gauss and Codazzi equations)

$\displaystyle R(g)-|k|^{2}+\tau^{2}$ $\displaystyle=$ $\displaystyle
2\Lambda+2E,$ (16) $\displaystyle\nabla_{j}k^{j}_{i}-\nabla_{i}\tau$
$\displaystyle=$ $\displaystyle-\mathbf{J}_{i},$ (17)

where $S=g^{ij}\mathbf{S}_{ij}$. The vanishing of the covariant divergence of
the stress energy tensor i.e., $\nabla_{\nu}T^{\mu\nu}=0$ is equivalent to the
continuity equation and equations of motion of the matter

$\displaystyle\frac{\partial E}{\partial t}$ $\displaystyle=$ $\displaystyle
L_{X}E+NE\tau-N\nabla_{i}\mathbf{J}^{i}-2\mathbf{J}^{i}\nabla_{i}N+N\mathbf{S}^{ij}k_{ij},$
(18) $\displaystyle\frac{\partial\mathbf{J}^{i}}{\partial t}$ $\displaystyle=$
$\displaystyle
L_{X}\mathbf{J}^{i}+N\tau\mathbf{J}^{i}-\nabla_{j}(N\mathbf{S}^{ij})+2Nk^{i}_{j}\mathbf{J}^{j}-E\nabla^{i}N.$
(19)

Here we want to study the Einstein-Euler coupled system. Let the ‘$n+1$’
velocity field of a perfect fluid be denoted by $\mathbf{u}$ which satisfies
the normalization condition $\hat{g}(\mathbf{u},\mathbf{u})=-1$. One may for
convenience decompose the $n+1$ velocity $\mathbf{u}$ into its component
parallel and perpendicular to constant $t$ hypersurface $M$ as follows

$\displaystyle\mathbf{u}=v-\hat{g}(\mathbf{u},\mathbf{n})\mathbf{n}$ (20)

where $v$ is a vector field parallel to the spatial manifold $M$ i.e.,
$v\in\mathfrak{X}(M)$. Importantly note that $\mathbf{u}_{i}=v_{i}$ but
$\mathbf{u}^{i}\neq v^{i}$ unless the shift vector field $X$ vanishes. The
stress energy tensor for a perfect fluid with $n+1$ velocity field
$\mathbf{u}$ reads

$\displaystyle\mathbf{T}=(P+\rho)\mathbf{u}\otimes\mathbf{u}+P\hat{g},$ (21)

where $P$ and $\rho$ are the pressure and the mass energy density,
respectively and $\hat{g}$ is the spacetime metric. After projecting onto the
spatial manifold $M$, several components of the stress energy tensor may be
computed as follows

$\displaystyle
E=(P+\rho)(N\mathbf{u}^{0})^{2}-P,\mathbf{J}^{i}=(P+\rho)\sqrt{1+g(v,v)}v^{i},$
(22) $\displaystyle\mathbf{S}=(P+\rho)v\otimes v+Pg.$

However, notice that the field equations do not close with the information of
the stress-energy tensor alone, and therefore one needs an equation of state
relating pressure $P$ and mass-energy density $\rho$. The choice of the
equation of state is non-trivial and there are numerous studies about the
equation of state alone in the cosmology literature (e.g., [51, 52, 53]). In
the standard model of cosmology ($\Lambda$CDM), one frequently uses a
barotropic equation of state of the type $P=(\gamma_{ad}-1)\rho$, where
$\gamma_{ad}$ is the adiabatic index. The speed of sound $C_{s}$ is defined as
$C^{2}_{s}:=\frac{\partial P}{\partial\rho}$, where the derivative is computed
at constant entropy. However, for the barotropic equation of state chosen
here, the entropy equation decouples from the field equations and therefore
plays no role in our analysis. The speed of sound is then a constant
$\sqrt{\gamma_{ad}-1}$. Due to causality, one must have $0\leq C^{2}_{s}\leq
1$ yielding $\gamma_{ad}\in[1,2]$. This is equivalent to the fact that the
sound cone is contained within the light cone in the tangent space. Here
$\gamma_{ad}=1$ corresponds to pressure-less fluid or ‘dust’ and
$\gamma_{ad}=2$ corresponds to a ultra-relativistic stiff fluid. We note an
important fact that choosing this equation of state, we ignore the rest energy
of the fluid which may be relevant in the late epoch of the universe
evolution. In other words, the equation of state $P=(\gamma_{ad}-1)(\rho-ne)$,
$n$ being the baryon number density and $e$ the fluid rest energy per
particle, may be more appropriate [54]. Here we will stick with the equation
of state $P=(\gamma_{ad}-1)\rho$ for simplicity and in such case, the baryon
number conservation equation $\nabla_{\mu}(nu^{\mu})=0$ is a simple
consequence of the Euler’s equations i.e., $\nabla_{\nu}T^{\mu\nu}=0$. With
the barotropic equation of state, the components of the stress energy tensor
are computed to be

$\displaystyle E=(P+\rho)(1+g(v,v))-P=\rho+\gamma_{ad}\rho g(v,v)$ (23)
$\displaystyle\mathbf{J}^{i}=\gamma_{ad}\rho\sqrt{1+g(v,v)}v^{i},\mathbf{S}=\gamma_{ad}\rho
v\otimes v+(\gamma_{ad}-1)\rho g.$

The complete system of evolution and constraint equations reads

$\displaystyle\partial_{t}\rho+\frac{\gamma\rho
N}{[1+g(v,v)]^{1/2}}\nabla_{i}v^{i}=L_{X}\rho-\frac{\gamma_{ad}\rho
g(\partial_{t}v,v)}{1+g(v,v)}+\frac{\gamma_{ad}\rho
Nk(v,v)}{1+g(v,v)}-\frac{2\gamma_{ad}\rho v_{i}L_{X}v^{i}}{1+g(v,v)}$
$\displaystyle-\frac{\gamma_{ad}N\rho
L_{v}N}{[1+g(v,v)]^{1/2}}+Nk^{i}_{i}\rho-\frac{NL_{v}\rho}{[1+g(v,v)]^{1/2}}.$

and

$\displaystyle\gamma_{ad}\rho\left[u^{0}\partial_{t}v^{i}-2Nu^{0}k^{i}_{j}v^{j}-u^{0}(L_{X}v)^{i}+v^{j}\nabla_{j}v^{i}+(1+g(v,v))\nabla^{i}N\right]$
(24)
$\displaystyle+(\gamma_{ad}-1)\left[\nabla^{i}\rho+v^{i}L_{v}\rho+u^{0}v^{i}(\partial_{t}\rho-
L_{X}\rho)\right]=0$ $\displaystyle\partial_{t}g_{ij}$ $\displaystyle=$
$\displaystyle-2Nk_{ij}+L_{X}g_{ij},$ (25) $\displaystyle\partial_{t}k_{ij}$
$\displaystyle=$ $\displaystyle-\nabla_{i}\nabla_{j}N+N\left\\{R_{ij}+\tau
k_{ij}-2k_{ik}k^{k}_{j}-\frac{2\Lambda}{n-1}g_{ij}\right.$
$\displaystyle\left.-\gamma_{ad}\rho v_{i}v_{j}+\frac{\gamma_{ad}-2}{n-1}\rho
g_{ij}\right\\}+L_{X}k_{ij},$ $\displaystyle R(g)-|k|^{2}+\tau^{2}$
$\displaystyle=$ $\displaystyle
2\Lambda+2\rho\left\\{1+\gamma_{ad}g(v,v)\right\\},$ (27)
$\displaystyle\nabla_{i}k^{ij}-\nabla^{j}\tau$ $\displaystyle=$
$\displaystyle-\gamma_{ad}\rho\sqrt{1+g(v,v)}v^{i},$ (28)

where $u^{0}$ is explicitly computed to be $u^{0}=\frac{[1+g(v,v)]^{1/2}}{N}$.

As clearly seen from the field equations, the gauge variables $N$ and $X$ do
not satisfy evolution equations as physically expected. Therefore, we need to
fix the gauge choice (or choose a gauge slice) to obtain equations for $N$ and
$X$. As we have explicitly mentioned before, we will use the constant mean
extrinsic curvature spatial harmonic gauge (CMCSH). Setting the mean extrinsic
curvature ($\tau:=\tr_{g}k$) of the hypersurface $\Sigma$ as a constant over
$\Sigma$ allows it to play the role of time. However, one is obligated to
address the issue related to the existence of a CMC hypersurface in the
spacetime $\hat{M}$ [55, 56, 57, 58]. Luckily, our background solution (1)
possesses constant mean extrinsic curvature slices, which is evident through
explicit calculations. Therefore simply setting the mean extrinsic curvature
of the perturbed spacetimes to be constant on each spatial slice seems to be a
reasonable choice

$\displaystyle\partial_{i}\tau$ $\displaystyle=$ $\displaystyle 0.$ (29)

Such an assumption can be made without loss of generality for solutions that
are close to the background solution [59, 60]. If one evaluates the time
derivative of $\tau$ using the evolution equations, one obtains the following
equation

$\displaystyle\frac{\partial\tau}{\partial
t}=\Delta_{\gamma}N+\left\\{|k|^{2}_{g}+[\frac{(n\gamma-2)}{n-1}+\gamma_{ad}g(v,v)]\rho-\frac{2\Lambda}{n-1}\right\\}N+L_{X}\tau$
(30)

which upon utilizing the CMC gauge condition (29) yields the elliptic equation
for the lapse function

$\displaystyle\Delta_{g}N+\left\\{|k|^{2}_{g}+[\frac{(n\gamma-2)}{n-1}+\gamma_{ad}g(v,v)]\rho-\frac{2\Lambda}{n-1}\right\\}N=$
$\displaystyle\frac{\partial\tau}{\partial t}.$ (31)

Note that we may obtain an elliptic equation for $\partial_{t}N$ by commuting
the lapse equation (31) with the vector field $\partial_{t}$. Since this
causes equivalent loss of one derivative of every field present in the
equation, care must be taken in order to estimate $\partial_{t}N$. We will do
so in a later section.

In addition to the temporal gauge fixing, we need to fix the spatial gauge as
well. We use the spatial harmonic gauge introduced by Andersson and Moncrief
[61] to study the local well-posedness of the vacuum Einstein equations. Let
$\chi:(\Sigma,g)\to(\Sigma,\gamma)$ be a harmonic map with Dirichlet energy
functional $\frac{1}{2}\int_{\Sigma}g^{ij}\frac{\partial\chi^{k}}{\partial
x^{i}}\frac{\partial\chi^{l}}{\partial x^{j}}\gamma_{jl}\mu_{g}$. Since the
harmonic maps are, by definition, critical points of this Dirichlet energy
functional, $\chi$ satisfies the following formal Euler-Lagrange equations

$\displaystyle
g^{ij}\left(\partial_{i}\partial_{j}\chi^{k}-\Gamma[g]_{ij}^{l}\partial_{l}\chi^{k}+\Gamma[\gamma]_{\alpha\beta}^{k}\partial_{i}\chi^{\alpha}\partial_{j}\chi^{\beta}\right)=0.$
(32)

Now the spatial harmonic gauge fixing amounts to setting $\chi=id$, which
yields a vanishing tension field

$\displaystyle-g^{ij}\left(\Gamma[g]^{k}_{ij}-\Gamma[\gamma]^{k}_{ij}\right)=0.$
(33)

A subtle issue that arises in the context of spatial gauge fixing is the
existence of a harmonic map between $(\Sigma,g)$ and $(\Sigma,\gamma)$. There
are a plethora of studies that address such questions [62]. The one study that
is most relevant to the current context is the one by Eells and Sampson [63]
who proved if $M$ is a compact manifold with non-positive Riemann curvature,
then any continuous map from a compact manifold into $M$ is homotopic to a
harmonic map. Andersson and Moncrief [61] applied this gauge to prove a local
existence theorem for the vacuum Einstein’s equations in arbitrary dimensions.
Later they have successfully implemented this gauge to prove a small data
nonlinear stability result of the ‘Milne’ model and its higher-dimensional
generalizations (so-called ‘Lorentz cone spacetimes’) [27, 28]. In addition to
the study by [27], there have been several other studies where the spatial
harmonic gauge is utilized to study the long-time existence issues of the
Einstein field equations coupled to matter sources. The implementation of the
spatial harmonic gauge condition (33) yields the following elliptic equation
for the shift vector field

$\displaystyle\Delta_{g}X^{i}-R^{i}_{j}X^{j}$ $\displaystyle=$
$\displaystyle-2k^{ij}\nabla_{j}N-2N\nabla_{j}k^{ij}+\tau\nabla^{i}N$
$\displaystyle+(2Nk^{jk}-2\nabla^{j}X^{k})(\Gamma[g]^{i}_{jk}-\Gamma[\gamma]^{i}_{jk})-g^{ij}\partial_{t}\Gamma[\gamma]_{ij}^{k}.$

We designate this gauge fixed system as the CMCSH-Einstein-Euler-$\Lambda$
system.

###  2.2 Re-scaled equations and local well-posedness

The background solution we are interested in is the constant negative
sectional spatial curvature FLRW model in ‘$3+1$’ dimensions. In higher
dimensions ($n>3$), however, the spatial metrics are negative Einstein (not
necessarily hyperbolic). The solution is expressible in the following warped
product form

${}^{n+1}g$ $\displaystyle=$ $\displaystyle-dt\otimes
dt+a(t)^{2}\gamma_{ij}dx^{i}\otimes dx^{j},$ (35)

where $R_{ij}(\gamma)=-\frac{1}{n}\gamma_{ij}$ and $t\in(0,\infty)$. Clearly
we have $N=1,~{}X^{i}=0,~{}g_{ij}=a(t)^{2}\gamma_{ij}$. These spacetimes are
globally foliated by constant mean curvature slices which is evident from the
following simple calculations

$\displaystyle k_{ij}$ $\displaystyle=$
$\displaystyle-\frac{1}{2N}\partial_{t}g_{ij}=-a\dot{a}\gamma_{ij},$ (36)
$\displaystyle\tau(t)$ $\displaystyle=$ $\displaystyle
k_{ij}g^{ij}=-a\dot{a}\gamma_{ij}\frac{1}{a(t)^{2}}\gamma^{ij}=-\frac{n\dot{a}}{a}.$

A standard linear algebra decomposition
$k_{ij}=k^{tr}_{ij}+\frac{\tr_{g}k}{n}g_{ij}$, where $k^{tr}_{ij}$ is the
trace-less part of $k_{ij}$ and $\tr_{g}k=\tau$ is its trace, yields
$k^{tr}=0$ for the background solution (35). Note that in order to study the
dynamics of the perturbations, we need to re-scale the field equations. We
want to re-scale in such a way as to make the background time independent and
then analyze the associated re-scaled dynamical system. For this particular
purpose, we will need the background density and $n-$velocity. Let us now
explicitly calculate the scale factor $a(t)$, the background density, and the
velocity. Note that we are only interested in the asymptotic behaviour of the
scale factor at $t\to\infty$ instead of its precise expression for our
stability analysis. The evolution equation for energy density of the
background solutions reads

$\displaystyle\frac{\partial\rho}{\partial t}$ $\displaystyle=$
$\displaystyle-n\gamma_{ad}\frac{\dot{a}}{a}\rho,$ (37)

integration of which yields

$\displaystyle\rho a(t)^{n\gamma_{ad}}=C_{\rho},$ (38)

for some finite time independent positive constant $C_{\rho}$. The momentum
constraint (28) equation yields

$\displaystyle v^{i}=0$ (39)

since $k^{tr}=0$ and $\partial_{i}\tau=0$ in CMCSH gauge. Utilizing the
Hamiltonian constraint, we obtain the following ODE for the scale factor

$\displaystyle\frac{1}{a}\frac{da}{dt}=\sqrt{\frac{1}{n(n-1)}\left(2\Lambda+\frac{1}{a^{2}}+\frac{2C_{\rho}}{a^{n\gamma_{ad}}}\right)},$
(40)

integration of which yields at $t\to\infty$

$\displaystyle a\sim e^{\alpha t}$ (41)

where $\alpha:=\sqrt{\frac{2\Lambda}{n(n-1)}}$. Taking $\Lambda=0$,
integration of (40) yields

$\displaystyle a\sim a(0)+\frac{t}{\sqrt{n(n-1)}}$ (42)

as $t\to\infty$.

Now let us focus on obtaining a re-scaled system of evolution and constraint
equations. Recall that the metric in local coordinates $(t,x^{i})$ may be
written as

$\displaystyle\hat{g}$ $\displaystyle=$ $\displaystyle-N^{2}dt\otimes
dt+g_{ij}(dx^{i}+X^{i}dt)\otimes(dx^{j}+X^{j}dt)$ (43)

and the background solutions are written as

${}^{n+1}g_{background}=-dt\otimes dt+a^{2}\gamma_{ij}dx^{i}\otimes dx^{j},$
(44)

where $R[\gamma]_{ij}=-\frac{1}{n}\gamma_{ij}$ and $t\in(0,\infty)$. If we
assign $a$ the dimension of length (natural) (similar to the re-scaling used
by [28]), the dimensions of the rest of the entities follow $t\sim
length,~{}x^{i}\sim(length)^{0},~{}g_{ij}\sim(length)^{2},~{}N\sim(length)^{0},~{}X^{i}\sim(length)^{-1},~{}\tau\sim(length)^{-1}$,
$\kappa^{tr}_{ij}\sim length$. Therefore, the following scaling follows
naturally

$\displaystyle\tilde{g}_{ij}=a^{2}g_{ij},\tilde{N}=N,\tilde{X}^{i}=\frac{1}{a}X^{i},\tilde{k^{tr}}_{ij}=ak^{tr}_{ij}.$
(45)

Here, we denote the dimension-full entities with a tilde sign and
dimensionless entities are left alone for simplicity. The Hamiltonian
constraint (27) and the momentum constraint (28) take the following forms

$\displaystyle R(g)-|k^{tr}|^{2}+1+2C_{\rho}a^{2-n\gamma}$ $\displaystyle=$
$\displaystyle
2a(t)^{2}\tilde{\rho}\left(1+\gamma_{ad}\tilde{g}(\tilde{v},\tilde{v})\right),$
(46) $\displaystyle\nabla_{j}k^{trij}$ $\displaystyle=$ $\displaystyle
a(t)^{3}\gamma_{ad}\tilde{\rho}\sqrt{1+\tilde{g}(\tilde{v},\tilde{v})}\tilde{v}^{i}.$

Notice that we are yet to re-scale the matter fields i.e., $(\rho,v)$. We have
already noticed that the background energy density satisfies
$\rho(t)=\frac{C_{\rho}}{a^{n\gamma_{ad}}}$ for some constant $C_{\rho}>0$
depending on the initial density. Therefore, we naturally scale the energy
density by $\frac{1}{a^{n\gamma_{ad}}}$ i.e.,

$\displaystyle\tilde{\rho}$ $\displaystyle=$
$\displaystyle\frac{1}{a(t)^{n\gamma_{ad}}}\rho$ (47)

so that the background density becomes time independent and an $O(1)$ entity.
Notice that $C_{\rho}$ is a constant with dimension of
$(length)^{n\gamma_{ad}-2}$. Now, since, the background $n-$velocity vanishes
for the background solution, we do not have a straightforward scaling.
However, we have the normalization condition for the $n+1-$velocity field
$\textbf{u}=-^{n+1}g(u,n)n+v$

$\displaystyle-(\tilde{N}\tilde{\textbf{u}}^{0})^{2}+\tilde{g}_{ij}\tilde{v}^{i}\tilde{v}^{j}$
$\displaystyle=$ $\displaystyle-1,$ (48)

and therefore scaling of the metric $\tilde{g}_{ij}=a^{2}g_{ij}$ yields the
following natural scaling of $v^{i}$

$\displaystyle\tilde{v}^{i}$ $\displaystyle=$ $\displaystyle\frac{1}{a}v^{i},$
(49)

since the right hand side is simply a constant. The complete evolution and
constraints of the gravity coupled fluid system may be expressed in the
following re-scaled form

$\displaystyle\frac{\partial g_{ij}}{\partial t}$ $\displaystyle=$
$\displaystyle
2\frac{\dot{a}}{a}\left(N-1\right)g_{ij}-\frac{2}{a}Nk^{tr}_{ij}+\frac{1}{a}(L_{X}g)_{ij},$
(50) $\displaystyle\frac{\partial k^{tr}_{ij}}{\partial t}$ $\displaystyle=$
$\displaystyle-\underbrace{\frac{\dot{a}}{a}\left((N-1)(n-2)+n-1\right)}_{should~{}be~{}<0~{}to~{}generate~{}decay}k^{tr}_{ij}-\frac{2}{a}Nk^{tr}_{ik}k^{trk}_{j}+\underbrace{\frac{1}{a}NR_{ij}}_{gauge~{}fixing~{}is~{}required}$
$\displaystyle-\left(\frac{a}{n^{2}}(\tau^{2}-\frac{2n\Lambda}{n-1})+\frac{n\gamma_{ad}-2}{n(n-1)}a^{1-n\gamma_{ad}}C_{\rho}-\frac{N(1+2C_{\rho}a^{2-n\gamma})}{a(n-1)}\right.$
$\displaystyle\left.+\frac{\gamma_{ad}-2}{n-1}a^{1-n\gamma_{ad}}N\rho
g_{ij}\right)g_{ij}-\frac{1}{a}\nabla_{i}\nabla_{j}N-a^{1-n\gamma_{ad}}N\gamma\rho
v_{i}v_{j}+\frac{1}{a}(L_{X}k^{tr})_{ij},$

$\displaystyle\partial_{t}\rho+\frac{\gamma_{ad}\rho
N\nabla_{i}v^{i}}{a[1+g(v,v)]^{1/2}}+\frac{(N-1)n\gamma_{ad}\rho\dot{a}}{a}$
(52) $\displaystyle=\frac{1}{a}L_{X}\rho-\frac{\gamma_{ad}\rho
g(\partial_{t}v,v)}{1+g(v,v)}-\frac{NL_{v}\rho}{a[1+g(v,v)]^{1/2}}+\frac{\gamma_{ad}\rho
N}{1+g(v,v)}\left(\frac{1}{a}k^{tr}(v,v)-\frac{\dot{a}}{a}g(v,v)\right)$
$\displaystyle+\frac{\gamma_{ad}\rho
v_{i}L_{X}v^{i}}{a[1+g(v,v)]}-\frac{\gamma_{ad}N\rho
L_{v}N}{a[1+g(v,v)]^{1/2}},$ $\displaystyle\gamma_{ad}\rho
u^{0}\partial_{t}v^{i}+\underbrace{\gamma_{ad}\rho
u^{0}[2N-1-n(\gamma_{ad}-1)]\frac{\dot{a}}{a}}_{should~{}be~{}>0~{}to~{}generate~{}decay}v^{i}-\frac{2\gamma_{ad}\rho
Nu^{0}k^{tri}_{j}v^{j}}{a}-\frac{\gamma_{ad}\rho u^{0}}{a}L_{X}v^{i}$
$\displaystyle+\frac{\gamma_{ad}\rho}{a}v^{j}\nabla_{j}v^{i}+\frac{\gamma_{ad}\rho[1+g(v,v)]}{a}\nabla^{i}N+\frac{(\gamma_{ad}-1)}{a}\left(\nabla^{i}\rho+v^{i}L_{v}\rho+au^{0}v^{i}\partial_{t}\rho\right.$
(53) $\displaystyle\left.-u^{0}v^{i}L_{X}\rho\right)=0,$ $\displaystyle
R(g)-|k^{tr}|^{2}-2a^{2-n\gamma_{ad}}\rho\\{1+\gamma_{ad}g(v,v)\\}+1+2C_{\rho}a^{2-n\gamma}=0,$
(54)
$\displaystyle\nabla_{j}k^{trij}=-\gamma_{ad}a^{2-n\gamma_{ad}}\rho\sqrt{1+g(v,v)}v^{i},$
(55)
$\displaystyle\Delta_{g}N+\left(\underbrace{|k^{tr}|^{2}+a(t)^{2}(\frac{\tau^{2}}{n}-\frac{2\Lambda}{n-1})+a^{2-n\gamma_{ad}}[\frac{(n\gamma_{ad}-2)}{n-1}+\gamma_{ad}g(v,v)]\rho}_{should~{}be~{}>0~{}for~{}elliptic~{}regularity}\right)N$
$\displaystyle=a^{2}\frac{\partial\tau}{\partial t},$ (56)
$\displaystyle\Delta_{g}X^{i}-R^{i}_{j}X^{j}=-n(1-\frac{2}{n})\dot{a}\nabla^{i}N-2\nabla^{j}Nk^{tri}_{j}-\underbrace{2N\nabla_{j}k^{trij}}_{use~{}momentum~{}constraint}$
(57)
$\displaystyle+(2Nk^{trjk}-2\nabla^{j}X^{k})\left(\Gamma[g]^{i}_{jk}-\Gamma[\gamma]^{i}_{jk}\right)-ag^{ij}\partial_{t}\Gamma[\gamma]_{ij}^{k}.$

In CMCSH gauge, the time function is obtained by setting $\tau=$monotonic
function of $t$ alone. Now, we have $\tau$ evaluated for the background
solution which is a monotonic function of $t$ alone since utilizing the
Hamiltonian constraint and the lapse equation at the background solution one
obtains

$\displaystyle\frac{\partial\tau}{\partial t}$ $\displaystyle=$
$\displaystyle\frac{1}{n}(\tau^{2}-\frac{2n\Lambda}{n-1})+\frac{n\gamma_{ad}-2}{n-1}a^{-n\gamma_{ad}}C_{\rho}$
$\displaystyle=$
$\displaystyle\frac{1}{a(t)^{2}(n-1)}\left(1+n\gamma_{ad}C_{\rho}a^{2-n\gamma_{ad}}\right)>0.$

We set the time function to be the solution of this equation

$\displaystyle\frac{\partial\tau}{\partial
t}=\frac{1}{a(t)^{2}(n-1)}\left(1+n\gamma_{ad}C_{\rho}a(t)^{2-n\gamma_{ad}}\right),$
(59)

i.e., $t=t(\tau)$ which settles the business of choosing a time coordinate.
One subtlety we need to address here is that regarding the structure of the
background solution space for dimensions strictly greater than $3+1$. We do so
in the next section. For now, we will state a local well-posedness theorem
which may be proved using the method developed by [61] or [64]. The proof will
require an energy current for the fluid sector, which will be defined in the
next section. However, since the proof of a local existence theorem follows in
a standard way by using techniques developed by [61], we will not state the
proof here. In addition, one needs to preserve the gauge conditions and
constraints along the solution curve as well. The constraints and gauges read

$\displaystyle
A_{1}:=\tr_{g}k-\int\left(\frac{1}{a(t)^{2}(n-1)}\left(1+n\gamma_{ad}C_{\rho}a(t)^{2-n\gamma_{ad}}\right)\right)dt,$
(60) $\displaystyle
A^{i}_{2}:g^{kl}(\Gamma[g]^{i}_{kl}-\Gamma[\gamma]^{i}_{kl}),$ (61)
$\displaystyle
A^{i}_{3}:=\nabla_{j}k^{trij}+\gamma_{ad}a^{2-n\gamma_{ad}}\rho\sqrt{1+g(v,v)}v^{i},$
(62) $\displaystyle
A_{4}:=R(g)-|k^{tr}|^{2}-2a^{2-n\gamma_{ad}}\rho\\{1+\gamma_{ad}g(v,v)\\}+1+2C_{\rho}a^{2-n\gamma}.$
(63)

Explicit computations utilizing the evolution equations leads to a hyperbolic
system for $(A_{1},A^{i}_{2},A^{i}_{3},A_{4})$ in CMCSH gauge. Using an energy
similar to the one presented in section $4$ of [61], one shows that if
$(A_{1},A^{i}_{2},A^{i}_{3},A_{4})=\textbf{0}$ initially then
$(A_{1},A^{i}_{2},A^{i}_{3},A_{4})\equiv\textbf{0}$ along the solution curve
$t\mapsto(g(t),k^{tr}(t),\rho(t),v(t),N(t),X(t))$. The local well-posedness
theorem is stated as follows
Local Well-posedness Theorem: Let $(g_{0},k^{tr}_{0},\rho_{0},v_{0})\in
H^{s}\times H^{s-1}\times H^{s-1}\times H^{s-1},~{}s>\frac{n}{2}+2$ be the
initial data for the Cauchy problem of the re-scaled Einstein-Euler-$\Lambda$
evolution equations (50-2.2) in Constant Mean Extrinsic Curvature Spatial
Harmonic (CMCSH) gauge satisfying the constraint equations (54-55). This CMCSH
Cauchy problem is well posed in $C\left([0,t^{*}];H^{s}\times H^{s-1}\times
H^{s-1}\times H^{s-1}\right)$. In particular, there exists a time $t^{*}>0$
dependent on
$||g_{0}||_{H^{s}},||k^{tr}_{0}||_{H^{s-1}},||\rho_{0}||_{H^{s-1}},||v_{0}||_{H^{s-1}}$
such that the solution map
$(g_{0},k^{tr}_{0},\rho_{0},v_{0})\mapsto(g(t),k^{tr}(t),\rho(t),v(t),N(t),X(t))$
is continuous as a map

$\displaystyle H^{s}\times H^{s-1}\times H^{s-1}\times H^{s-1}$ (64)
$\displaystyle\to H^{s}\times H^{s-1}\times H^{s-1}\times H^{s-1}\times
H^{s+1}\times H^{s+1}.$

Let $t^{*}$ be the maximal time of existence of a solution to the CMCSH Cauchy
problem with data $(g_{0},k^{tr}_{0},\rho_{0},v_{0})$, then either
$t^{*}=\infty$ or

$\displaystyle\lim_{t\to
t^{*}}\sup\max\left(C(\Omega_{t}),||\nabla[\gamma]g(t)||_{L^{\infty}},||k^{tr}||_{L^{\infty}},||\nabla\rho(t)||_{L^{\infty}},||\rho(t)||_{L^{\infty}},\right.$
(65)
$\displaystyle\left.||\nabla[\gamma]v(t)||_{L^{\infty}},||v(t)||_{L^{\infty}}\right)=\infty,$

where $C(\Omega_{t})$ is a constant dependent on the ellipticity constant of
the metric $\hat{g}$ defined to be
$\Omega_{t^{*}}:=\sup_{t\in[0,t^{*}]}\Omega[\hat{g}(t)]$,
$\Omega[\hat{g}]:=\Omega[g]+||N||_{L^{\infty}}+||N^{-1}||_{L^{\infty}}+||X||_{L^{\infty}}$,
and $\Omega^{-1}\gamma(\xi,\xi)\leq
g(\xi,\xi)\leq\Omega\gamma(\xi,\xi),~{}\forall\xi\in sections\\{T\Sigma\\}$.
In addition the gauges and constraints are preserved along the solution curve
$t\mapsto(g(t),k^{tr}(t),\rho(t),v(t),N(t),X(t))$ if they are satisfied
initially.

###  2.3 Center Manifold of the dynamics

The following mathematics concerning the center manifold dynamics in the
setting of the Einstein flow was first established by the work of Andersson-
Moncrief [27]. We now present some of this work and adapt it to the Einstein-
Euler-$\Lambda$ setting. A center manifold of a dynamical system is associated
with fixed points of the dynamics. Essentially a center manifold of a fixed
point corresponds to the nearby solutions (in phase space) that do not exhibit
exponential growth or decay. This may be zero-dimensional or a subspace of the
phase space and admits a manifold structure (so the name center ‘manifold’).
For a more precise and mathematically rigorous definition, we refer the reader
to [65, 66]. In a finite-dimensional setting, the existence of a center
manifold roughly corresponds to the case when the linearization of the flow
vector field has purely imaginary spectra. In an infinite-dimensional setting,
moduli spaces play the role of center manifolds (in an appropriate sense of
course). As we shall see, the center manifold in this particular occasion is
played by the Einstein moduli space, which consists of non-isolated fixed
points of the Einstein-Euler-$\Lambda$ flow and is characterized by the
negative Einstein metrics modulo gauge transformations.

Unlike the $n=3$ dimensional case where the negative Einstein spaces are
hyperbolic (and the center manifold of the gravitational dynamics consists of
a point), the higher dimensional case is more subtle. The matter degrees of
freedom corresponding to the background solution (in its re-scaled version) in
CMCSH gauge is essentially described by $v^{i}=0,\rho=\rho_{B.G}$. It becomes
more interesting when we consider the gravitational degrees of freedom. In the
case of $n=3$, a negative Einstein structure implies hyperbolic structure
through the Mostow rigidity theorem and therefore the manifold describing the
fixed point of the dynamics is essentially zero-dimensional. In other words,
the particular solution
$(\\{\gamma|Ricci[\gamma]=-\frac{1}{3}\gamma,k^{tr}=0,\rho=C_{\rho},v^{i}=0)$
serves as an isolated fixed point of the Einstein-Euler dynamics. However, in
the higher dimensional cases ($n>3$), a new possibility arises of having non-
hyperbolic negative Einstein spaces. When we linearize about any member of
such a non-isolated family of negative Einstein metrics, the linearized
equations will always admit a finite-dimensional space of neutral modes.
Naturally, this represents (modulo the matter degrees of freedom) the tangent
space to the background spacetimes (1). These smooth families of background
spacetimes determined by the corresponding families of negative Einstein
metrics and zero-dimensional matter degrees of freedom
($v^{i}=0,\rho=C_{\rho}$) form the ‘center manifold’ for the dynamical system
defined by the re-scaled Einstein-Euler-$\Lambda$ equations. The spatial
metric satisfies

$\displaystyle R_{ij}(\gamma)=-\frac{1}{n}\gamma_{ij}.$ (66)

Let us denote the space of metrics satisfying equation (66) by
$\mathcal{E}in_{-\frac{1}{n}}$.

Let $\gamma^{*}\in\mathcal{E}in_{-\frac{1}{n}}$ and $\mathcal{V}$ be its
connected component. Also consider $\mathcal{S}_{\gamma}$ to be the harmonic
slice of the identity diffeomorphism i.e., the set of
$\gamma\in\mathcal{E}in_{-\frac{1}{n}}$ for which the identity map
$id:(\Sigma,\gamma)\to(\Sigma,\gamma^{*})$ is harmonic (since any metric
$\gamma\in\mathcal{E}in_{-\frac{1}{n}}$ verifies the fixed point criteria
$R[\gamma]_{ij}=-\frac{1}{n}\gamma_{ij}$, it should also satisfy the CMCSH
gauge condition with respect to a background metric and in this case the
background metric is simply chosen to be $\gamma^{*}$, another element of
$\mathcal{E}in_{-\frac{1}{n}}$ and $||\gamma-\gamma^{*}||<\delta,\delta>0$
sufficiently small). This condition is equivalent to the vanishing of the
tension field $-V^{k}$ that is

$\displaystyle-V^{k}=-\gamma^{ij}(\Gamma[\gamma]^{k}_{ij}-\Gamma[\gamma^{*}]^{k}_{ij})=0.$
(67)

For $\gamma\in\mathcal{E}in_{-\frac{1}{n}}$, $\mathcal{S}_{\gamma}$ is a
submanifold of $\mathcal{M}_{\Sigma}$ for $\gamma$ sufficiently close to
$\gamma^{*}$ (easily proven using a standard procedure and therefore we omit
the proof, see [28] for details). The deformation space $\mathcal{N}$ of
$\gamma^{*}\in\mathcal{E}in_{-\frac{1}{n}}$ is defined as the intersection of
the $\gamma^{*}-$connected component
$\mathcal{V}\subset\mathcal{E}in_{-\frac{1}{n}}$ and the harmonic slice
$\mathcal{S}_{\gamma}$ i.e.,

$\displaystyle\mathcal{N}:=\mathcal{V}\cap\mathcal{S}_{\gamma}.$ (68)

$\mathcal{N}$ is assumed to be smooth (i.e., equipped with $C^{\infty}$
topology). In the case of $n=3$, following the Mostow rigidity theorem, the
negative Einstein structure is rigid [4]. This rigid Einstein structure in
$n=3$ corresponds to the hyperbolic structure up to isometry. For higher genus
Riemann surfaces $\Sigma_{genus}$ ($genus>1$), the space of distinct
hyperbolic structures is the classical Teichmüller space diffeomorphic to
$\mathbf{R}^{6genus-6}$. The tangent space to any point in the Teichmüller
space corresponds to the ‘neutral modes’ of $2+1$ gravity defined by
transverse-traceless tensors. However, $2+1$ gravity is fundamentally
different from the higher dimensional cases since the former is devoid of
gravitational waves degrees of freedom. For $n>3$, the centre manifold
$\mathcal{N}$ is a finite-dimensional submanifold of $\mathcal{M}_{\Sigma}$
(and in particular of $\mathcal{M}_{-1}$). Following the analysis of [28], the
formal tangent space $T_{\gamma}\mathcal{N}$ in local coordinates is
expressible as

$\displaystyle\frac{\partial\gamma}{\partial q^{a}}$ $\displaystyle=$
$\displaystyle l^{TT||}_{a}+L_{W^{||}_{a}}\gamma,$ (69)

where
$l^{TT||}_{a}\in\ker(\mathcal{L}_{\gamma,\gamma})=C^{TT||}(S^{2}\Sigma)\subset
C^{TT}(S^{2}\Sigma)$, $W^{||}\in\mathfrak{X}(\Sigma)$ satisfies through the
time derivative of the CMCSH gauge condition
$-\gamma^{ij}(\Gamma[\gamma]^{k}_{ij}-\Gamma[\gamma^{*}]^{k}_{ij})=0$

$\displaystyle-[\nabla[\gamma]^{m}\nabla[\gamma]_{m}W^{||i}+R[\gamma]^{i}_{m}W^{||m}]+(l^{TT||}+L_{W^{||}}\gamma)^{mn}(\Gamma[\gamma]^{k}_{mn}-\Gamma[\gamma^{*}]^{k}_{mn})$
$\displaystyle=0,$ (70)

and, $\\{q^{a}\\}_{a=1}^{dim(\mathcal{N})}$ is a local chart on $\mathcal{N}$,
$\mathfrak{X}(\Sigma)$ is the space of vector fields on $\Sigma$ (in a
suitable function space setting), $C^{TT}(S^{2}\Sigma)$ is the space of
$\gamma-$fransverse-traceless $2-$tensors on $\Sigma$. $\gamma^{*}$ is close
to $\gamma$ i.e., $||\gamma-\gamma^{*}||<\delta,\delta>0$ sifficiently small
(the norm is arbitrary since both $\gamma$ and $\gamma^{*}$ are elements of a
finite dimensional vector space $\mathcal{N}$). An important thing to note is
that all known examples of closed negative Einstein spaces have integrable
deformation spaces and that the deformation spaces are stable i.e.,
$Spec\\{\mathcal{L}_{\gamma,\gamma}\\}\geq 0$. The spectrum of the operator
$\mathcal{L}_{\gamma,\gamma}$ plays an important role determining the decay
rates in the pure vacuum case considered by [27, 28]. We stress that the
equation (2.3) plays a crucial role in the energy argument.

###  2.4 Kinematics of the perturbations and shadow gauge

As discussed in the previous section, the Einstein moduli space is finite-
dimensional for $n>3$ and serves as the center manifold of the dynamics (the
deformation space $\mathcal{N}$ to be precise). On the other hand, for the
$n=3$ case, following Mostow rigidity, negative Einstein structure implies
hyperbolicity and therefore the moduli space reduces to a point. We will
consider the general case when the moduli space is finite-dimensional since
the zero-dimensional moduli space case does not require special treatment. Let
$\mathcal{N}$ be the deformation space of $\gamma^{*}$. Now suppose, we
perturb the metric about $\gamma^{*}$. The perturbation is defined to be
$u:=g-\gamma^{*}$ such that $||u||_{H^{s}}<\delta$ for a sufficiently small
$\delta>0$. Due to the finite dimensionality of $\mathcal{N}$, $g$ in fact
could be an element of $\mathcal{N}$ lying within a $\delta-$ball of
$\gamma^{*}$ i.e., while perturbing one may choose directions tangent to
$\mathcal{N}$. But if $g$ lies on $\mathcal{N}$, then corresponding solution
satisfies the fixed point condition $Ric(g)=-\frac{1}{n}g$. Therefore these
‘trivial’ perturbations do not see the dynamics. The question is how to treat
such perturbations which are tangent to $\mathcal{N}$. This was precisely
handled by invoking a shadow gauge condition in the study of the fully
nonlinear stability problem of the vacuum Einstein’ equations by [28] and
later used in [46] to handle the case with a positive cosmological constant.
The shadow gauge reads

$\displaystyle(g-\gamma)\perp\mathcal{N}$ (71)

for $\gamma\in\mathcal{N}$, where the orthogonality is in the $L^{2}$ sense
with respect to the metric $\gamma$. This precisely takes care of the motion
of the orthogonal perturbations by studying the time evolution of
$h:=g-\gamma$, while the motion of the tangential perturbations is taken care
of by evolving $\gamma\in\mathcal{N}$. This does not affect the
characteristics of the background spacetimes as $\gamma$ does not leave the
moduli space i.e., it is still a negative Einstein metric with Einstein
constant $-\frac{1}{n}$. However, if we write the perturbation as
$u_{ij}=g_{ij}-\gamma_{ij}$, then $\frac{\partial g_{ij}}{\partial
t}=\frac{\partial u_{ij}}{\partial t}+\frac{\partial\gamma_{ij}}{\partial t}$
and therefore, $\frac{\partial\gamma_{ij}}{\partial t}$ enters into the
evolution equation for the metric. Therefore this additional contribution
arising as a result of finite dimensionality of the moduli space $\mathcal{N}$
needs to be estimated properly since this will inevitably show up in the time
derivative of the energy functional as well.

Figure 1: Kinematics of the perturbations for the Einstein-Euler system.
Instability materializes in $\mathcal{N}$ iff $\gamma(t)$ leaves every compact
set of $\mathcal{N}$ i.e., approaches $\partial\mathcal{N}$. Stability occurs
when $\gamma(t)$ evolves from one point of $\mathcal{N}$ to simply another
point of $\mathcal{N}$. The instability of $g(t)$ is naturally defined by the
unbounded growth of $g(t)-\gamma(t)$ in a suitable function space.

The shadow gauge condition is explicitly written as follows

$\displaystyle\int_{\Sigma}(g_{ij}-\gamma_{ij})\frac{\partial\gamma^{ij}}{\partial
q^{a}}\mu_{\gamma}=0,$ (72)

where $\frac{\partial\gamma_{ij}}{\partial q^{a}}$ is a local basis for the
tangent space of $\mathcal{N}$ in local coordinates
$\\{q^{a}\\}_{a=1}^{\dim(\mathcal{N})}$ and
$\frac{\partial\gamma^{ij}}{\partial
q^{a}}=-\gamma^{ik}\gamma^{jl}\frac{\partial\gamma_{kl}}{\partial q^{a}}$.
Technically, this shadow gauge condition is equivalent to projecting $g$ onto
$\gamma\in\mathcal{N}$. In another words, there is a projection map
$\mathcal{P}:\mathcal{M}\to\mathcal{N}$ such that

$\displaystyle\gamma=\mathcal{P}[g].$ (73)

Since the deformation space $\mathcal{N}$ (finite dimensional) is assumed to
have a $C^{\infty}$ structure, this projection map $\mathcal{P}$ can also be
interpreted as a smoothing map. Through explicit computations under the
assumption of $||g-\gamma||_{H^{s}}<\delta$, where $\delta>0$ is sufficiently
small, one may explicitly show that the following estimates hold

$\displaystyle||\frac{\partial\gamma}{\partial t}||\leq C||\frac{\partial
g}{\partial t}||_{H^{s-1}},$ (74)

where $C$ is a fixed constant depending only on the dimension and the norm on
the left hand side may be arbitrary since $\mathcal{N}$ is finite dimensional
yielding an equivalence between Sobolev norms of all orders.

##  3 Commuted equations for perturbations and Energy current

In this section, we obtain evolution equations for the perturbations and
construct the commuted equations for the fluid variables. We explicitly
construct an energy current motivated by the energy current that was
constructed by Christodoulou [37] in the context of special relativistic
perfect fluid dynamics. This current has a remarkable property that while
contracted with a time-like vector field it yields positive definite energy.
Since the Euler’s equations are theoretically derivable from a Lagrangian (see
[67] for a symmetric hyperbolic formulation of Euler’s equations and [47, 48,
49] for a Hamiltonian formulation of a relativistic perfect fluid coupled to
gravity), one naturally expects a compatible current (see Christodoulou [68])
associated with Euler equations. Due to the reasons mentioned in the
introduction, we will not use the energy current constructed by [37] but
rather define a current by hand which satisfies the necessary properties. Let
us define the following perturbation entities

$\displaystyle
u:=g-\gamma,~{}k^{tr}:=k^{tr}-0,~{}\omega:=N-1,~{}X:=X-0,~{}\delta\rho:=\rho-
C_{\rho},$ (75) $\displaystyle~{}v:=v-0.$

The evolution equations for the perturbations to the gravitational degrees of
freedom are as follows

$\displaystyle\frac{\partial u_{ij}}{\partial t}$ $\displaystyle=$
$\displaystyle 2\frac{\dot{a}}{a}\omega
g_{ij}-\frac{2}{a}Nk^{tr}_{ij}+h^{TT||}_{ij}+\frac{1}{a}(L_{X}u)_{ij}+(L_{Y^{||}+\frac{1}{a}X}\gamma)_{ij},$
(76) $\displaystyle\frac{\partial k^{tr}_{ij}}{\partial t}$ $\displaystyle=$
$\displaystyle-\frac{\dot{a}}{a}\left(\omega(n-2)+n-1\right)k^{tr}_{ij}-\frac{2}{a}Nk^{tr}_{ik}k^{trk}_{j}$
$\displaystyle+\left(a^{1-n\gamma_{ad}}\frac{\omega\gamma_{ad}C_{\rho}+N(\gamma_{ad}-2)\delta\rho}{n-1}+\frac{\omega}{a(n-1)}\right)g_{ij}$
$\displaystyle-\frac{1}{a}\nabla_{i}\nabla_{j}\omega-a^{1-n\gamma_{ad}}N\gamma\rho
v_{i}v_{j}+\frac{1}{a}(L_{X}k^{tr})_{ij}+\frac{1}{a}N\left(\frac{1}{2}\mathcal{L}_{g,\gamma}u_{ij}+\mathfrak{J}_{ij}\right),$

where we have used the fact that
$\partial_{t}\gamma_{ij}=-h^{TT||}_{ij}-L_{Y^{||}}\gamma_{ij}$ (69)
($h^{TT||}$ is parallel to $\mathcal{N}$ i.e., belongs to the kernel of
$\mathcal{L}_{\gamma,\gamma}$ and $Y^{||}$ satisfies the equation (2.3) with
$l^{TT||}$ replaced by $h^{TT||}$) and $\mathfrak{J}$ contains first
derivatives of $u$, is quasi-linear and satisfies
$||\mathfrak{J}||_{H^{s-1}}\leq C||u||_{H^{s}},~{}s>\frac{n}{2}+2$ which
directly follows from its expression

$\displaystyle\mathfrak{J}_{ij}=\frac{1}{2}\left(u_{il}u^{mn}R[\gamma]^{l}_{mjn}+u_{jl}u^{mn}R[\gamma]^{l}_{min}\right)+\frac{1}{2}g^{mn}\gamma^{ls}\left(\nabla[\gamma]_{j}u_{ns}\nabla[\gamma]_{l}u_{im}\right.$
$\displaystyle\left.+\nabla[\gamma]_{i}u_{lm}\nabla[\gamma]_{s}u_{jn}-\frac{1}{2}\nabla[\gamma]_{j}u_{ns}\nabla[\gamma]_{i}u_{lm}+\nabla[\gamma]_{m}u_{il}\nabla[\gamma]_{n}u_{js}\right.$
$\displaystyle\left.-\nabla[\gamma]_{m}u_{il}\nabla[\gamma]_{s}u_{jn}\right),$

where $u^{ij}:=\gamma^{ik}\gamma^{jl}u_{kl}$ and SH gauge condition is
implemented. From the structure of the above equations (76-3), it is clear
that these are hyperbolic evolution equations since the map
$u\mapsto\mathcal{L}_{g,\gamma}u$ is elliptic. Since the scenario is slightly
different for the case of Euler’s equations, we will first define the commuted
entities and obtain their evolution equations.

Multiple applications of the operator
$\mathcal{D}:=\langle\nabla[\gamma]\rangle=(1+\Delta_{\gamma})^{1/2}$ on the
fields would produce the desired higher order energy since
$W^{s,p}(\Sigma)=\mathcal{D}^{-s}L^{p}(\Sigma),~{}1<p<\infty,~{}s\in\mathbb{R}$,
($\mathcal{D}^{s}$ acts on the spectral side) with norm defined as

$\displaystyle||v||_{W^{s,p}(\Sigma)}=||\mathcal{D}^{s}v||_{L^{p}(\Sigma)}.$
(78)

Here $\Delta_{\gamma}$ is defined as
$\Delta_{\gamma}:=-\gamma^{ij}\nabla[\gamma]_{i}\nabla[\gamma]_{j}$. If $s$ is
a non-negative integer, then $W^{s,p},~{}1\leq p\leq\infty$ is the closure of
$C^{\infty}(\Sigma)$ with respect to the equivalent norm $\sum_{|I|\leq
s}||\nabla[\gamma]^{I}v||_{L^{p}(\Sigma)}$. Let us define the following
entities corresponding to the fluid variables

$\displaystyle\dot{\rho}:=\mathcal{D}^{s-1}(\rho-
C_{\rho}),~{}\dot{v}^{i}:=\mathcal{D}^{s-1}v^{i},~{}~{}s>\frac{n}{2}+2.$ (79)

The evolution equations for the commuted entities read

$\displaystyle\partial_{t}\dot{\rho}+\frac{\gamma_{ad}\rho
N\nabla_{i}\dot{v}^{i}}{a[1+g(v,v)]^{1/2}}=\frac{1}{a}L_{X}\dot{\rho}-\frac{\gamma_{ad}\rho
g(\partial_{t}\dot{v},v)}{1+g(v,v)}-\frac{NL_{v}\dot{\rho}}{a[1+g(v,v)]^{1/2}}$
(80) $\displaystyle+\frac{\gamma_{ad}\rho
v_{i}L_{X}\dot{v}^{i}}{a[1+g(v,v)]}+\mathcal{T}_{1}$

$\displaystyle\gamma_{ad}\rho
u^{0}\partial_{t}\dot{v}^{i}-\frac{\gamma_{ad}\rho
u^{0}}{a}L_{X}\dot{v}^{i}+\frac{\gamma_{ad}\rho}{a}v^{j}\nabla_{j}\dot{v}^{i}+\frac{(\gamma_{ad}-1)}{a}\left(\nabla^{i}\dot{\rho}+v^{i}L_{v}\dot{\rho}+au^{0}v^{i}\partial_{t}\dot{\rho}\right.$
$\displaystyle\left.-u^{0}v^{i}L_{X}\dot{\rho}\right)=\mathcal{T}^{i}_{2}$

where $\mathcal{T}_{1}$ and $\mathcal{T}_{2}$ are defined as follows

$\displaystyle\mathcal{T}_{1}:=-[\mathcal{D}^{s-1},\partial_{t}]\rho-[\mathcal{D}^{s-1},\frac{\gamma_{ad}\rho
N\nabla_{i}}{a[1+g(v,v)]^{1/2}}]v^{i}+\frac{1}{a}[\mathcal{D}^{s-1},L_{X}]\rho$
(81) $\displaystyle-[\mathcal{D}^{s-1},\frac{\gamma_{ad}\rho
g_{ij}v^{j}\partial_{t}}{1+g(v,v)}]v^{i}-\frac{1}{a}[\mathcal{D}^{s-1},\frac{NL_{v}}{[1+g(v,v)]^{1/2}}]\rho+\gamma_{ad}[\mathcal{D}^{s-1},\frac{\rho
v_{i}L_{X}}{a[1+g(v,v)]}]v^{i}$
$\displaystyle+\mathcal{D}^{s-1}\left[\frac{\gamma_{ad}\rho
N}{1+g(v,v)}\left(\frac{1}{a}k^{tr}(v,v)-\frac{\dot{a}}{a}g(v,v)\right)\right]-\mathcal{D}^{s-1}\left[\frac{\gamma_{ad}N\rho
L_{v}N}{a[1+g(v,v)]^{1/2}}\right]$
$\displaystyle-\mathcal{D}^{s-1}\left[\frac{(N-1)n\gamma_{ad}\dot{a}}{a}\rho\right],$
$\displaystyle\mathcal{T}^{i}_{2}:=-\underbrace{\gamma_{ad}\rho
u^{0}[2N-1-n(\gamma_{ad}-1)]\frac{\dot{a}}{a}\dot{v}^{i}}_{one~{}of~{}the~{}most~{}important~{}term:~{}generates~{}decay}$
(82) $\displaystyle+\gamma_{ad}[\mathcal{D}^{s-1},\rho
u^{0}\partial_{t}]v^{i}-\gamma_{ad}\left[\mathcal{D}^{s-1},\rho
u^{0}[2N-1-n(\gamma_{ad}-1)]\right]v^{i}+\frac{2\gamma_{ad}}{a}\mathcal{D}^{s-1}[\rho
Nu^{0}k^{tri}_{j}v^{j}]$
$\displaystyle+\frac{\gamma_{ad}}{a}[\mathcal{D}^{s-1},\rho
u^{0}L_{X}]v^{i}-\frac{\gamma_{ad}}{a}[\mathcal{D}^{s-1},\rho
v^{j}\nabla_{j}]v^{i}+\frac{\gamma_{ad}}{a}\mathcal{D}^{s-1}\left[\frac{\rho[1+g(v,v)]}{N}\nabla^{i}N\right]$
$\displaystyle+\frac{(\gamma_{ad}-1)}{a}\left\\{[\mathcal{D}^{s-1},\nabla^{i}]\rho+[\mathcal{D}^{s-1},v^{i}L_{v}]\rho+a[\mathcal{D}^{s-1},u^{0}v^{i}\partial_{t}]\rho-[\mathcal{D}^{s-1},u^{0}v^{i}L_{X}]\rho\right\\}.$

Notice that we have only retained the principal and decay generating terms in
the evolution equations for commuted entities and collected remaining terms as
$\mathcal{T}_{1}$ and $\mathcal{T}_{2}$. These terms need to be estimated
carefully, which we will do in a later section. Now we define the energy
current. The components of the energy current $\mathcal{C}$ associated with
the commuted fluid variables are defined as follows

$\displaystyle\mathcal{C}^{0}:=\frac{(\gamma_{ad}-1)}{\gamma_{ad}\rho
N}[1+g(v,v)]^{1/2}\dot{\rho}^{2}+2(\gamma_{ad}-1)\frac{g_{ij}v^{i}\dot{v}^{j}}{N[1+g(v,v)]^{1/2}}\dot{\rho}$
$\displaystyle+\frac{\gamma_{ad}\rho[1+g(v,v)]^{1/2}}{N}\left(-\frac{g(v,\dot{v})^{2}}{1+g(v,v)}+g_{ij}\dot{v}^{i}\dot{v}^{j}\right)$

and

$\displaystyle\mathcal{C}^{i}:=\frac{(\gamma_{ad}-1)v^{i}\dot{\rho}^{2}}{\gamma_{ad}\rho}+2(\gamma_{ad}-1)\dot{v}^{i}\dot{\rho}+\gamma_{ad}\rho
v^{i}(-\frac{g(v,\dot{v})^{2}}{1+g(v,v)}+g_{ij}\dot{v}^{i}\dot{v}^{j}).$ (83)

The energy current vector field is
$\mathcal{C}:=C^{0}\partial_{t}+\frac{1}{a}\mathcal{C}^{i}\partial_{i}$. Now
we show how this energy current may be used to obtain an energy. The energy is
naturally defined by contracting the current with the unit time-like Cauchy
hypersurface orthogonal vector field
$\mathbf{n}:=\frac{1}{N}\partial_{t}-\frac{1}{aN}X$

$\displaystyle
E_{fluid}:=-\int_{\Sigma}\mathcal{C}^{\mu}\mathbf{n}_{\mu}\mu_{g},$ (84)

positivity, coercivity, and convexity of which will be shown in lemma (2). We
need to define an energy associated with the gravitational degrees of freedom
as well. Motivated by the energy associated with the linear stability analysis
[23], we define a natural wave equation type of energy (which can be read off
from the evolution equations) for the gravitational degrees of freedom
introduced by [28]

$\displaystyle E_{Ein}=\sum_{I=1}^{s}\left(\frac{1}{2}\langle
k^{tr}|\mathcal{L}^{I-1}_{g,\gamma}k^{tr}\rangle_{L^{2}}+\frac{1}{8}\langle
u|\mathcal{L}^{I}_{g,\gamma}u\rangle_{L^{2}}\right).$ (85)

One may explicitly write down the lowest order gravitational energy as follows

$\displaystyle E^{1}_{Ein}$ (86) $\displaystyle\
=\frac{1}{8}\int_{M}(\nabla[\gamma]_{k}u_{ij}\nabla[\gamma]_{l}u_{mn}g^{kl}\gamma^{im}\gamma^{jn}+(R[\gamma]^{k}~{}_{ij}~{}^{m}+R[\gamma]^{k}~{}_{ji}~{}^{m})u_{km}u_{ln}\gamma^{il}\gamma^{jn})\mu_{g}$
$\displaystyle+\frac{1}{2}\int_{M}k^{tr}_{ij}k^{tr}_{kl}\gamma^{ik}\gamma^{jl}\mu_{g}.$

We want to establish the coercive property of these energies. In other words
we would like to show that these energies control the desired norm of the data
$(u,k^{tr},\dot{\rho},\dot{v})$. Let us define the total energy as follows

$\displaystyle E_{total}:=E_{fluid}+E_{Ein}$ (87)
$\displaystyle=\int_{\Sigma}\left(\frac{(\gamma_{ad}-1)}{\gamma_{ad}\rho}[1+g(v,v)]^{1/2}\dot{\rho}^{2}+2(\gamma_{ad}-1)\frac{g_{ij}v^{i}\dot{v}^{j}}{[1+g(v,v)]^{1/2}}\dot{\rho}\right.$
$\displaystyle\left.+\gamma_{ad}\rho[1+g(v,v)]^{1/2}\left(-\frac{g(v,\dot{v})^{2}}{1+g(v,v)}+g_{ij}\dot{v}^{i}\dot{v}^{j}\right)\right)\mu_{g}$
$\displaystyle+\sum_{i=1}^{s}\left(\frac{1}{2}\langle
k^{tr}|\mathcal{L}^{i-1}_{g,\gamma}k^{tr}\rangle_{L^{2}}+\frac{1}{8}\langle
u|\mathcal{L}^{i}_{g,\gamma}u\rangle_{L^{2}}\right).$

The following two lemmas establish the coercive property of the energy for
small data. Let us write $u=g-\gamma$ as follows

$\displaystyle u=g-\gamma=u^{||}+u^{\perp}+z,$ (88)

where $u^{||}\in\ker(\mathcal{L}_{\gamma,\gamma})$,
$u^{\perp}\perp\ker(\mathcal{L}_{\gamma,\gamma})$ and
$||z||_{H^{s}}\lesssim||u^{||}||^{2}_{H^{s}}+||u^{\perp}||^{2}_{H^{s}}$. In
particular, $u^{||},u^{\perp}\in C^{TT}(S^{2}\Sigma)$.
Lemma 1: Let ($\gamma,g,k^{tr}$) be such that $(g,\gamma)$ satisfy the shadow
gauge and $g-\gamma=u$. There exists $\delta>0$ sufficiently small, and a
constant $C=C(\delta)>0$ such that if
$(u,k^{tr})\in\mathcal{B}_{\delta}(\mathbf{0})\in H^{s}\times
H^{s-1}$,$s>\frac{n}{2}+2$ the following estimate holds

$\displaystyle||u^{||}||_{H^{s}}\leq
C\left(||u^{\perp}||^{2}_{H^{s}}+||k^{tr}||^{2}_{H^{s-1}}\right).$ (89)

Proof: Following the shadow gauge ($u$ is $L^{2}$-orthogonal to $\mathcal{N}$)
(72), we may write

$\displaystyle\langle(g-\gamma),h^{TT||}+L_{Y^{||}}\gamma\rangle_{L^{2}}$
$\displaystyle=$ $\displaystyle 0$ (90) $\displaystyle\langle
u^{||}+u^{\perp}+z,h^{TT||}+L_{Y^{||}}\gamma\rangle_{L^{2}}$ $\displaystyle=$
$\displaystyle 0$ $\displaystyle\langle u^{||},h^{TT||}\rangle_{L^{2}}+\langle
z,h^{TT||}+L_{Y^{||}}\gamma\rangle_{L^{2}}$ $\displaystyle=$ $\displaystyle 0$

where we have used the facts that $\langle
u^{\perp},h^{TT||}\rangle_{L^{2}}=0$ and $\langle
u^{||}+u^{\perp},L_{Y^{||}}\gamma\rangle_{L^{2}}=0$. Using the relation
obtained, we may say that $u^{||}$ is a smooth function of $z$ which satisfies

$\displaystyle||z||_{H^{s}}\leq
C(||u^{||}||^{2}_{H^{s}}+||u^{\perp}||^{2}_{H^{s}}),$ (91)

and therefore, $u^{||}$ satisfies the following estimate

$\displaystyle||u^{||}||_{H^{s}}\leq
C\left(||u^{\perp}||^{2}_{H^{s}}+||k^{tr}||^{2}_{H^{s-1}}\right).~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\Box$
(92)

Lemma 2: Let $s>\frac{n}{2}+2$, $\gamma\in\mathcal{N}$, and $E_{total}$ be the
the total energy defined in (87) and $\gamma_{ad}\in(1,\frac{n+1}{n})$. Then
$\exists~{}\delta>0,C=C(\delta)>0$, such that
$\forall(u,k^{tr},\delta\rho,v)\in\mathcal{B}_{\delta}(\mathbf{0})\in
H^{s}\times H^{s-1}\times H^{s-1}\times H^{s-1}$, the following estimate holds

$\displaystyle||u^{\perp}||^{2}_{H^{s}}+||v||^{2}_{H^{s-1}}+||\delta\rho||^{2}_{H^{s-1}}+||v||^{2}_{H^{s-1}}\leq
CE_{total}.$ (93)

Proof: Firstly notice that the following holds by definition

$\displaystyle||\delta\rho||^{2}_{H^{s-1}}=||\dot{\rho}||^{2}_{L^{2}},~{}||v||^{2}_{H^{s-1}}=||\dot{v}||^{2}_{L^{2}}.$
(94)

First we show that the total energy is positive definite. Assuming
$spec\\{\mathcal{L}_{\gamma,\gamma}\\}>0$, it follows that the gravitational
energy is non-negative definite i.e.,

$\displaystyle E_{Ein}=\sum_{i=1}^{s}\left(\frac{1}{2}\langle
k^{tr}|\mathcal{L}^{i-1}_{g,\gamma}k^{tr}\rangle_{L^{2}}+\frac{1}{8}\langle
u|\mathcal{L}^{i}_{g,\gamma}u\rangle_{L^{2}}\right)\geq 0$ (95)

for $u\notin\ker(\mathcal{L}_{\gamma,\gamma})$ and under the small data
assumption. The equality holds only for $u\equiv 0,~{}k^{tr}\equiv 0$. The
non-negativity of the fluid energy follows due to the following calculations.
Cauchy-Schwartz yields $g(v,\dot{v})^{2}\leq g(v,v)g(\dot{v},\dot{v})$ and
therefore

$\displaystyle-\frac{g(v,\dot{v})^{2}}{1+g(v,v)}+g_{ij}\dot{v}^{i}\dot{v}^{j}\geq
g(\dot{v},\dot{v})-\frac{g(v,v)g(\dot{v},\dot{v})}{1+g(v,v)}=\frac{g(\dot{v},\dot{v})}{1+g(v,v)}$
(96)

and

$\displaystyle-\frac{g(v,\dot{v})^{2}}{1+g(v,v)}+g_{ij}\dot{v}^{i}\dot{v}^{j}\leq
g(\dot{v},\dot{v})+\frac{g(v,v)g(\dot{v},\dot{v})}{1+g(v,v)}=\frac{g(\dot{v},\dot{v})(1+2g(v,v))}{1+g(v,v)}.$

Now consider the following

$\displaystyle
2\frac{g_{ij}v^{i}\dot{v}^{j}}{[1+g(v,v)]^{1/2}}\dot{\rho}\leq\frac{\dot{\rho}^{2}[1+g(v,v)]^{1/2}}{\alpha(n)\gamma_{ad}\rho}+\frac{\alpha(n)\gamma_{ad}\rho
g(v,\dot{v})^{2}}{(1+g(v,v))^{3/2}}$ (97)
$\displaystyle\leq\frac{\dot{\rho}^{2}[1+g(v,v)]^{1/2}}{\alpha(n)\gamma_{ad}\rho}+\frac{\alpha(n)\gamma_{ad}\rho
g(v,v)g(\dot{v},\dot{v})}{(1+g(v,v))^{3/2}}$

which yields

$\displaystyle E_{fluid}\geq$ (98)
$\displaystyle\int_{\Sigma}\left(\frac{\gamma_{ad}-1}{\gamma_{ad}\rho}(1-\frac{1}{\alpha(n)})[1+g(v,v)]^{1/2}\dot{\rho}^{2}+\frac{\gamma_{ad}\rho
g(\dot{v},\dot{v})}{[1+g(v,v)]^{1/2}}\right.$
$\displaystyle\left.-\frac{\alpha(n)\gamma_{ad}(\gamma_{ad}-1)\rho
g(v,v)g(\dot{v},\dot{v})}{(1+g(v,v))^{3/2}}\right)\mu_{g}$
$\displaystyle=\int_{\Sigma}\left(\frac{\gamma_{ad}-1}{\gamma_{ad}\rho}(1-\frac{1}{\alpha(n)})[1+g(v,v)]^{1/2}\dot{\rho}^{2}\right.$
$\displaystyle\left.+\frac{\gamma_{ad}\rho
g(\dot{v},\dot{v})}{(1+g(v,v))^{3/2}}\left(1+g(v,v)-\alpha(n)(\gamma_{ad}-1)g(v,v)\right)\right)\mu_{g}\geq
0$

provided $\alpha(n)(\gamma_{ad}-1)\leq 1$ and $\alpha(n)>1$. The assumption
$\gamma_{ad}<\frac{n+1}{n}$ yields a solution $\alpha(n)$ depending only on
the dimension $n$.

We now compute the Hessian of $E_{total}$ at
$(u,k^{tr},\dot{\rho},\dot{v})=(0,0,0,0)$ for variations $(p,q,r,y)$ and
obtain

$\displaystyle
D^{2}E_{total}|_{(0,0,0,0)}\cdot((p,q,r,y),(p,q,r,y))=\sum_{i=1}^{s}\left(\frac{1}{4}\langle
p,\mathcal{L}^{i}_{\gamma,\gamma}p\rangle_{L^{2}}+\langle
q,\mathcal{L}^{i-1}_{\gamma,\gamma}q\rangle_{L^{2}}\right)$
$\displaystyle+\int_{\Sigma}(\frac{\gamma_{ad}-1}{\gamma_{ad}\rho}r^{2}+\gamma_{ad}C_{\rho}\gamma(y,y))\mu_{\gamma}\geq
0$

and equality holds if and only if $(q,r,y)=(0,0,0)$ and
$p\in\ker(\mathcal{L}_{\gamma,\gamma})$. Therefore,
$D^{2}E_{total}|_{(0,0,0,0)}:(H^{s}-\ker(\mathcal{L}_{\gamma,\gamma}))\times
H^{s-1}\times L^{2}\times L^{2}\to Image(D^{2}E_{total}|_{(0,0,0,0)})$ is an
isomorphism yielding

$\displaystyle||u^{\perp}||^{2}_{H^{s}}+||k^{tr}||^{2}_{H^{s-1}}+||\dot{\rho}||^{2}_{L^{2}}+||\dot{v}||^{2}_{L^{2}}$
(99) $\displaystyle\leq
CD^{2}E_{total}|_{(0,0,0,0)}\cdot((p,q,r,y),(p,q,r,y))$

which upon utilizing the Hilbert space version of the Morse lemma (or Taylor’s
theorem) yields

$\displaystyle||u^{\perp}||^{2}_{H^{s}}+||v||^{2}_{H^{s-1}}+||\delta\rho||^{2}_{H^{s-1}}+||v||^{2}_{H^{s-1}}\leq
CE_{total}.~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\Box$
(100)

Now using the previous lemma (lemma 1), we obtain the following coercive
property of the total energy

$\displaystyle||u||^{2}_{H^{s}}+||v||^{2}_{H^{s-1}}+||\delta\rho||^{2}_{H^{s-1}}+||v||^{2}_{H^{s-1}}\leq
CE_{total}.$ (101)

This is one of the most important properties that is necessary for the global
existence proof. Now we will obtain a differential equation for the total
energy. The time derivative of the total energy reads

$\displaystyle\frac{dE_{total}}{dt}=\frac{dE_{Ein}}{dt}+\frac{dE_{fluid}}{dt}.$
(102)

Here the gravitational part $\frac{dE_{Ein}}{dt}$ can be evaluated the same
way as in [28] and only the matter coupling terms need to be kept track of.
The delicate part is the computation of the time derivative of the fluid
energy. Therefore we focus on evaluating $\frac{dE_{fluid}}{dt}$ first. The
following lemma yields a differential equation for the fluid energy.

Lemma 3: Let the adiabatic index $\gamma_{ad}$ lie in the interval
$(1,\frac{n+1}{n})$ and $a(t)$, $t\in(0,\infty)$ be the scale factor. Assume
$(u,k^{tr},\dot{\rho},\dot{v})$ solve the Einstein’s and commuted Euler’s
equations (50-2.2). Then the fluid energy
$E_{fluid}:=-\int_{\Sigma}\mathcal{C}^{\mu}\mathbf{n}_{\mu}\mu_{g}$ satisfies
the following differential equation in time

$\displaystyle\frac{dE_{fluid}}{dt}=\frac{n\dot{a}}{a}\int_{\Sigma}N(N-1)\mathcal{C}^{0}\mu_{g}+\int_{\Sigma}\left(\partial_{t}\mathcal{C}^{0}+\frac{1}{a}\nabla_{i}\mathcal{C}^{i}+(\frac{1}{a}\nabla_{i}X^{i}+\frac{1}{N}\partial_{t}N)\mathcal{C}^{0}\right.$
$\displaystyle\left.+\frac{1}{aN}\mathcal{C}^{i}\nabla_{i}N\right)N\mu_{g}.$

Proof: The proof is a consequence of an elementary divergence identity for the
energy current $\mathcal{C}$. Integrating the divergence identity (on the
scale-free geometry)

$\displaystyle\nabla_{\mu}\mathcal{C}^{\mu}=\frac{1}{N\mu_{g}}\partial_{\mu}(N\mu_{g}\mathcal{C}^{\mu})$
(103)

over a Cauchy hypersurface $\Sigma$ yields

$\displaystyle\int_{\Sigma}\nabla_{\mu}\mathcal{C}^{\mu}N\mu_{g}=\int_{\Sigma}\partial_{\mu}(N\mu_{g}\mathcal{C}^{\mu})=\int_{\Sigma}\partial_{t}(N\mu_{g}\mathcal{C}^{0})+\frac{1}{a}\int_{\Sigma}\partial_{i}(N\mu_{g}\mathcal{C}^{i})$
(104)
$\displaystyle=\partial_{t}\int_{\Sigma_{t}}N\mathcal{C}^{0}\mu_{g}=-\frac{d}{dt}\int_{\Sigma_{t}}\mathcal{C}^{\mu}\mathbf{n}_{\mu}\mu_{g}=\frac{dE_{fluid}}{dt},$

where we have used Stokes’s theorem on closed manifold to remove the total
derivative term. Explicit computation of $\nabla_{\mu}\mathcal{C}^{\mu}$
yields the following differential equation for the fluid energy $E_{fluid}$

$\displaystyle\frac{dE_{fluid}}{dt}=\frac{n\dot{a}}{a}\int_{\Sigma}N(N-1)\mathcal{C}^{0}\mu_{g}+\int_{\Sigma}\left(\partial_{t}\mathcal{C}^{0}+\frac{1}{a}\nabla_{i}\mathcal{C}^{i}+(\frac{1}{a}\nabla_{i}X^{i}+\frac{1}{N}\partial_{t}N)\mathcal{C}^{0}\right.$
$\displaystyle\left.+\frac{1}{aN}\mathcal{C}^{i}\nabla_{i}N\right)N\mu_{g}.~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\Box$

Notice that this energy equation contains the top order term
$\partial_{t}\mathcal{C}^{0}+\frac{1}{a}\nabla_{i}\mathcal{C}^{i}$ which from
the definition of the energy current, will consist of terms that are
derivatives of the commuted variables $\dot{\rho}$ and $\dot{v}^{i}$. This is
dangerous since the energy does not contain any differentiated $\dot{\rho}$
and $\dot{v}^{i}$ terms. However, by using explicit computations, we will show
that such dangerous terms cancel point-wise with their negative counter-parts
or yield terms that are of total derivative forms and therefore vanish by
integration. This cancellation property exhibits the hyperbolic nature of
Euler’s equations.

Claim: The entity
$\partial_{t}\mathcal{C}^{0}+\frac{1}{a}\nabla_{i}\mathcal{C}^{i}$ contains
purely algebraic terms of $\dot{\rho}$ and $\dot{v}$ modulo terms involving
Lie derivatives with respect to the shift vector field $X$.

Proof: We compute the term
$\partial_{t}\mathcal{C}^{0}+\frac{1}{a}\nabla_{i}\mathcal{C}^{i}$ explicitly.
First we compute $\partial_{t}\mathcal{C}^{0}$ and $\nabla_{i}\mathcal{C}^{i}$
separately. We write the following

$\displaystyle\partial_{t}\mathcal{C}^{0}=[\partial_{t}\mathcal{C}^{0}]_{principal}+[\partial_{t}\mathcal{C}^{0}]_{remainder},$
(105)
$\displaystyle\nabla_{i}\mathcal{C}^{i}=[\nabla_{i}\mathcal{C}^{i}]_{principal}+[\nabla_{i}\mathcal{C}^{i}]_{remainder},$
(106)

where the subscript $principal$ denotes the terms that involve highest
derivatives of $\dot{\rho}$ and $\dot{v}$ and $remainder$ denotes the terms
that are purely algebraic in $(\dot{\rho},\dot{v})$. Explicit calculation
yields

$\displaystyle[\partial_{t}\mathcal{C}^{0}]_{principal}=\frac{2(\gamma_{ad}-1)}{\gamma_{ad}\rho
N}[1+g(v,v)]^{1/2}\dot{\rho}\partial_{t}\dot{\rho}+2(\gamma_{ad}-1)\left(\frac{\dot{\rho}g_{ij}v^{i}}{N[1+g(v,v)]^{1/2}}\partial_{t}\dot{v}^{j}\right.$
$\displaystyle\left.+\frac{g_{ij}v^{i}\dot{v}^{j}}{N[1+g(v,v)]^{1/2}}\partial_{t}\dot{\rho}\right)+\frac{\gamma_{ad}\rho[1+g(v,v)]^{1/2}}{N}\left[(2g_{ij}(\partial_{t}\dot{v}^{i})\dot{v}^{j}-\frac{2g(v,\dot{v})g_{ij}v^{i}\partial_{t}\dot{v}^{j}}{1+g(v,v)}\right],$

$\displaystyle[\nabla_{i}\mathcal{C}^{i}]_{principal}=\frac{2(\gamma_{ad}-1)}{\gamma_{ad}\rho}\dot{\rho}v^{i}\nabla_{i}\dot{\rho}+2(\gamma_{ad}-1)\left[\dot{\rho}\nabla_{i}\dot{v}^{i}+\dot{v}^{i}\nabla_{i}\dot{\rho}\right]$
(107) $\displaystyle+\gamma_{ad}\rho
v^{i}\left[2\dot{v}_{k}\nabla_{i}\dot{v}^{k}-\frac{2g(v,\dot{v})g(v,\nabla_{i}\dot{v})}{1+g(v,v)}\right],$

$\displaystyle[\partial_{t}\mathcal{C}^{0}]_{remainder}$ (108)
$\displaystyle=\frac{\gamma_{ad}-1}{\gamma_{ad}}\left(\frac{[1+g(v,v)]^{-1/2}}{2\rho
N}(\partial_{t}g(v,v)+2g(\partial_{t}v,v))-\frac{[1+g(v,v)]^{1/2}}{\rho^{2}N^{2}}[\rho\partial_{t}N+N\partial_{t}\rho]\right)\dot{\rho}^{2}$
$\displaystyle+2(\gamma_{ad}-1)\dot{\rho}\dot{v}^{j}\left(-\frac{\partial_{t}Ng_{ij}v^{i}}{N^{2}[1+g(v,v)]^{1/2}}-\frac{g_{ij}v^{i}}{N[1+g(v,v)]^{3/2}}[\partial_{t}g(v,v)+2g(\partial_{t}v,v)]\right.$
$\displaystyle\left.+\frac{1}{N[1+g(v,v)]^{1/2}}(\partial_{t}g_{ij}v^{i}+g_{ij}\partial_{t}v^{i})\right)+\gamma_{ad}\left[g(\dot{v},\dot{v})-\frac{g(v,\dot{v})^{2}}{1+g(v,v)}\right]\left(\frac{[1+g(v,v)]^{1/2}}{N}\partial_{t}\rho\right.$
$\displaystyle\left.+\frac{\rho[1+g(v,v)]^{-1/2}}{2N}(\partial_{t}g(v,v)+2g(\partial_{t}v,v))-\frac{\rho[1+g(v,v)]^{1/2}}{N^{2}}\partial_{t}N\right)$
$\displaystyle+\frac{\gamma_{ad}\rho[1+g(v,v)]^{1/2}}{N}\left[(\partial_{t}g_{ij})\dot{v}^{i}\dot{v}^{j}+\frac{g(v,\dot{v})^{2}}{[1+g(v,v)]^{2}}(\partial_{t}g(v,v)+2g(\partial_{t}v,v))\right.$
$\displaystyle\left.-\frac{g(v,\dot{v})}{1+g(v,v)}(\partial_{t}g_{ij}v^{i}\dot{v}^{j}+g_{ij}\dot{v}^{i}\partial_{t}v^{j})\right],$
$\displaystyle[\nabla_{i}\mathcal{C}^{i}]_{remainder}$ (109)
$\displaystyle=-\frac{\gamma_{ad}-1}{\gamma_{ad}}\frac{\dot{\rho}^{2}}{\rho^{2}}v^{i}\nabla_{i}\rho+\frac{\gamma_{ad}-1}{\gamma_{ad}}\frac{\dot{\rho}^{2}}{\rho}\nabla_{i}v^{i}+\gamma_{ad}[v^{i}\nabla_{i}\rho+\rho\nabla_{i}v^{i}]$
$\displaystyle\left[-\frac{g(v,\dot{v})^{2}}{1+g(v,v)}+g(\dot{v},\dot{v})\right]+\gamma_{ad}\rho
v^{i}\left(\frac{2g(v,\nabla_{i}v)g(v,\dot{v})^{2}}{[1+g(v,v)]^{2}}-\frac{2g(v,\dot{v})g(\dot{v},\nabla_{i}v)}{1+g(v,v)}\right).$

Now notice that there are several terms involving derivatives of the commuted
fluid variables $(\dot{\rho},\dot{v})$ in the expressions of
$[\partial_{t}\mathcal{C}^{0}]_{principal}$ and
$[\nabla_{i}\mathcal{C}^{i}]_{principal}$ that appear in a complicated
fashion. Therefore we need to carefully combine these terms to obtain
cancellation of the dangerous terms. Using the commuted equation (80), we
expand the first term $\frac{2(\gamma_{ad}-1)}{\gamma_{ad}\rho
N}[1+g(v,v)]^{1/2}\dot{\rho}\partial_{t}\dot{\rho}$ of
$[\partial_{t}\mathcal{C}^{0}]_{principal}$ as follows

$\displaystyle\frac{2(\gamma_{ad}-1)}{\gamma_{ad}\rho
N}[1+g(v,v)]^{1/2}\dot{\rho}\partial_{t}\dot{\rho}$ (110)
$\displaystyle=\frac{2(\gamma_{ad}-1)}{\gamma_{ad}\rho
N}[1+g(v,v)]^{1/2}\dot{\rho}\left(-\frac{\gamma_{ad}\rho
N\nabla_{i}\dot{v}^{i}}{a[1+g(v,v)]^{1/2}}+\frac{1}{a}L_{X}\dot{\rho}-\frac{\gamma_{ad}\rho
g(\partial_{t}\dot{v},v)}{1+g(v,v)}\right.$
$\displaystyle\left.-\frac{NL_{v}\dot{\rho}}{a[1+g(v,v)]^{1/2}}+\frac{\gamma_{ad}\rho
v_{i}L_{X}\dot{v}^{i}}{a[1+g(v,v)]}+\mathcal{T}_{1}\right).$

Now collection of the principal terms together with the previous expression
yields

$\displaystyle\left(\partial_{t}\mathcal{C}^{0}+\frac{1}{a}\nabla_{i}\mathcal{C}^{i}\right)_{principal}$
(111) $\displaystyle=\frac{2(\gamma_{ad}-1)}{\gamma_{ad}\rho
N}[1+g(v,v)]^{1/2}\dot{\rho}\left(-\frac{\gamma_{ad}\rho
N\nabla_{i}\dot{v}^{i}}{a[1+g(v,v)]^{1/2}}+\frac{1}{a}L_{X}\dot{\rho}-\underbrace{\frac{\gamma_{ad}\rho
g(\partial_{t}\dot{v},v)}{1+g(v,v)}}_{B}\right.$
$\displaystyle\left.-\frac{NL_{v}\dot{\rho}}{a[1+g(v,v)]^{1/2}}+\frac{\gamma_{ad}\rho
v_{i}L_{X}\dot{v}^{i}}{a[1+g(v,v)]}+\mathcal{T}_{1}\right)+\underbrace{\frac{2(\gamma_{ad}-1)\dot{\rho}g(v,\partial_{t}\dot{v})}{N[1+g(v,v)]^{1/2}}}_{-B}$
$\displaystyle+\underbrace{\frac{2(\gamma_{ad}-1)g(v,\dot{v})\partial_{t}\dot{\rho}}{N[1+g(v,v)]^{1/2}}+\frac{2\gamma_{ad}\rho[1+g(v,v)]^{1/2}}{N}g(\dot{v},\partial_{t}\dot{v})-\frac{2\gamma_{ad}\rho
g(v,\dot{v})g(v,\partial_{t}\dot{v})}{N[1+g(v,v)]^{1/2}}}_{1A:~{}Combine~{}and~{}use~{}equations~{}of~{}variations}$
$\displaystyle+\frac{1}{a}\left(\frac{2(\gamma_{ad}-1)}{\gamma_{ad}\rho}\dot{\rho}v^{i}\nabla_{i}\dot{\rho}+2(\gamma_{ad}-1)[\dot{\rho}\nabla_{i}\dot{v}^{i}+\dot{v}^{i}\nabla_{i}\dot{\rho}]+\gamma_{ad}\rho
v^{i}[2\dot{v}_{k}\nabla_{i}\dot{v}^{k}\right.$
$\displaystyle\left.-\frac{2g(v,\dot{v})g(v,\nabla_{i}\dot{v})}{1+g(v,v)}]\right).$

We have numbered the terms here so that we can keep track of each dangerous
term and how they cancel with their negative counterpart. First notice that
the terms denoted by $B$ and $-B$ cancel each other (we denote them
accordingly). Now we consider the combination of three terms denoted by $1A$

$\displaystyle
1A=\frac{2(\gamma_{ad}-1)g(v,\dot{v})\partial_{t}\dot{\rho}}{N[1+g(v,v)]^{1/2}}+\frac{2\gamma_{ad}\rho[1+g(v,v)]^{1/2}}{N}g(\dot{v},\partial_{t}\dot{v})-\frac{2\gamma_{ad}\rho
g(v,\dot{v})g(v,\partial_{t}\dot{v})}{N[1+g(v,v)]^{1/2}}$
$\displaystyle=\frac{2g(v,\dot{v})}{N[1+g(v,v)]^{1/2}}\left((\gamma_{ad}-1)\partial_{t}\dot{\rho}-\gamma_{ad}\rho
g(v,\partial_{t}\dot{v})\right)+\frac{2\gamma_{ad}\rho[1+g(v,v)]^{1/2}}{N}g(\dot{v},\partial_{t}\dot{v})$
$\displaystyle=\frac{2g(v,\dot{v})}{N[1+g(v,v)]^{1/2}}\left((\gamma_{ad}-1)\partial_{t}\dot{\rho}-g(v,\frac{\gamma_{ad}\rho}{a}L_{X}\dot{v}-\frac{\gamma_{ad}\rho}{au^{0}}v^{k}\nabla_{k}\dot{v}-\frac{\gamma_{ad}-1}{au^{0}}\nabla\dot{\rho}\right.$
$\displaystyle\left.-\frac{\gamma_{ad}-1}{au^{0}}vL_{v}\dot{\rho}-(\gamma_{ad}-1)v\partial_{t}\dot{\rho}+\frac{\gamma_{ad}-1}{a}vL_{X}\dot{\rho}+\frac{\mathcal{T}_{2}}{u^{0}})\right)$
$\displaystyle+\frac{2[1+g(v,v)]^{1/2}}{N}g\left(\dot{v},\frac{\gamma_{ad}\rho}{a}L_{X}\dot{v}-\frac{\gamma_{ad}\rho}{au^{0}}v^{k}\nabla_{k}\dot{v}-\frac{\gamma_{ad}-1}{au^{0}}\nabla\dot{\rho}-\frac{\gamma_{ad}-1}{au^{0}}vL_{v}\dot{\rho}\right.$
$\displaystyle\left.-(\gamma_{ad}-1)v\partial_{t}\dot{\rho}+\frac{\gamma_{ad}-1}{a}vL_{X}\dot{\rho}+\frac{\mathcal{T}_{2}}{u^{0}}\right)$
$\displaystyle=\underbrace{\frac{2g(v,\dot{v})}{N[1+g(v,v)]^{1/2}}\left((\gamma_{ad}-1)\partial_{t}\dot{\rho}(1+g(v,v))\right)-\frac{2[1+g(v,v)]^{1/2}}{N}(\gamma_{ad}-1)g(v,\dot{v})\partial_{t}\dot{\rho}}_{point-
wise~{}cancellation}$
$\displaystyle+\frac{2g(v,\dot{v})}{N[1+g(v,v)]^{1/2}}\left(-\frac{\gamma_{ad}\rho}{a}g(v,L_{X}\dot{v})+\frac{\gamma_{ad}\rho}{au^{0}}g(v,v^{k}\nabla_{k}\dot{v})+\frac{\gamma_{ad}-1}{au^{0}}L_{v}\dot{\rho}\right.$
$\displaystyle\left.+\frac{\gamma_{ad}-1}{au^{0}}g(v,v)L_{v}\dot{\rho}-\frac{\gamma_{ad}-1}{a}g(v,v)L_{X}\dot{\rho}-g(v,\frac{\mathcal{T}_{2}}{u^{0}})\right)$
$\displaystyle+\frac{2[1+g(v,v)]^{1/2}}{N}\left(\frac{\gamma_{ad}\rho}{a}g(\dot{v},L_{X}\dot{v})-\frac{\gamma_{ad}\rho}{au^{0}}g(\dot{v},v^{k}\nabla_{k}\dot{v})-\frac{\gamma_{ad}-1}{au^{0}}\dot{v}^{k}\nabla_{k}\dot{\rho}\right.$
$\displaystyle\left.-\frac{\gamma_{ad}-1}{au^{0}}g(v,\dot{v})L_{v}\dot{\rho}+\frac{\gamma_{ad}-1}{a}g(v,\dot{v})L_{X}\dot{\rho}+g(\dot{v},\frac{\mathcal{T}_{2}}{u^{0}})\right)+\frac{2(\gamma_{ad}-1)[1+g(v,v)]^{1/2}}{\gamma_{ad}\rho
N}\dot{\rho}\mathcal{T}_{1},$
$\displaystyle=\frac{2g(v,\dot{v})}{N[1+g(v,v)]^{1/2}}\left(-\frac{\gamma_{ad}\rho}{a}g(v,L_{X}\dot{v})+\frac{\gamma_{ad}\rho}{au^{0}}g(v,v^{k}\nabla_{k}\dot{v})+\frac{\gamma_{ad}-1}{au^{0}}L_{v}\dot{\rho}\right.$
$\displaystyle\left.+\frac{\gamma_{ad}-1}{au^{0}}g(v,v)L_{v}\dot{\rho}-\frac{\gamma_{ad}-1}{a}g(v,v)L_{X}\dot{\rho}-g(v,\frac{\mathcal{T}_{2}}{u^{0}})\right)$
$\displaystyle+\frac{2[1+g(v,v)]^{1/2}}{N}\left(\frac{\gamma_{ad}\rho}{a}g(\dot{v},L_{X}\dot{v})-\frac{\gamma_{ad}\rho}{au^{0}}g(\dot{v},v^{k}\nabla_{k}\dot{v})-\frac{\gamma_{ad}-1}{au^{0}}\dot{v}^{k}\nabla_{k}\dot{\rho}\right.$
$\displaystyle\left.-\frac{\gamma_{ad}-1}{au^{0}}g(v,\dot{v})L_{v}\dot{\rho}+\frac{\gamma_{ad}-1}{a}g(v,\dot{v})L_{X}\dot{\rho}+g(\dot{v},\frac{\mathcal{T}_{2}}{u^{0}})\right)+\frac{2(\gamma_{ad}-1)[1+g(v,v)]^{1/2}}{\gamma_{ad}\rho
N}\dot{\rho}\mathcal{T}_{1}.$

Substitution of $1A$ in the expression of
$(\partial_{t}\mathcal{C}^{0}+\frac{1}{a}\nabla_{i}\mathcal{C}^{i})_{principal}$
(after carrying out the possible point-wise cancellations) yields

$\displaystyle(\partial_{t}\mathcal{C}^{0}+\frac{1}{a}\nabla_{i}\mathcal{C}^{i})_{principal}=\underbrace{-\frac{2(\gamma_{ad}-1)}{a}\dot{\rho}\nabla_{i}\dot{v}^{i}}_{I}+\frac{2(\gamma_{ad}-1)(1+g(v,v))^{1/2}}{a\gamma_{ad}\rho
N}\dot{\rho}L_{X}\dot{\rho}$
$\displaystyle-\underbrace{\frac{2(\gamma_{ad}-1)}{a\gamma_{ad}\rho}\dot{\rho}L_{v}\dot{\rho}}_{II}+\underbrace{\frac{2(\gamma_{ad}-1)}{aN[1+g(v,v)]^{1/2}}\dot{\rho}g(v,L_{X}\dot{v})}_{III}$
$\displaystyle+\frac{1}{a}\left(\underbrace{\frac{2(\gamma_{ad}-1)}{\gamma_{ad}\rho}\dot{\rho}v^{i}\nabla_{i}\dot{\rho}}_{-II}+2(\gamma_{ad}-1)[\underbrace{\dot{\rho}\nabla_{i}\dot{v}^{i}}_{-I}+\underbrace{\dot{v}^{i}\nabla_{i}\dot{\rho}}_{IV}]+\gamma_{ad}\rho
v^{i}[\underbrace{2\dot{v}_{k}\nabla_{i}\dot{v}^{k}}_{V}\right.$
$\displaystyle\left.-\underbrace{\frac{2g(v,\dot{v})g(v,\nabla_{i}\dot{v})}{1+g(v,v)}]}_{VI}\right)+\frac{2g(v,\dot{v})}{N[1+g(v,v)]^{1/2}}\left(-\frac{\gamma_{ad}\rho}{a}g(v,L_{X}\dot{v})\right.$
$\displaystyle\left.+\underbrace{\frac{\gamma_{ad}\rho}{au^{0}}g(v,v^{k}\nabla_{k}\dot{v})}_{-VI}+\underbrace{\frac{\gamma_{ad}-1}{au^{0}}L_{v}\dot{\rho}+\frac{\gamma_{ad}-1}{au^{0}}g(v,v)L_{v}\dot{\rho}}_{VII}\right.$
$\displaystyle\left.-\underbrace{\frac{\gamma_{ad}-1}{a}g(v,v)L_{X}\dot{\rho}}_{III}-g(v,\frac{\mathcal{T}_{2}}{u^{0}})\right)$
$\displaystyle+\frac{2[1+g(v,v)]^{1/2}}{N}\left(\frac{\gamma_{ad}\rho}{a}g(\dot{v},L_{X}\dot{v})-\underbrace{\frac{\gamma_{ad}\rho}{au^{0}}g(\dot{v},v^{k}\nabla_{k}\dot{v})}_{-V}-\underbrace{\frac{\gamma_{ad}-1}{au^{0}}\dot{v}^{k}\nabla_{k}\dot{\rho}}_{-IV}\right.$
$\displaystyle\left.-\underbrace{\frac{\gamma_{ad}-1}{au^{0}}g(v,\dot{v})L_{v}\dot{\rho}}_{-VII}+\underbrace{\frac{\gamma_{ad}-1}{a}g(v,\dot{v})L_{X}\dot{\rho}}_{III}+g(\dot{v},\frac{\mathcal{T}_{2}}{u^{0}})\right)$
$\displaystyle+\frac{2(\gamma_{ad}-1)[1+g(v,v)]^{1/2}}{\gamma_{ad}\rho
N}\dot{\rho}\mathcal{T}_{1}.$

Here the terms that are denoted by the same number belong to the same
category, that is, the combination of these terms yields lower-order innocuous
terms. The remaining terms are denoted by numbers and their negative
counterparts are denoted by the negative of their numbers. Note that there are
three terms that fall under the category $III$ that requires the following few
lines of additional calculations to generate lower-order innocuous terms

$\displaystyle\dot{\rho}g(v,L_{X}\dot{v})=\dot{\rho}g(v,X^{k}\nabla_{k}\dot{v}-\dot{v}^{k}\nabla_{k}X)$
(112)
$\displaystyle=\dot{\rho}X^{k}\nabla_{k}g(v,\dot{v})-\dot{\rho}g(X^{k}\nabla_{k}v,\dot{v})-\dot{\rho}g(v,\dot{v}^{k}\nabla_{k}X)$
$\displaystyle=X^{k}\nabla_{k}(\dot{\rho}g(v,\dot{v}))-g(v,\dot{v})L_{X}\dot{\rho}-\dot{\rho}g(X^{k}\nabla_{k}v,\dot{v})-\dot{\rho}g(v,\dot{v}^{k}\nabla_{k}X)$
$\displaystyle=-g(v,\dot{v})L_{X}\dot{\rho}-L_{X}(\dot{\rho}g(v,\dot{v}))-\dot{\rho}g(X^{k}\nabla_{k}v,\dot{v})-\dot{\rho}g(v,\dot{v}^{k}\nabla_{k}X).$

After all the possible cancellations, we collect the remaining terms to obtain
the following final expression for
$\partial_{t}\mathcal{C}^{0}+\frac{1}{a}\nabla_{i}\mathcal{C}^{i}$

$\displaystyle\int_{\Sigma}(\partial_{t}\mathcal{C}^{0}+\frac{1}{a}\nabla_{i}\mathcal{C}^{i})N\mu_{g}$
(113)
$\displaystyle=\int_{\Sigma}\left(\underbrace{\frac{2(\gamma_{ad}-1)(1+g(v,v))^{1/2}}{a\gamma_{ad}\rho
N}\dot{\rho}L_{X}\dot{\rho}}_{I}-\underbrace{\frac{2\gamma_{ad}\rho
g(v,\dot{v})}{aN[1+g(v,v)]^{1/2}}g(v,L_{X}\dot{v})}_{II}\right.$
$\displaystyle\left.-\frac{2g(v,\dot{v})}{N[1+g(v,v)]^{1/2}}g(v,\frac{\mathcal{T}_{2}}{u^{0}})+\underbrace{\frac{2\gamma_{ad}\rho[1+g(v,v)]^{1/2}}{aN}g(\dot{v},L_{X}\dot{v})}_{III}+\frac{2[1+g(v,v)]^{1/2}}{N}g(\dot{v},\frac{\mathcal{T}_{2}}{u^{0}})\right.$
$\displaystyle\left.-\frac{2(\gamma_{ad}-1)}{aN[1+g(v,v)]^{1/2}}\left(\underbrace{L_{X}(\dot{\rho}g(v,\dot{v}))}_{IV}+\dot{\rho}g(X^{k}\nabla_{k}v,\dot{v})+\dot{\rho}g(v,\dot{v}^{k}\nabla_{k}X)\right)\right.$
$\displaystyle\left.+\frac{2(\gamma_{ad}-1)[1+g(v,v)]^{1/2}}{\gamma_{ad}\rho
N}\dot{\rho}\mathcal{T}_{1}\right.$
$\displaystyle\left.+\frac{\gamma_{ad}-1}{\gamma_{ad}}\left(\frac{[1+g(v,v)]^{-1/2}}{2\rho
N}(\partial_{t}g(v,v)+2g(\partial_{t}v,v))-\frac{[1+g(v,v)]^{1/2}}{\rho^{2}N^{2}}[\rho\partial_{t}N+N\partial_{t}\rho]\right)\dot{\rho}^{2}\right.$
$\displaystyle\left.+2(\gamma_{ad}-1)\dot{\rho}\dot{v}^{j}\left(-\frac{\partial_{t}Ng_{ij}v^{i}}{N^{2}[1+g(v,v)]^{1/2}}-\frac{g_{ij}v^{i}}{N[1+g(v,v)]^{3/2}}[\partial_{t}g(v,v)+2g(\partial_{t}v,v)]\right.\right.$
$\displaystyle\left.\left.+\frac{1}{N[1+g(v,v)]^{1/2}}(\partial_{t}g_{ij}v^{i}+g_{ij}\partial_{t}v^{i})\right)+\gamma_{ad}\left[g(\dot{v},\dot{v})-\frac{g(v,\dot{v})^{2}}{1+g(v,v)}\right]\left(\frac{[1+g(v,v)]^{1/2}}{N}\partial_{t}\rho\right.\right.$
$\displaystyle\left.\left.+\frac{\rho[1+g(v,v)]^{-1/2}}{2N}(\partial_{t}g(v,v)+2g(\partial_{t}v,v))-\frac{\rho[1+g(v,v)]^{1/2}}{N^{2}}\partial_{t}N\right)\right.$
$\displaystyle\left.+\frac{\gamma_{ad}\rho[1+g(v,v)]^{1/2}}{N}\left[(\partial_{t}g_{ij})\dot{v}^{i}\dot{v}^{j}+\frac{g(v,\dot{v})^{2}}{[1+g(v,v)]^{2}}(\partial_{t}g(v,v)+2g(\partial_{t}v,v))\right.\right.$
$\displaystyle\left.\left.-\frac{g(v,\dot{v})}{1+g(v,v)}(\partial_{t}g_{ij}v^{i}\dot{v}^{j}+g_{ij}\dot{v}^{i}\partial_{t}v^{j})\right]+\frac{1}{a}\left(-\frac{\gamma_{ad}-1}{\gamma_{ad}}\frac{\dot{\rho}^{2}}{\rho^{2}}v^{i}\nabla_{i}\rho+\frac{\gamma_{ad}-1}{\gamma_{ad}}\frac{\dot{\rho}^{2}}{\rho}\nabla_{i}v^{i}\right.\right.$
$\displaystyle\left.\left.+\gamma_{ad}[v^{i}\nabla_{i}\rho+\rho\nabla_{i}v^{i}]\left[-\frac{g(v,\dot{v})^{2}}{1+g(v,v)}+g(\dot{v},\dot{v})\right]+\gamma_{ad}\rho
v^{i}\left(\frac{2g(v,\nabla_{i}v)g(v,\dot{v})^{2}}{[1+g(v,v)]^{2}}\right.\right.\right.$
$\displaystyle\left.\left.\left.-\frac{2g(v,\dot{v})g(\dot{v},\nabla_{i}v)}{1+g(v,v)}\right)\right)\right)N\mu_{g}.$

The terms that are numbered $I,II,III$, and $IV$ contain derivatives of
$\dot{\rho}$ and $\dot{v}$. However, these are innocuous since these
derivatives involve Lie derivatives with respect to the shift vector field
$X$. We will reduce these terms to lower order terms through an integration by
parts operation which will be discussed in the next section. The remaining
terms are purely algebraic in nature. This completes the proof of the claim.
Recall from equation (102) and lemma 3 that the time derivative of the total
energy reads

$\displaystyle\frac{dE_{total}}{dt}=\frac{dE_{Ein}}{dt}+\int_{\Sigma}\left(\partial_{t}\mathcal{C}^{0}+\frac{1}{a}\nabla_{i}\mathcal{C}^{i}+(\frac{1}{a}\nabla_{i}X^{i}+\frac{1}{N}\partial_{t}N)\mathcal{C}^{0}\right.$
$\displaystyle\left.+\frac{1}{aN}\mathcal{C}^{i}\nabla_{i}N\right)N\mu_{g}+\frac{n\dot{a}}{a}\int_{\Sigma}N(N-1)\mathcal{C}^{0}\mu_{g}.$

Using the reduced expression for
$\partial_{t}\mathcal{C}^{0}+\frac{1}{a}\nabla_{i}\mathcal{C}^{i}$ (113), we
will estimate $\frac{dE_{total}}{dt}$.

##  4 Energy estimates

In this section, we will obtain a differential inequality for the total
energy. For this purpose, we will need the following set of inequalities.

###  4.1 Important inequalities

1. [69] Let $f_{1}\in H^{t_{1}},~{}f_{2}\in H^{t_{2}}$ and $k\leq\min(t_{1},t_{2},t_{1}+t_{2}-\frac{n}{2}),~{}t_{i}\geq 0$ and some $t_{i}>0$, then the following estimate holds

$\displaystyle||f_{1}f_{2}||_{H^{k}}\leq
C||f_{1}||_{H^{t_{1}}}||f_{2}||_{H^{t_{2}}}$ (114)

2. [70] For $\mathcal{P}\in\mathcal{OP}^{s},~{}s\in\mathbb{R}$, the following holds

$\displaystyle||\mathcal{P}f||_{H^{r}}\leq
C||f||_{H^{s+r}},~{}~{}r\in\mathbb{R}.$ (115)

Here $\mathcal{OP}^{s}$ denotes the space of pseudo-differential operators
with symbol in Hormander’s class $S^{s}_{1,0}$ (see [70]). For example
$\nabla[\gamma]\in\mathcal{OP}^{1}$.
3.[70] Assume $\mathcal{P}\in\mathcal{OP}^{s},~{}s>0$, then

$\displaystyle||[P,f_{1}]f_{2}||_{H^{\sigma}}\leq
C\left(||\nabla[\gamma]f_{1}||_{L^{\infty}}||f_{2}||_{H^{s-1+\sigma}}+||f_{1}||_{H^{s+\sigma}}||f_{2}||_{L^{\infty}}\right)$
(116)

4.[70] Let $1\leq s<\mu$ and $F\in C^{\mu}(\mathbb{R})$ with $F(0)=0$, where
$C^{\mu}(\mathbb{R})$ denotes the Hölder space. Then for $f\in H^{s}\cap
L^{\infty}$, the following inequality holds

$\displaystyle||F(f)||_{H^{s}}\leq
C||f||_{H^{s}}(1+||f||^{\mu-1}_{L^{\infty}})$ (117)

5.[71] If $s>0$, then $H^{s}\cap L^{\infty}$ is an algebra i.e.,

$\displaystyle||f_{1}f_{2}||_{H^{s}}\leq
C(||f_{1}||_{L^{\infty}}||f_{2}||_{H^{s}}+||f_{1}||_{H^{s}}||f_{2}||_{L^{\infty}})$
(118)

and for $s>\frac{n}{2}+1$

$\displaystyle||f_{1}f_{2}||_{H^{s}}\leq C||f_{1}||_{H^{s}}||f_{2}||_{H^{s}}$
(119)

due to Sobolev embedding.

In addition to these inequalities, we will also use Hölder’s and Minkowski’s
inequality whenever necessary. First we obtain the necessary elliptic
estimates

###  4.2 Elliptic estimates

In this section, we will obtain estimates for the perturbation to the lapse
function $\omega=N-1$, the shift vector field $X$, and the shifted shift
vector field $Y^{||}+\frac{1}{a}X$. These estimates will be obtained by
utilizing the elliptic equations that are obtained by means of gauge fixing.
Before moving to the analysis of the elliptic equations, we need to obtain an
estimate for $\partial_{t}\gamma$ describing the motion parallel to the center
manifold $\mathcal{N}$. The following lemma provides an estimate for
$\partial_{t}\gamma$

Lemma 4: Let $\gamma\in\mathcal{N}$ be the shadow of g such that
$||g-\gamma||_{H^{s}}<\delta$ for a suitable $\delta>0$, $s>\frac{n}{2}+2$.
The $\mathcal{N}-$tangential velocity $\partial_{t}\gamma$ verifies the
following estimate

$\displaystyle||\partial_{t}\gamma||_{H^{s-1}}$ (120) $\displaystyle\leq
C\left(\frac{\dot{a}}{a}||\omega||_{H^{s+1}}(1+||u||_{H^{s}})+\frac{1}{a}(1+||\omega||_{H^{s+1}})||k^{tr}||_{H^{s-1}}\right.$
$\displaystyle\left.+\frac{1}{a}||X||_{H^{s+1}}||u||_{H^{s}}+\frac{1}{a}||X||_{H^{s+1}}\right),$

for a constant $C$ that depends on a fixed background geometry.
Proof: In the light of the shadow gauge condition (72), the
$\mathcal{N}-$tangential velocity $\partial_{t}\gamma$ satisfies

$\displaystyle||\partial_{t}\gamma||\leq C||\partial_{t}g||_{H^{s-1}},$ (121)

where the norm in the left hand side refers to a Sobolev norm of arbitrary
order (since $\partial_{t}\gamma$ is an element of a finite dimensional vector
space). Utilizing the evolution equation (50) and Minkowski’s inequality
together with the fact that $H^{s_{1}}\hookrightarrow H^{s_{2}}$ for
$s_{1}>s_{2}$, we obtain

$\displaystyle||\partial_{t}\gamma||$ (122) $\displaystyle\leq
C\left(\frac{\dot{a}}{a}||\omega||_{H^{s+1}}(1+||u||_{H^{s}})+\frac{1}{a}(1+||\omega||_{H^{s+1}})||k^{tr}||_{H^{s-1}}\right.$
$\displaystyle\left.+\frac{1}{a}||X||_{H^{s+1}}||u||_{H^{s}}+\frac{1}{a}||X||_{H^{s+1}}\right),$

and in particular

$\displaystyle||\partial_{t}\gamma||_{H^{s-1}}$ (123) $\displaystyle\leq
C\left(\frac{\dot{a}}{a}||\omega||_{H^{s+1}}(1+||u||_{H^{s}})+\frac{1}{a}(1+||\omega||_{H^{s+1}})||k^{tr}||_{H^{s-1}}\right.$
$\displaystyle\left.+\frac{1}{a}||X||_{H^{s+1}}||u||_{H^{s}}+\frac{1}{a}||X||_{H^{s+1}}\right)$

with a possible modification of the constant
$C$.$~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\Box$

Now we will obtain the estimates for $\omega$, $\partial_{t}N$ and $X$
utilizing the elliptic equations obtained through gauge fixing. The following
lemma provides an estimate for the perturbation to the elapse function.

Lemma 5: Let $B_{\delta}(\mathbf{0})\subset H^{s}\times H^{s-1}\times
H^{s-1}\times H^{s-1}$ be a ball for sufficiently small $\delta>0$ and
$(u,k^{tr},\delta\rho,v)\in B_{\delta}(\mathbf{0})$ solve the re-scaled
Einstein-Euler-$\Lambda$ system in CMCSH gauge. Then the following estimate
holds for the perturbation to the lapse function ($\omega=N-1$)

$\displaystyle||\omega||_{H^{s+1}}\leq$ (124) $\displaystyle
C(||k^{tr}||^{2}_{H^{s-1}}+a^{2-n\gamma_{ad}}||\delta\rho||_{H^{s-1}}+a^{2-n\gamma_{ad}}||v||^{2}_{H^{s-1}}(1+||\delta\rho||_{H^{s-1}})).$

Proof: The elliptic equation for the perturbation $\omega=N-1$ to the lapse
function $N$ reads

$\displaystyle\Delta_{g}\omega+\underbrace{\left(|k^{tr}|^{2}+\frac{2C_{\rho}a^{2-n\gamma_{ad}}}{n-1}+\frac{1}{n-1}+a^{2-n\gamma_{ad}}[\frac{(n\gamma_{ad}-2)}{n-1}+\gamma_{ad}g(v,v)](C_{\rho}+\delta\rho)\right)}_{>0}\omega$
$\displaystyle=-|k^{tr}|^{2}_{g}-\frac{n\gamma_{ad}-2}{n-1}a^{2-n\gamma_{ad}}\delta\rho-\gamma_{ad}C_{\rho}a^{2-n\gamma_{ad}}g(v,v)-\gamma_{ad}a^{2-n\gamma_{ad}}g(v,v)\delta\rho.$

Now the elliptic operator
$\Delta_{g}+\left(|k^{tr}|^{2}+\frac{2C_{\rho}a^{2-n\gamma_{ad}}}{n-1}+\frac{1}{n-1}+a^{2-n\gamma_{ad}}[\frac{(n\gamma_{ad}-2)}{n-1}+\gamma_{ad}g(v,v)](C_{\rho}+\delta\rho)\right)id:H^{s+1}\to
H^{s-1}$ is an isomorphism yielding the following estimate

$\displaystyle||\omega||_{H^{s+1}}\leq$ (125) $\displaystyle
C(||k^{tr}||^{2}_{H^{s-1}}+a^{2-n\gamma_{ad}}||\delta\rho||_{H^{s-1}}+a^{2-n\gamma_{ad}}||v||^{2}_{H^{s-1}}(1+||\delta\rho||_{H^{s-1}})).$

A short proof by contradiction may be as follows. Let us denote the elliptic
operator in question by $\mathcal{Q}$. Assume that an estimate of the type
$||\omega||_{H^{s+1}}\lesssim||\mathcal{Q}\omega||_{H^{s-1}}$ does not hold.
Then there exists a sequence $\\{\omega_{i}\\}_{i=1}^{\infty}$ with
$||\omega_{i}||_{H^{s+1}}=1$ and $||\mathcal{Q}\omega_{i}||_{H^{s-1}}\to 0$ as
$i\to\infty$. Now $\Sigma$ is compact and therefore $H^{s+1}$ is compactly
embedded into $H^{s-1}$ (notice $s>\frac{n}{2}+2$). This yields a sub-sequence
$\\{\omega_{i_{j}}\\}$ converging to $\omega^{*}\in H^{s+1}$ strongly in
$H^{s-1}$, which by construction satisfies $\mathcal{Q}\omega^{*}=0$. This
contradicts the fact that the operator $\Omega$ is injective. Therefore an
estimate of the type $||\omega||_{H^{s+1}}\lesssim||A\omega||_{H^{s-1}}$ holds
$\forall s\geq 1$ (which is satisfied here since $s>\frac{n}{2}+2$).

Estimation of the shift vector field is slightly more involved since the CMCSH
gauge is known to exhibit the so called Gribov degeneracy [72, 73]. The shift
equation reads

$\displaystyle\Delta_{g}X^{i}-R^{i}_{j}X^{j}=-n(1-\frac{2}{n})\dot{a}\nabla^{i}N-2\nabla^{j}Nk^{tri}_{j}-2N\nabla_{j}k^{trij}$
(126)
$\displaystyle+(2Nk^{trjk}-2\nabla^{j}X^{k})\left(\Gamma[g]^{i}_{jk}-\Gamma[\gamma]^{i}_{jk}\right)-ag^{ij}\partial_{t}\Gamma[\gamma]_{ij}^{k}.$

In order to estimate $X$ in terms of the remaining variables, we want to show
that the map
$X\mapsto\Delta_{g}X^{i}-R^{i}_{j}X^{j}+2\nabla^{j}X^{k}\left(\Gamma[g]^{i}_{jk}-\Gamma[\gamma]^{i}_{jk}\right)$
is an isomorphism. However, apart from the operator $\Delta_{g}$, the
remaining two terms do not have a definite sign. Therefore, it is not obvious
that such an isomorphism property holds. However, here we will explicitly
prove the isomorphism property of this operator using the small data
assumption ($||g-\gamma||_{H^{s}}<\delta$). The following lemma provides the
necessary estimate for the shift vector field.

Lemma 6: Let $B_{\delta}(\mathbf{0})\subset H^{s}\times H^{s-1}\times
H^{s-1}\times H^{s-1}$ be a ball for sufficiently small $\delta>0$ and
$(u,k^{tr},\delta\rho,v)\in B_{\delta}(\mathbf{0})$ solve the re-scaled
Einstein-Euler-$\Lambda$ system in CMCSH gauge. Then the following estimate
holds for the shift vector field

$\displaystyle||X||_{H^{s+1}}\leq
C\left(\dot{a}||\omega||_{H^{s+1}}+||\omega||_{H^{s+1}}||k^{tr}||_{H^{s-1}}+a^{2-n\gamma_{ad}}\right.$
(127)
$\displaystyle\left.\\{||v||^{3}_{H^{s-1}}(1+||u||_{H^{s}})+||v||_{H^{s-1}}\\}+||k^{tr}||_{H^{s-1}}||u||_{H^{s}}\right)$
$\displaystyle+C\left(\dot{a}||\omega||_{H^{s+1}}(1+||u||_{H^{s}})+||k^{tr}||_{H^{s-1}}\right)$

Proof: The first step in obtaining the estimate for $X$ is to prove that the
map
$X\mapsto\Delta_{g}X^{i}-R^{i}_{j}X^{j}+2\nabla^{j}X^{k}\left(\Gamma[g]^{i}_{jk}-\Gamma[\gamma]^{i}_{jk}\right)$
is an isomorphism between suitable function spaces. Let us denote the flow of
$X$ on a constant time hypersurface by $\varphi_{s}$ ($s$ being the flow
parameter). Let us pull back $V^{\flat}$ (where
$V^{i}:=g^{kl}(\Gamma[g]^{i}_{kl}-\Gamma[\gamma]^{i}_{kl})$ is the negative of
the tension vector field) by the flow of $X$ to obtain the following

$\displaystyle[(\varphi^{*}_{s}V^{\flat})^{\sharp}]^{i}=(\varphi^{*}_{s}g)^{jk}(\varphi^{*}_{s}(\Gamma[g]-\Gamma[\gamma]))^{i}_{kl}=(\varphi^{*}_{s}g)^{jk}((\Gamma[\varphi^{*}_{s}g]-\Gamma[\varphi^{*}_{s}\gamma]))^{i}_{kl}$

which holds because of the fact that the difference of connections transforms
as a section of a suitable vector bundle. Differentiating both sides with
respect to $s$ and setting $s=0$ yields

$\displaystyle
L_{X}V^{i}=(L_{X}g)^{jk}(\Gamma[g]^{i}_{jk}-\Gamma[\gamma]^{i}_{jk})+g^{jk}\frac{d}{ds}\left((\Gamma[\varphi^{*}_{s}g]-\Gamma[\varphi^{*}_{s}\gamma])^{i}_{jk}\right)|_{s=0}$
$\displaystyle=(L_{X}g)^{jk}(\Gamma[g]^{i}_{jk}-\Gamma[\gamma]^{i}_{jk})+\frac{1}{2}g^{jk}g^{il}\left[\nabla_{j}(L_{X}g_{lk})+\nabla_{k}(L_{X}g_{jl})\right.$
$\displaystyle\left.-\nabla_{l}(L_{X}g_{jk})\right]-\frac{1}{2}g^{jk}\gamma^{il}\left[\nabla[\gamma]_{j}(L_{X}\gamma_{lk})+\nabla[\gamma]_{k}(L_{X}\gamma_{jl})-\nabla[\gamma]_{l}(L_{X}\gamma_{jk})\right]$

Explicit computation in CMCSH gauge ($V^{i}=0$) yields

$\displaystyle
g^{kl}\nabla_{k}\nabla_{l}X^{i}+R^{i}_{j}X^{j}+2\nabla^{k}X^{l}(\Gamma[g]^{i}_{kl}-\Gamma[\gamma]^{i}_{kl})$
(128)
$\displaystyle=g^{kl}\nabla[\gamma]_{k}\nabla[\gamma]_{l}X^{i}_{\gamma}+g^{kl}R[\gamma]_{mk}~{}^{i}~{}_{l}X^{m}_{\gamma},$

where $X^{i}_{\gamma}:=X_{j}\gamma^{ij}$. Now the map
$P:X\mapsto-g^{kl}\nabla[\gamma]_{k}\nabla[\gamma]_{l}X^{i}_{\gamma}-g^{kl}R[\gamma]_{mk}~{}^{i}~{}_{l}X^{m}_{\gamma}$
is an isomorphism which may be shown as follows. Let’s consider that $\ker(P)$
is not empty and $0\neq Z\in\ker(P)$ i.e.,

$\displaystyle-g^{kl}\nabla[\gamma]_{k}\nabla[\gamma]_{l}Z^{i}-g^{kl}R[\gamma]_{mk}~{}^{i}~{}_{l}Z^{m}=0.$
(129)

Multiplying both sides of the above equation by $Z$ and integrating over
$\Sigma$ yields

$\displaystyle\int_{\Sigma}\left(-g^{kl}Z_{i}\nabla[\gamma]_{k}\nabla[\gamma]_{l}Z^{i}-g^{kl}Z_{i}R[\gamma]_{mk}~{}^{i}~{}_{l}Z^{m}\right)\mu_{\gamma}=0$
(130)
$\displaystyle\int_{\Sigma}\left(\gamma^{kl}\nabla[\gamma]_{k}Z_{i}\nabla[\gamma]_{l}Z^{i}-\gamma^{kl}R[\gamma]_{mkil}Z^{i}Z^{m}-(g^{kl}-\gamma^{kl})Z_{i}\nabla[\gamma]_{k}\nabla[\gamma]_{l}Z^{i}\right.$
$\displaystyle\left.-(g^{kl}-\gamma^{kl})R[\gamma]_{mkil}Z^{i}Z^{m}\right)\mu_{\gamma}=0$

Now for $||g-\gamma||_{H^{s}}<\delta$ sufficiently small, we have

$\displaystyle\left(\gamma^{kl}\nabla[\gamma]_{k}Z_{i}\nabla[\gamma]_{l}Z^{i}-\gamma^{kl}R[\gamma]_{mkil}Z^{i}Z^{m}-(g^{kl}-\gamma^{kl})Z_{i}\nabla[\gamma]_{k}\nabla[\gamma]_{l}Z^{i}\right.$
(131)
$\displaystyle\left.-(g^{kl}-\gamma^{kl})R[\gamma]_{mkil}Z^{i}Z^{m}\right)>0$

since $Ricci[\gamma](Z,Z)<0$ yielding

$\displaystyle Z=0$ (132)

i.e., $\ker(P)=\\{0\\}$. Since $\Sigma$ is compact and $P$ is second order
elliptic, showing injectivity is enough to ensure its isomorphism property.
Therefore in CMCSH gauge, we establish the isomorphism of the map $X\mapsto
g^{kl}\nabla_{k}\nabla_{l}X^{i}+R^{i}_{j}X^{j}+2\nabla^{k}X^{l}(\Gamma[g]^{i}_{kl}-\Gamma[\gamma]^{i}_{kl})$.
Using the equation for the shift (126), this isomorphism yields the following
estimate for the shift vector field

$\displaystyle||X||_{H^{s+1}}\leq
C||\left((1-\frac{2}{n})a\tau\nabla^{i}N-2\nabla^{j}Nk^{tri}_{j}\right.$
$\displaystyle\left.-2N\left(\gamma_{ad}a^{2-n\gamma_{ad}}\rho(Nu^{0}-1)v^{i}+\gamma_{ad}a^{2-n\gamma_{ad}}\rho
v^{i}\right)+2Nk^{trjk}(\Gamma[g]^{i}_{jk}-\Gamma[\gamma]^{i}_{jk})\right.$
$\displaystyle
r\left.-ag^{ij}\partial_{t}\Gamma[\gamma]_{ij}^{k}\right)||_{H^{s-1}}.$ (133)

The source term that is present in the right-hand side of the previous
estimate (4.2) contains several terms that are multiplied by $\dot{a}$ which
may be potentially dangerous. We need to keep track of these terms carefully.
The most problematic term $ag^{ij}\partial_{t}\Gamma[\gamma]_{ij}^{k}$ needs
to be estimated (notice that it is multiplied by $a$). Since
$\partial_{t}\gamma$ is tangent to $\mathcal{N}$, it may be written as follows

$\displaystyle\partial_{t}\gamma$ $\displaystyle=$
$\displaystyle-h^{TT||}-L_{Y^{||}}\gamma,$ (134)

with $Y^{||}$ satisfying (2.3)

$\displaystyle[\nabla[\gamma]^{m}\nabla[\gamma]_{m}Y^{||i}+R[\gamma]^{||i}_{m}Y^{||m}]-(h^{TT||}+L_{Y^{||}}\gamma)^{mn}(\Gamma[\gamma]^{i}_{mn}$
(135) $\displaystyle-\Gamma[\gamma^{*}]^{i}_{mn})=0$

and $h^{TT||}\in\ker(\mathcal{L}_{\gamma,\gamma})$. Now
$\gamma^{mn}\partial_{t}\Gamma^{i}_{mn}[\gamma]=\gamma^{mn}D\Gamma^{i}_{mn}[\gamma]\cdot\partial_{t}\gamma$
for $\partial_{t}\gamma\in T_{\gamma}\mathcal{N}$ may be written as

$\displaystyle\gamma^{mn}D\Gamma[\gamma]^{i}_{mn}\cdot\partial_{t}\gamma$
$\displaystyle=$
$\displaystyle\gamma^{mn}\gamma^{ik}\left(\nabla[\gamma]_{m}\partial_{t}\gamma_{kn}+\nabla[\gamma]_{n}\partial_{t}\gamma_{mk}-\nabla[\gamma]_{k}\partial_{t}\gamma_{mn}\right)$
$\displaystyle=$ $\displaystyle
2\gamma^{ik}\nabla[\gamma]^{n}\partial_{t}\gamma_{kn}-\nabla[\gamma]^{i}(tr_{\gamma}\partial_{t}\gamma),$

from which $\gamma^{mn}D\Gamma^{i}_{mn}[\gamma]h^{TT||}=0$ follows immediately
since $\nabla[\gamma]^{i}h^{TT||}_{ij}=0,~{}\tr_{\gamma}h^{TT||}=0$. The
remaining terms lead to the following

$\displaystyle-\gamma^{mn}D\Gamma[\gamma]^{i}_{mn}\cdot\partial_{t}\gamma=2\gamma^{ik}\nabla[\gamma]^{n}\left(\nabla[\gamma]_{k}Y^{||}_{n}+\nabla[\gamma]_{n}Y^{||}_{k}\right)-2\nabla[\gamma]^{i}(\nabla[\gamma]_{m}Y^{||m})$
$\displaystyle=2\gamma^{ik}(\nabla[\gamma]_{k}\nabla[\gamma]_{n}Y^{n}+R[\gamma]_{mk}Y^{||m})+2\nabla[\gamma]^{n}\nabla[\gamma]_{n}Y^{||i}-2\nabla[\gamma]^{i}(\nabla[\gamma]_{m}Y^{||m})$
$\displaystyle=2(\nabla[\gamma]^{n}\nabla[\gamma]_{n}Y^{||i}+R[\gamma]^{i}_{m}Y^{||m})$
$\displaystyle=2(L_{Y^{||}}\gamma+h^{TT})^{mn}(\Gamma[\gamma]^{i}_{mn}-\Gamma[\gamma^{*}]^{i}_{mn})$
$\displaystyle=2(\partial_{t}\gamma)^{mn}(\Gamma[\gamma]^{i}_{mn}-\Gamma[\gamma^{*}]^{i}_{mn}).$

This yields

$\displaystyle||\gamma^{mn}D\Gamma[\gamma]^{i}_{mn}\partial_{t}\gamma||_{H^{s-1}}\leq
C||\gamma-\gamma^{*}||_{H^{s}}||\partial_{t}\gamma||_{H^{s-1}}$ (136)

Therefore, using the triangle inequality, we have the following

$\displaystyle||g^{mn}D\Gamma[\gamma]^{i}_{mn}\partial_{t}\gamma||_{H^{s-1}}\leq
C(||g-\gamma||_{H^{s}}+||\gamma-\gamma^{*}||_{H^{s}})||\partial_{t}\gamma||_{H^{s-1}}.$
(137)

Utilizing the estimate of $\partial_{t}\gamma$ from lemma 4, we obtain the
following estimate for the shift vector field $X$

$\displaystyle||X||_{H^{s+1}}\left(1-C(||u||_{H^{s}}+||\gamma-\gamma^{*}||_{H^{s}})(1+||u||_{H^{s}})\right)$
(138) $\displaystyle\leq
C\left(\dot{a}||\omega||_{H^{s+1}}+||\omega||_{H^{s+1}}||k^{tr}||_{H^{s-1}}+a^{2-n\gamma_{ad}}(1+||\delta\rho||_{H^{s-1}})\right.$
$\displaystyle\left.(||v||^{2}_{H^{s-1}}(1+||u||_{H^{s}})||v||_{H^{s-1}}+||v||_{H^{s-1}})+||k^{tr}||_{H^{s-1}}||u||_{H^{s}}\right)$
$\displaystyle+C\left(\dot{a}||\omega||_{H^{s+1}}(1+||u||_{H^{s}})+(1+||\omega||_{H^{s+1}})||k^{tr}||_{H^{s-1}}\right)$

which for sufficiently small $||u||_{H^{s}}+||\gamma-\gamma^{*}||_{H^{s}}$
yields

$\displaystyle||X||_{H^{s+1}}\leq
C\left(\dot{a}||\omega||_{H^{s+1}}+||\omega||_{H^{s+1}}||k^{tr}||_{H^{s-1}}+a^{2-n\gamma_{ad}}(1+||\delta\rho||_{H^{s-1}})\right.$
$\displaystyle\left.(||v||^{2}_{H^{s-1}}(1+||u||_{H^{s}})||v||_{H^{s-1}}+||v||_{H^{s-1}})+||k^{tr}||_{H^{s-1}}||u||_{H^{s}}\right)$
$\displaystyle+C\left(\dot{a}||\omega||_{H^{s+1}}(1+||u||_{H^{s}})+||k^{tr}||_{H^{s-1}}\right).~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\Box$

Here several terms contain the factor $\dot{a}$ which may be potentially