Einstein flow with matter sources: stability and convergence

Viewed on a sufficiently coarse-grained scale the portion of our universe that is accessible to observation appears to be spatially homogeneous and isotropic. If, as is usually imagined, one should be able to extrapolate these features to (a suitably coarse-grained model of) the universe as a whole then only a handful of spatial manifolds need be considered in cosmology — the familiar Friedmann-Lemaître-Robertson-Walker (FLRW) archetypes of constant positive, vanishing or negative curvature [1, 3]. These geometries consist, up to an overall, time-dependent scale factor, of the 3-sphere, $\mathbb{S}^{3}$, with its canonical ‘round’ metric, Euclidean 3-space, $\mathbb{E}^{3}$, hyperbolic 3-space, $\mathbb{H}^{3}$ and the quotient space $\mathbb{RP}(3)=\mathbb{S}^{3}/\mathbb{Z}_{2}$ obtainable from $\mathbb{S}^{3}$ by the identification of antipodal points [4]. Of these possibilities only the sphere and its 2-fold quotient $\mathbb{RP}^{3}$ are closed and thus compatible with a universe model of finite extent. It is not known of course whether the actual universe is spatially closed or not but, to simplify the present discussion, we shall limit our attention herein to models that are. More precisely we shall focus on spacetimes admitting Cauchy hypersurfaces that are each diffeomorphic to a smooth, connected 3-manifold that is compact, orientable and without boundary.

On the other hand if one takes literally the cosmological principle that only manifolds supporting a globally homogeneous and isotropic metric should be considered in models for the actual universe then, within the spatially compact setting considered here, only the 3-sphere and $\mathbb{RP}^{3}$ would remain. But the astronomical observations which motivate this principle are necessarily limited to a (possibly quite small) fraction of the entire universe and are compatible with models admitting metrics that are only locally, but not necessarily globally, spatially homogeneous and isotropic. As is well-known there are spatially compact variants of all of the basic Friedmann-Lemaître- Robertson-Walker cosmological models, mathematically constructable (in the cases of vanishing or negative curvature) by taking suitable compact quotients of Euclidean 3-space $\mathbb{E}^{3}$ or of hyperbolic 3-space $\mathbb{H}^{3}$. One can also take infinitely many possible quotients of $\mathbb{S}^{3}$ to obtain the so-called spherical space forms that are locally compatible with the FLRW constant positive curvature geometry but are no longer diffeomorphic to the 3-sphere.

Still more generally though we shall find that there is a dynamical mechanism at work within the Einstein ‘flow’, suitably viewed in terms of the evolution of 3-manifolds to develop 4-dimensional, globally hyperbolic spacetimes, and extended to include suitable matter sources and a positive cosmological constant $\Lambda$, that strongly suggests that even manifolds that do not admit a locally homogeneous and isotropic metric at all will nevertheless evolve in such a way as to be asymptotically compatible with the observed homogeniety and isotropy. This reflects an argument which we shall sketch that, under Einsteinian evolution, the summands making up $M$ (in a connected sum decomposition) that do support locally homogeneous and isotropic metrics will tend to overwhelmingly dominate the spatial volume asymptotically as the universe model continues to expand and furthermore that the actual evolving (inhomogeneous, non-isotropic) metric on M will naturally tend to flow towards a homogeneous, isotropic one on each of these asymptotically volume-dominating summands.

We do not claim that this mechanism is yet so compelling, either mathematically or physically, as to convince one that the actual universe has a more exotic topology but only that such a possibility is not strictly excluded by current observations. However, it is intriguing to investigate the possibility that there may be a dynamical reason, provided by Einstein’s equations, for the observed fact that the universe seems to be at least locally homogeneous and isotropic and that this mechanism may therefore allow an attractive logical alternative to simply extrapolating observations of necessarily limited scope to the universe as a whole.

But what are the (compact, connected, orientable) 3-manifolds available for consideration? This question has been profoundly clarified in recent years by the dramatic progress on lower dimensional topology made possible through the advancements in Ricci flow [5, 17]. One now knows for example that, since the Poincaré conjecture has finally been proven, any such 3-manifold M that is in fact simply connected must be diffeomorphic to the ordinary 3-sphere $\mathbb{S}^{3}$. Setting aside this so-called ‘trivial’ manifold the remaining possibilities consist of an infinite list of nontrivial manifolds, each of which is diffeomorphic (designated herein by $\approx$) to a finite connected sum of the following form:

	$\displaystyle M\approx$		(1)
	$\displaystyle\mathbf{S}^{3}/\Gamma_{1}\#..\#\mathbf{S}^{3}/\Gamma_{k}\#(\mathbf{S}^{2}\times\mathbf{S}^{1})_{1}\#..\#(\mathbf{S}^{2}\times\mathbf{S}^{1})_{l}\#\mathbf{K}(\pi,1)_{1}\#..\#\mathbf{K}(\pi,1)_{m}.$

Here $k,l$ and $m$ are integers $\geq 0,k+l+m\geq 1$ and if either $k$, $l$ or $m$ is 0 then terms of that type do not occur. The connected sum $M\#N$ of two closed connected, oriented n-manifolds is constructed by removing the interiors of an embedded closed n-ball in each of $M$ and $N$ and then identifying the resulting $\mathbb{S}^{n-1}$ boundary components by an orientation-reversing diffeomorphism. The resulting $n-$manifold will be smooth, connected, closed and consistently oriented with the original orientations of $M$ and $N$. The above decomposition of $M$ is only uniquely defined provided we set aside $\mathbb{S}^{3}$ since $M\#\mathbb{S}^{3}\approx M$ for any 3-manifold $M$.

In the above formula if $k\geq 1$, then each $\Gamma_{i}$, $1\leq i\leq k$ is a finite, nontrivial ($\Gamma_{i}\neq I$) subgroup of SO(4) acting freely and orthogonally on $\mathbb{S}^{3}$. The individual summands $\mathbb{S}^{3}/\Gamma_{i}$ are the spherical space forms alluded to previously and, by construction, each is compatible with an FLRW metric of constant positive spatial curvature (i.e., $k=+1$ models in the usual notation). The individual ‘handle’ summands $\mathbb{S}^{2}\times\mathbb{S}^{1}$ admit metrics of the Kantowski-Sachs type that are homogeneous but not isotropic and so not even locally of FLRW type.

The remaining summands in the above ‘prime decomposition’ theorem [6] are the $K(\pi,1)$ manifolds of Eilenberg-MacLane type wherein, by definition $\pi=\pi_{1}(M)$, the fundamental group of $M$ and all of the higher homotopy groups are trivial, that is $\pi_{i}(M)=0$ for $i>1$. Equivalently, the universal covering space of $M$ is contractible and, in this case, known to be diffeomorphic to $\mathbb{R}^{3}$ [6]. Since the higher homotopy groups, $\pi_{i}(M)$ for $i>1$ can be interpreted as the homotopy classes of continuous maps $\mathbb{S}^{i}\to M$, each such map must be homotopic to a constant map. For this reason $K(\pi,1)$ manifolds are said to be aspherical.

This general class of $K(\pi,1)$ manifolds includes, as special cases, the 3-torus and five additional manifolds, finitely covered by the torus, that are said to be of ‘flat type’ since they are the only compact, connected, orientable 3-manifolds that each, individually, admits a flat metric and thus supports spatially compactified versions of the FLRW spaces of flat type (i.e., $k=0$ models).

Other $K(\pi,1)$ spaces include the vast set of compact hyperbolic manifolds $\mathbb{H}^{3}/\Gamma$, where here $\Gamma$ is a discrete torsion-free (i.e., no nontrivial element has finite order) co-compact subgroup of the Lie group Isom⁺($\mathbb{H}^{3}$) of orientation preserving isometries of $\mathbb{H}^{3}$ that, in fact, is Lie-group isomorphic to the proper orthochronous Lorentz group SO⁺(3, 1). Each of these, individually, supports spatially compactified versions of the FLRW spacetimes of constant negative (spatial) curvature (i.e., $k=-1$ models).

Additional $K(\pi,1)$ manifolds include the trivial circle bundles over higher genus surfaces $\Sigma_{p}$ for $p\geq 2$ (where $\Sigma_{p}$ designates a compact, connected, orientable surface of genus $p$) and nontrivial circle bundles over $\Sigma_{p}$ for $p\geq 1$. Note that the trivial circle bundles $\mathbb{S}^{2}\times\mathbb{S}^{1}$ and $\mathbb{T}^{2}\times\mathbb{S}^{1}\approx\mathbb{T}^{3}$ are already included among the previous prime factors discussed and that nontrivial circle bundles over $\mathbb{S}^{2}$ are included among the spherical space forms $\mathbb{S}^{3}/\Gamma$ for suitable choices of $\Gamma$. Still further examples of $K(\pi,1)$ manifolds are compact 3-manifolds that fiber nontrivially over the circle with fiber $\Sigma_{p}$ for $p\geq 1$. Any such manifold is obtained by identifying the boundary components of $[0,1]\times\Sigma_{p}$ with a (nontrivial) orientation-reversing diffeomorphism of $\Sigma_{p}$.

It is known however that every prime $K(\pi,1)$ manifold is decomposible into a (possibly trivial but always finite) collection of (complete, finite volume) hyperbolic and graph manifold components. The possibility of such a (nontrivial) decomposition arises whenever the $K(\pi,1)$ manifold under study admits a nonempty family $\{T_{i}\}$ of disjoint embedded incompressible two-tori. An embedded two-torus $\mathbb{T}^{2}$ is said to be incompressible if every incontractible loop in the torus remains incontractible when viewed as a loop in the ambient manifold. A closed oriented 3-manifold G (possibly with boundary) is a graph manifold if there exists a finite collection $\{T^{{}^{\prime}}_{i}\}$ of disjoint embedded incompressible tori $\{T^{{}^{\prime}}_{i}\}\subset$G such that each component $G_{j}$ of $G-\cup T^{{}^{\prime}}_{i}$ is a Seifert-fibered space (A Seifert-fibered space is a 3-manifold foliated by circular fibers in such a way that each fiber has a tubular neighborhood (characterized by a pair of co-prime integers) of the special type known as a standard fibered torus). Thus a graph manifold is a union of Seifert-fibered spaces glued together by toral automorphisms along toral boundary components. The collection of tori is allowed to be empty so that, in particular, a Seifert-fibered manifold itself is a graph manifold. Decomposing a 3-manifold by cutting along essential two-spheres (to yield its prime factors) and then along incompressible tori, when present, are the basic operations that reduce a manifold to its ‘geometric’ constituents [6]. The Thurston conjecture that every such 3-manifold can be reduced in this way has now been established via arguments employing Ricci flow [5, 17].

It may seem entirely academic to consider such general, ‘exotic’ 3-manifolds as the composite (i.e., nontrivial connected sum) ones described above as arenas for general relativity when essentially all of the explicitly known solutions of Einstein’s equations (in this spatially compact setting) involve only individual, ‘prime factors’. As we shall see however some rather general conclusions are derivable concerning the behaviors of solutions to the field equations on such exotic manifolds and astronomical observations do not logically exclude the possibility that the actual universe could have such a global topological structure. It is furthermore conceivable that the validity of central open issues in general relativity like the cosmic censorship conjecture could depend crucially upon the spatial topology of the spacetime under study.

In our study, the spacetime is described by an ‘n+1’ dimensional globally hyperbolic Lorentzian manifold $\tilde{M}\approx$ $\mathbb{R}\times M$ with each level set $\{t\}\times M$ of the time function $t$ being an orientable n-manifold diffeomorphic to a Cauchy hypersurface and equipped with a Riemannian metric. Throughout our analysis, we will consider $M$ to be of negative Yamabe type (see appendix for detailed description) since this will suffice to ensure that an expanding universe model can never achieve a maximal hypersurface and then ‘recollapse’. Such a split $\mathbb{R}\times M$ may be executed by introducing a lapse function $N$ and shift vector field $X$ belonging to suitable function spaces and defined such that

$\displaystyle\partial_{t}$

$\displaystyle=$

$\displaystyle N\hat{n}+X$

(2)

with $t$ and $\hat{n}$ being the time function and a hypersurface orthogonal future directed timelike unit vector, respectively (so that in particular $\hat{g}(\hat{n},\hat{n})=-1$). The above splitting expresses the spacetime metric $\tilde{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)$

(3)

and the stress-energy tensor as

$\displaystyle\mathbf{T}$

$\displaystyle=$

$\displaystyle E\mathbf{n}\otimes\mathbf{n}+\mathbf{J}\otimes\mathbf{n}+\mathbf{n}\otimes\mathbf{J}+\mathbf{S},$

(4)

where $\mathbf{J}\in\mathfrak{X}(M)$, $\mathbf{S}\in S^{2}(M)$. Here, $\mathfrak{X}(M)$ and $S^{2}(M)$ are the space of vector fields and the space of symmetric covariant 2-tensors, respectively. Under such an $n+1$ decomposition, the Einstein equations,

$\displaystyle R_{\mu\nu}-\frac{1}{2}Rg_{\mu\nu}+\Lambda g_{\mu\nu}$

$\displaystyle=$

$\displaystyle T_{\mu\nu},$

(5)

take the form (in units for which $8\pi G=c=1$) of evolution equations,

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

along with constraints (Gauss and Codazzi equations)

	$\displaystyle R(g)-|k|^{2}+\tau^{2}$	$\displaystyle=$	$\displaystyle 2\Lambda+2E,$		(8)
	$\displaystyle\nabla_{j}k^{j}_{i}-\nabla_{i}\tau$	$\displaystyle=$	$\displaystyle-J_{i},$		(9)

where $S=g^{ij}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 motions of the matter

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

Here we impose ‘constant mean extrinsic curvature’ (CMC) as the temporal gauge condition which yields an elliptic equation for the lapse function (notice that the lapse and shift do not have evolution equations; a consequence of the diffeomorphism invariance of the Einstein equations). Later on, we will select a suitable, complementary spatial gauge. CMC gauge entails

$\displaystyle\tau=tr_{g}k=monotonic~{}function~{}of~{}t~{}alone$

(11)

so $\tau$ is thus constant throughout each $t=constant$ hypersurface and therefore can play the role of time. Using the evolution and constraint equations, one may derive the following equation for the lapse function

$\displaystyle\frac{\partial\tau}{\partial t}$

$\displaystyle=$

$\displaystyle\Delta_{g}N+\{|k|^{2}+\frac{S}{n-1}+\frac{n-2}{n-1}E-\frac{2\Lambda}{n-1}\}N+L_{X}\tau$

(12)

which after implementing CMC gauge ($\partial_{i}\tau=0$) yields

$\displaystyle\frac{\partial\tau}{\partial t}$

$\displaystyle=$

$\displaystyle\Delta_{g}N+\{|k|^{2}+\frac{S}{n-1}+\frac{n-2}{n-1}E-\frac{2\Lambda}{n-1}\}N,$

where $|k|^{2}=k_{ij}k^{ij}$ and the Laplacian is defined as $\Delta_{g}=-\nabla[g]^{i}\nabla[g]_{i}$ and therefore has positive spectrum on compact connected manifolds.

Utilizing the ‘CMC’ condition, the momentum constraint (9) simplifies to

$\displaystyle\nabla[g]_{j}k^{j}_{i}$

$\displaystyle=$

$\displaystyle-J_{i},$

(13)

the solution of which may be written as

$\displaystyle k^{i}_{j}=K^{tri}_{j}+\frac{\tau}{n}\delta^{i}_{j}$

(14)

with $K^{tr}$ being traceless with respect to $g$ i.e., $K^{tr}_{ij}g^{ij}=0$ (and so with respect to any metric conformal to $g$). Note that $K^{tr}$ is obtained by solving the following equation

$\displaystyle\nabla[g]_{j}K^{trj}_{i}$

$\displaystyle=$

$\displaystyle-J_{i}.$

(15)

The Hamiltonian constraint (8) may be written using the solution of the momentum constraint (13) as follows

$\displaystyle R(g)$

$\displaystyle=$

$\displaystyle|K^{tr}|^{2}+2E-\frac{n-1}{n}(\tau^{2}-\frac{2n\Lambda}{n-1}).$

Since we are primarily interested in the case $\Lambda>0$ (the simplest model for ‘dark energy’), it might appear that $\tau^{2}-\frac{2n\Lambda}{n-1}$ could be negative. But then the Hamiltonian constraint would imply that if the energy density $E$ is non-negative, then $R(g)\geq 0$ everywhere on $M$ which is impossible for a manifold of negative Yamabe type (see appendix for details about Yamabe classification). Let’s consider the energy condition and establish an allowed range for the constant extrinsic mean curvature $\tau$ which is now playing the role of time. The weak energy condition yields

$\displaystyle\mathbf{T}(\mathbf{n},\mathbf{n})\geq 0,~{}E\geq 0$

(16)

so that to any time-like observer the energy density is nonnegative. Physically relevant classical matter sources are expected to satisfy this energy condition. We will only consider matter sources with point-wise nonnegative energy density throughout the spacetime. Therefore a universe filled with matter satisfying the weak energy condition will always have

$\displaystyle\tau^{2}-\frac{2n\Lambda}{n-1}>0,$

(17)

and for expanding models (so that $\frac{\partial\tau}{\partial t}>0$),

$\displaystyle-\infty<\tau<-\sqrt{\frac{2n\Lambda}{n-1}}.$

(18)

Since we are primarily interested in the asymptotic behavior of an ‘expanding’ universe model (the physically relevant case), we set the range of $\tau$ to be $(-\infty,-\sqrt{\frac{2n\Lambda}{n-1}})$ once and for all. We turn our attention to the lapse equation in an expanding universe model by setting (to specify the time gauge precisely)

$\displaystyle\frac{\partial\tau}{\partial t}$

$\displaystyle=$

$\displaystyle\frac{n}{2(n-1)}(\tau^{2}-\frac{2n\Lambda}{n-1})^{\frac{n}{2}}>0$

(19)

whose solution $\tau=\tau(t)$ plays the role of time from now onwards. The lapse equation for this choice of time co-ordinate is expressible as

	$\displaystyle\Delta_{g}N+(|K^{tr}|^{2}+(\frac{\tau^{2}}{n}-\frac{2\Lambda}{n-1})+\frac{1}{n-1}(S+(n-2)E))N$		(20)
	$\displaystyle=\frac{n}{2(n-1)}(\tau^{2}-\frac{2n\Lambda}{n-1})^{\frac{n}{2}}.$

a unique solution of which is implied by

$\displaystyle(n-2)E+S\geq 0.$

(21)

Note that matter satisfying the strong energy condition

	$\displaystyle(T_{\mu\nu}-\frac{1}{2}Tg_{\mu\nu})n^{\mu}n^{\nu}\geq 0,$		(22)
	$\displaystyle E+S\geq 0$

automatically satisfies the condition (21) for $n\geq 3$ since we are considering $E\geq 0$. In a sense the matter fields satisfying both weak and strong energy conditions form a subset in the space of mater fields satisfying $E\geq 0$ and $S+(n-2)E\geq 0$ for $n\geq 3$ and therefore such sources are allowed for our analysis. Of course, most physically relevant sources do satisfy the property of non-negative energy density (weak energy condition) and the attractive nature of gravity (strong energy condition). Several known sources of physical interest satisfy both the weak and strong energy conditions. These include for example perfect fluids and Vlasov matter (see [23, 21, 19, 7, 22] for details about Vlasov matter).

We have established the existence of a unique solution of the lapse equation provided that the matter sources in the universe satisfy a suitable energy condition. A standard maximum principle argument for the corresponding elliptic equation yields the following estimate for the lapse function

	$\displaystyle 0$	$\displaystyle<$	$\displaystyle\frac{\frac{n}{2(n-1)}(\tau^{2}-\frac{2n\Lambda}{n-1})^{\frac{n}{2}}}{\left((\frac{\tau^{2}}{n}-\frac{2\Lambda}{n-1})+\sup(|K^{tr}|^{2}+\frac{1}{n-1}(S+(n-2)E))\right)}\leq N\leq$
			$\displaystyle\frac{n^{2}}{2(n-1)}(\tau^{2}-\frac{2n\Lambda}{n-1})^{\frac{n}{2}-1}.$

The important thing here is to note that the lapse function is positive for an expanding universe model. The results obtained so far will be sufficient to study the dynamical behavior in terms of a weak Lyapunov functional.

A Lyapunov functional is indispensable while studying the stability properties of a dynamical system. Construction of a Lyapunov functional is however a challenging task. In order to study the dynamical properties of the Einstein-matter system (possibly including a positive cosmological constant), we require a Lyapunov functional for this field theoretic problem. Using a conformal technique, [18, 2] constructed a reduced phase space as the cotangent bundle of the higher dimensional analogue of the Teichmüller space of $2+1$ dimensional gravity and obtained the following true Hamiltonian of the dynamics³³3 wherein $g$ is expressed in terms of a ‘conformal’ metric $\gamma$ and a conformal factor $\varphi$ that satisfies the associated Lichnerowicz equation

$\displaystyle H_{reduced}$

$\displaystyle:=$

$\displaystyle\frac{2(n-1)}{n}\int_{M}\frac{\partial\tau}{\partial t}\mu_{g}.$

(24)

In the particular case of a vacuum limit, it is indeed possible to construct such a reduced Hamiltonian in ‘CMC’ gauge, which also acts as a Lyapunov functional. Motivated by their study, we consider the rescaled volume functional as a Lyapunov functional and call it $L(g,k)$

$\displaystyle L(g,k)$

$\displaystyle=$

$\displaystyle\frac{2(n-1)}{n}\int_{M}\frac{\partial\tau}{\partial t}\mu_{g},$

(25)

where $g$ is expressed as $g_{ij}=\varphi^{4/n-2}\gamma_{ij},~{}\varphi:M\to\mathbb{R}_{>0},~{}R(\gamma)=-1$. We call $L(g,k)$ a weak Lyapunov functional because, following the expression of $\frac{\partial\tau}{\partial t}$ and the corresponding Lichnerowicz equation for $\varphi$, it controls only the $H^{1}\times L^{2}$ norm of the data $(g,k)$ while the desired norm would more likely be $H^{s}\times H^{s-1},s>\frac{n}{2}+\epsilon$ for some $\epsilon\geq 1$. Note that we do not, at this point, have a local existence theorem for the Cauchy problem of the Einstein system with arbitrary matter sources satisfying the desired energy conditions. At the moment, let us focus on the time evolution of this weak Lyapunov functional. Due to its weak character, we will not be able to state a theorem concerning the stability (either Lyapunov or asymptotic) of the spacetime on the basis of its time evolution. Nevertheless, we will be able to obtain some important physical results related to the asymptotic behavior of the spacetime (in the expanding direction) utilizing the limiting behavior of the Lyapunov function. The time derivative of this functional may be obtained as

	$\displaystyle\frac{dL(g,k)}{dt}$		(26)
	$\displaystyle=n\tau(\tau^{2}-\frac{2n\Lambda}{n-1})^{\frac{n}{2}-1}\int_{M}(|K^{tr}|^{2}+\frac{1}{n-1}(S+(n-2)E))N\mu_{g},$		(27)

where we have used the field equations (6-7) and the identity

	$\displaystyle\frac{\partial^{2}\tau}{\partial t^{2}}$	$\displaystyle=$	$\displaystyle\frac{n}{2(n-1)}\tau(\tau^{2}-\frac{2n\Lambda}{n-1})^{\frac{n}{2}-1}\left(n\Delta_{g}N+n|K^{tr}|^{2}N\right.$
			$\displaystyle\left.+(\tau^{2}-\frac{2n\Lambda}{n-1})N+\frac{n}{n-1}(S+(n-2)E)N\right),$

which was obtained using the lapse equation and the CMC gauge condition ($\partial_{i}\tau=0$). We have also used Stokes’ theorem to eliminate the covariant divergence terms in the integral over the compact manifold $M$. Along the solution curve in the expanding direction ($\frac{\partial\tau}{\partial t}>0$ and $-\infty<\tau<-\sqrt{\frac{2n\Lambda}{n-1}}$), therefore, the Lyapunov functional monotonically decays i.e.,

$\displaystyle\frac{dL(g,k)}{dt}<0$

(29)

and can only attain its infimum (i.e., $\frac{dL}{dt}=0$) if the following conditions are met (perhaps only asymptotically)

$\displaystyle K^{tr}$

$\displaystyle=$

$\displaystyle 0,~{}S+(n-2)E=0.$

(30)

Substitution of the first condition into the momentum constraint (13) immediately yields

$\displaystyle J^{i}\equiv 0$

(31)

everywhere on $M$. In addition to satisfying $E,S+(n-2)E\geq 0$, if the matter sources also satisfy the physically relevant strong energy condition i.e., $(S+E)\geq 0$, then for $n\geq 3$ (the cases of primary interest), the second condition for the infimum of the Lyapunov function translates to

$\displaystyle S+(n-2)E=0,~{}E=0$

(32)

and therefore

$\displaystyle S=0.$

(33)

This result, therefore, states that the weak Lyapunov functional is monotonically decaying in the direction of cosmological expansion and approaches its infimum only in the limit that the matter sources be asymptotically ‘turned off’. In the limit of such ‘extinguished’ matter, one may compute $\frac{d^{2}L}{dt^{2}}=0$ and $\frac{d^{3}L}{dt^{3}}=2n\tau(\tau^{2}-\frac{2n\Lambda}{n-1})^{\frac{n}{2}-1}\int_{M}|\partial_{t}K^{tr}|^{2}\mu_{g}<0$ unless $\partial_{t}K^{tr}=0$ as well. At this limit, one may simplify the evolution and constraint equations by substituting $K^{tr}=\partial_{t}K^{tr}=0$ to obtain the background warped product spacetimes

	$\displaystyle\hat{g}$		(34)
	$\displaystyle=-\frac{n^{2}}{(\tau^{2}-\frac{2n\Lambda}{n-1})^{2}}d\tau\otimes d\tau+\frac{n}{(n-1)(\tau^{2}-\frac{2n\Lambda}{n-1})}\gamma_{ij}dx^{i}\otimes dx^{j},\tau\in(-\infty,-\frac{2n\Lambda}{n-1}),$

with $Ric(\gamma)=-\frac{1}{n}\gamma$. Each of these spacetimes admits a globally defined future directed time like conformal Killing field $Y$ given by $Y:=Y^{\mu}\partial_{\mu}=\left(\tau^{2}-\frac{2n\Lambda}{n-1}\right)^{1/2}\partial_{\tau}$ with $L_{Y}\hat{g}=\frac{-2\tau}{(\tau^{2}-\frac{2n\Lambda}{n-1})^{1/2}}\hat{g}$. Therefore we designate these spacetimes as ‘conformal’ spacetimes. Notice that if one turns off the cosmological constant, one retrieves the so called ‘Lorentz cone’ spacetimes analyzed in [13]

$\displaystyle\hat{g}=-\frac{n^{2}}{\tau^{4}}d\tau\otimes d\tau+\frac{n}{(n-1)\tau^{2}}\gamma_{ij}dx^{i}\otimes dx^{j},\tau\in(-\infty,0),$

(35)

which in the ‘$3+1$’ dimensional case is also known as the ’Milne’ spacetime. In this special case $Y$ reduces to the homothetic Killing field $Y=-\tau\partial_{\tau}$ with $L_{Y}\hat{g}=2\hat{g}$. The vital question is whether the Lyapunov functional $L$ ever attains its infimum. In an expanding universe model, if the matter sources do not re-collapse to form a singularity, then the matter density falls off and, in the limiting case, may be considered to be negligible. Observing the monotonic decay of the Lyapunov functional along a solution curve, one may be tempted to conjecture the asymptotic stability of a matter filled spacetime. However, we remind the reader again that such a property only provides a weak notion of stability if the matter sources do not develop singularities. Even a vacuum spacetime might be able to go singular whereby the pure gravitational energy density can collapse to form a singularity before the spatial volume of the universe reaches infinity. These questions are related to the Cosmic Censorship conjecture and are part of active research.

A physical cosmology is often built on the assumption that a perfect fluid is the source for Einstein’s equations. It is well known that perfect fluids possess the pathological property of forming shock singularities in finite time (such a result was established by D. Christodoulou [20] for the case of a perfect fluid evolving on the Minkowski spacetimes). However, we will not address such issues here but rather will observe in a later section that on a sufficiently rapidly expanding spacetime one may avoid shocks by imposing a certain smallness condition on the data. The perfect fluid stress-energy tensor

$\displaystyle\mathbf{T}^{\mu\nu}$

$\displaystyle=$

$\displaystyle(P+\rho)u^{\mu}u^{\nu}+P\hat{g}^{\mu\nu}$

(36)

yields

	$\displaystyle E=\mathbf{T}(\mathbf{n},\mathbf{n})=(P+\rho)(Nu^{0})^{2}-P\geq 0,$		(37)
	$\displaystyle S_{ij}:=\mathbf{T}(\partial_{i},\partial_{j})=(P+\rho)u_{i}u_{j}+Pg_{ij}$		(38)

with the equation of state $P=(\gamma_{a}-1)\rho$. Now assuming $1\leq\gamma_{a}\leq 2$ and a positive mass energy density $\rho$, we obtain that the pressure satisfies $0\leq P\leq\rho$. We use the $n+1$ split of the velocity $u=v-g(u,\mathbf{n})\mathbf{n}$ to express the entities in terms of the spatial ($M$ tangential) velocity $v$. The energy density and momentum flux density tensor $S_{ij}$ read

$\displaystyle E=(P+\rho)(1+g_{ij}v^{i}v^{j})-P,~{}S_{ij}=(P+\rho)v_{i}v_{j}+Pg_{ij}.$

(39)

The energy conditions that are required by our analysis

$\displaystyle E\geq 0,~{}(n-2)E+S\geq 0$

(40)

take the forms of the following inequalities in the case of a perfect fluid

	$\displaystyle(P+\rho)g_{ij}v^{i}v^{j}+\rho\geq 0,$		(41)
	$\displaystyle(n-1)(P+\rho)g_{ij}v^{i}v^{j}+nP+(n-2)\rho\geq 0$		(42)

which are trivially satisfied under the assumption of $\rho\geq 0$ and $\gamma_{a}\in[1,2]$. In order to extract more dynamical information from our weak Lyapunov functional, we use a conformal technique similar to [13]. Using the momentum constraint, the second fundamental form may be written as follows (14)

$\displaystyle K^{ij}=K^{trij}+\frac{\tau}{n}g^{ij},$

(43)

where $g_{ij}K^{trij}=0$, but $\nabla_{j}K^{trij}=-J^{i}$ under ths CMC condition ($\partial_{i}\tau=0$). We use the conformal transformation as described in [13, 2, 1] only replacing the momenta conjugate to the metric $g$ by the second fundamental form through the inverse of a Legendre transformation $\pi=-\mu_{g}(k-(tr_{g}k)g)$. The conformal transformation reads

$\displaystyle(g_{ij},K^{trij})=(\varphi^{\frac{4}{n-2}}\gamma_{ij},\varphi^{-\frac{2(n+2)}{n-2}}\kappa^{trij}),$

(44)

where $\gamma$ and $\kappa^{tr}$ are the scale-free fields satisfying $\gamma^{ij}\kappa^{tr}_{ij}=0$ with $\gamma\in\mathcal{M}_{-1}$ and $\varphi:M\to\mathbf{R}_{>0}$. Here $\mathcal{M}_{-1}$ is defined as $\mathcal{M}_{-1}=\{\gamma$ is a Riemannian metric on M $|R(\gamma)=-1\}$. In reality, the fully reduced dynamics assumes aan equivalence class lying in the orbit space $\mathcal{M}_{-1}/D_{0}$, $D_{0}$ being the group of diffeomorphisms of $M$ isotopic to the identity. To avoid technical complexities, the calculations will simply be restricted to $\mathcal{M}_{-1}$ as the entities we are interested in (such as $\int_{M}\mu_{\gamma}=\int_{M}\sqrt{\det{\gamma_{ij}}}dx^{1}\wedge dx^{2}...\wedge dx^{n}$) are $D_{0}$ invariant. Note that $K^{tr}$ may be decomposed into a transverse-traceless part (with respect to $g$) and a conformal Killing tensor part

$\displaystyle K^{trij}$

$\displaystyle=$

$\displaystyle K^{TTij}+\varphi^{\frac{-2n}{n-2}}(L_{Y}g-\frac{2}{n}\nabla_{m}Y^{m}g)^{ij},$

(45)

where $Y\in\mathfrak{X}(M)$. Under this conformal transformation, the components of the conformal Killing tensor $(L_{g}Y-\frac{2}{n}\nabla_{m}Y^{m}g)$ transform as follows

$\displaystyle(L_{g}Y-\frac{2}{n}\nabla_{m}Y^{m}g)^{ij}=\varphi^{\frac{-4}{n-2}}(L_{\gamma}Y-\frac{2}{n}\nabla_{m}Y^{m}\gamma)^{ij}$

(46)

yielding a consistent transformation of the transverse-traceless tensor $K^{TT}$, that is, $K^{TTij}=\varphi^{-\frac{2(n+2)}{n-2}}\kappa^{TTij}$. Therefore, we may write

$\displaystyle\kappa^{trij}=\kappa^{TT}+(L_{\gamma}Y-\frac{2}{n}\nabla_{m}Y^{m}\gamma)^{ij}.$

(47)

Under this conformal transformation, the Hamiltonian constraint (8) becomes the following semi-linear elliptic equation for the conformal factor

	$\displaystyle\Delta_{\gamma}\varphi-\frac{n-2}{4(n-1)}\varphi-\frac{(n-2)}{4(n-1)}\varphi^{\frac{-3n+2}{n-2}}|\kappa^{TT}+(L_{Y}\gamma-\frac{2}{n\mu_{\gamma}}\nabla_{m}Y^{m}\gamma)|^{2}$		(48)
	$\displaystyle-\frac{(n-2)E_{g}}{2(n-1)}\varphi^{\frac{n+2}{n-2}}+\frac{n-2}{4n}(\tau^{2}-\frac{2n\Lambda}{n-1})\varphi^{\frac{n+2}{n-2}}=0,$

where $E_{g}$ is the matter energy density that needs scaling under the conformal transformation (44). We will scale the matter degrees of freedom in the following way. First, the momentum constraint after the re-scaling of the traceless second fundamental form reads

$\displaystyle\nabla[\gamma]_{j}\kappa^{trij}=-\varphi^{\frac{2(n+2)}{n-2}}J_{g}^{i},$

(49)

where $J_{g}^{i}$ is the momentum density without conformal scaling. Now of course, if one chooses the following scaling for the momentum density (York scaling [25])

$\displaystyle J_{g}^{i}=\varphi^{\frac{-2(n+2)}{n-2}}J_{\gamma}^{i},$

(50)

the momentum constraint becomes decoupled from the Lichnerowicz equation (48) i.e.,

$\displaystyle\nabla[\gamma]_{j}\kappa^{trij}=-J_{\gamma}^{i}.$

(51)

Now upon substituting $\kappa^{tr}$ from equation (47) into equation (51) one arrives at an elliptic equation for $Y$

$\displaystyle-\Delta_{\gamma}Y^{i}+R[\gamma]^{i}_{m}Y^{m}+(1-\frac{2}{n})\nabla[\gamma]^{i}(\nabla[\gamma]_{m}Y^{m})$

$\displaystyle=$

$\displaystyle J[\gamma]^{i}.$

(52)

The vector field $Y$ therefore depends solely on the metric $\gamma$ and the matter source $J_{\gamma}$, not on the conformal factor $\varphi$. Now, utilizing the normalization condition $g_{\mu\nu}u^{\mu}u^{\nu}=-1$, we obtain

$\displaystyle-(Nu^{0})^{2}+g_{ij}v^{i}v^{j}=-1.$

(53)

Since, the right hand side of this equation is a constant (and therefore conformally invariant), each term on the left hand side may be chosen to be conformally invariant leading to the following scalings of $Nu^{0}$ and $v^{i}$

$\displaystyle(Nu^{0})_{g}$

$\displaystyle=$

$\displaystyle(Nu^{0})_{\gamma},v^{i}_{g}=\varphi^{-\frac{2}{n-2}}v^{i}_{\gamma},$

(54)

which, together with the expressions for matter energy and momentum density, yields the following scalings for $P,\rho,$ and $E_{g}$

$\displaystyle(P+\rho)_{g}$

$\displaystyle=$

$\displaystyle\varphi^{-\frac{2(n+1)}{n-2}}(P+\rho)_{\gamma},~{}E_{g}=\varphi^{-\frac{2(n+1)}{n-2}}E_{\gamma}.$

(55)

Notice that these scalings are only valid for a perfect fluid source since we have explicitly made use of the normalization relation (53). Rescaling of the matter energy density $E_{g}$ yields the following Lichnerowicz equation

	$\displaystyle\Delta_{\gamma}\varphi-\frac{n-2}{4(n-1)}\varphi-\frac{(n-2)}{4(n-1)}\varphi^{\frac{-3n+2}{n-2}}|\kappa^{TT}+(L_{Y}\gamma-\frac{2}{n\mu_{\gamma}}\nabla_{m}Y^{m}\gamma)|^{2}$		(56)
	$\displaystyle-\frac{(n-2)E_{\gamma}}{2(n-1)}\varphi^{-\frac{n}{n-2}}+\frac{n-2}{4n}(\tau^{2}-\frac{2n\Lambda}{n-1})\varphi^{\frac{n+2}{n-2}}=0,$

a unique solution of which may be obtained either by the standard sub and super solution technique of [25, 26] or by a direct method such as the variational technique of [27]. Note that the exponent of $\varphi$ in the term $\frac{(n-2)E_{\gamma}}{2(n-1)}\varphi^{-\frac{n}{n-2}}$ is crucial in proving the uniqueness and existence of solutions to the Lichnerowicz equation and therefore, the rescaling of the energy density $E_{g}$ is necessary. A straightforward maximum principle argument, applied to the Lichnerowicz equation, shows that the unique positive solution $\varphi=\varphi(\tau,\gamma,\kappa^{TT},Y,E_{\gamma})$ satisfies

$\displaystyle\varphi^{\frac{4}{n-2}}\geq\frac{n}{n-1}\frac{1}{(\tau^{2}-\frac{2n\Lambda}{n-1})}$

(57)

with equality holding everywhere on $M$ iff

$\displaystyle\kappa^{tr}=\kappa^{TT}+\{L_{Y}\gamma-\frac{2}{n}\nabla_{m}Y^{m}\gamma\}\equiv 0,E_{\gamma}\equiv 0$

(58)

on $M$. Therefore the asymptotic analysis becomes straightforward as the infimum of the Lyapunov functional is attained precisely when the previous two conditions are met. The Lyapunov functional expressed in terms of the conformal variables reads

$\displaystyle L(\tau,\gamma,\kappa^{tr},E_{\gamma})=\int_{M}(\tau^{2}-\frac{2n\Lambda}{n-1})^{n/2}\varphi^{2n/(n-2)}(\tau,\gamma,\kappa^{tr},E_{\gamma})\mu_{\gamma}.$

(59)

Utilizing the lower bound (57), the infimum of the Lyapunov functional over the space $(\mathcal{M}_{-1}/D_{0})\times S^{tr}_{2}(M)\times\mathcal{E}$ ($\mathcal{E}$ is the space of re-scaled energy densities, $S^{tr}_{2}(M)$ is the space of symmetric covariant traceless 2-tensors with respect to the metric $\gamma$) may be computed as

	$\displaystyle\inf_{\mathcal{M}_{-1}/D_{0}\times S^{tr}_{2}(M)\times\mathcal{E}}L(\tau,\gamma,\kappa^{tr},E_{\gamma})$		(60)
	$\displaystyle=\inf_{\mathcal{M}_{-1}/D_{0}}\int_{M}(\tau^{2}-\frac{2n\Lambda}{n-1})^{n/2}\varphi^{2n/(n-2)}(\tau,\gamma,\kappa^{tr}=0,E_{\gamma}=0)\mu_{\gamma},$
	$\displaystyle=(\frac{n}{n-1})^{n/2}\inf_{\mathcal{M}_{-1}/D_{0}}\int_{M}\mu_{\gamma},$
	$\displaystyle=(\frac{n}{n-1})^{n/2}\left\{-\sigma(M)\right\}^{n/2},$

where $\sigma(M)$($<0$) is a topological invariant (the higher dimensional analog of the Euler characteristic of a higher genus surface) of the manifold $M$ (of negative Yamabe type considered here).

The most interesting case here is that of the physical universe (i.e., the $3+1$ dimensional case). Utilizing Ricci-flow techniques, the $\sigma$ constant (and therefore the infimum of the weak Lyapunov functional) of the most general compact 3-manifold of negative Yamabe type has been computed and as such is given by

$\displaystyle|\sigma(M)|$

$\displaystyle=$

$\displaystyle(vol_{-1}H)^{2/3},$

(61)

where $vol_{-1}H$ is the volume of the hyperbolic part of $M$ computed with respect to the hyperbolic metric normalized to have scalar curvature $-1$ [33, 34]. Therefore, apart from the hyperbolic family of $\mathbf{K}(\pi,1)$ manifolds, the remaining components of $M$ i.e., wormholes ($\mathbb{S}^{2}\times\mathbb{S}^{1}$), spherical space forms ($\mathbb{S}^{3}/\Gamma_{k}$), and the graph manifolds (non-hyperbolic parts of $\mathbf{K}(\pi,1)$ manifolds) do not contribute to the $\sigma$ constant. Hence, they do not contribute to the infimum of the weak Lyapunov functional as well. Since the weak Lyapunov functional is geometrically the rescaled volume of $M$ (25), following its monotonic decay towards an infimum dominated only by the hyperbolic components of the spatial manifold, one is led to the natural conclusion that the Einstein flow, in the presence of suitable matter sources and a positive cosmological constant, drives the universe towards an asymptotic state that is volume dominated by hyperbolic components. In other words, the dynamical mechanism at work within the Einstein flow (with or without matter) drives the physical universe towards an asymptotic state which is characterized by a locally homogeneous and isotropic spatial metric.

In order to understand the extent to which the weak Lyapunov functional decays to its infimum (or may be obstructed from doing so), it is necessary to study the global properties of the solutions. This is indeed a very difficult open problem of classical general relativity. The recent progress on this so called ‘large data’ problem is only limited to a few highly symmetric cases (e.g., Gowdy spacetimes [28, 29, 31]) or in lower dimensions (e.g., ‘2+1’ vacuum gravity [32] where special techniques have been implemented). The Lyapunov functional (25) may always be used to study large data problems. However, as we have mentioned previously, it is not sufficient to control the desired norm of the data ($g,k$ and matter degrees of freedom) where a local existence theorem may be expected to hold. A solution of the Einstein-matter system may break down via curvature concentration or breakdown of the matter field itself before the Lyapunov functional achieves its infimum. In order to investigate the dynamical behaviour of the matter coupled Einstein equations and conclude to what extent the Lyapunov functional decays, it is necessary to study the stability properties of the solutions on which the later achieves its infimum and the matter source vanishes.

Before tackling this extremely difficult ‘large’ data problem (which is not even resolved in the case of pure vacuum gravity let alone allowing for an arbitrary matter source), it is natural to ask whether solutions sufficiently close to the solutions ((34) or (35)) do indeed remain sufficiently close or converge to these backgrounds. First we consider the pure vacuum limit. In the pure vacuum case, [24] studied the stability of the $3+1$ dimensional Lorentz cone spacetimes (35) utilizing the Bel-Robinson energy and its higher order generalizations (which controlled the required $H^{3}\times H^{2}$ norm of the data $(g,k)$). Later these authors [36] generalized the results for the case of $n\geq 3$ utilizing a wave equation type energy and (its higher order generalization) which played the role of a Lyapunov functional for sufficiently small data. This indeed proved that the non-isolated fixed points playing the role of a center manifold of the suitably re-scaled Einstein dynamics on manifolds that do admit negative Einstein metrics are actual attractors for sufficiently small perturbations. Following these earlier studies several additional ones have been performed addressing the stability issue on manifolds admitting negative spatial Einstein metrics (hyperbolic manifolds for the $3+1$ dimensional case) including source terms in the case of small data. These include the $\Lambda-$vacuum ($\Lambda>0$) [38, 37], Klein-Gordon field [40, 39] (with $\Lambda=0$), dust [30], and Vlasov matter [44] cases, in particular. The stability of these special solutions incorporating several matter sources that satisfy the energy condition we prescribe for the analysis to hold true provides some ‘weak’ support for our conclusion.

Here we will study the stability properties of a special type of background solution that is isometric to the solution given in the earlier section ((34) or (35)). This is done simply because the Milne model (or its equivalent ‘conformal’ spacetime when $\Lambda>0$) is not compatible with physical observation since it is devoid of any matter content. The solutions that we consider are variants of FLRW solutions where the spatial component is given by a compact negative Einstein space (hyperbolic for the $3+1$ dimensional case). These models are the special solutions of the Einstein field equations coupled to a perfect fluid matter source and are explicitly given as follows

	${}^{n+1}g$	$\displaystyle=$	$\displaystyle-dt\otimes dt+a^{2}(t)\gamma_{ij}dx^{i}\otimes dx^{j},t\in[0,\infty)$		(62)
	$\displaystyle R[\gamma]_{ij}$	$\displaystyle=$	$\displaystyle-\frac{1}{n}\gamma_{ij},a\sim t~{}~{}as~{}~{}t\to\infty~{}~{}for~{}~{}\Lambda=0,$		(63)
	$\displaystyle a(t)$	$\displaystyle\sim$	$\displaystyle\sinh(\alpha t)\sqrt{a^{2}(0)+\frac{1}{2\Lambda}}+a(0)\cosh(\alpha t)~{}~{}as~{}~{}t\to\infty,~{}~{}\Lambda>0,$		(64)
	$\displaystyle a(0)$	$\displaystyle:=$	$\displaystyle a(t=0),\alpha:=\sqrt{\frac{2\Lambda}{n(n-1)}}.$

In addition, the density and the fluid $n-$velocity satisfy $\rho(t)=\frac{C}{a^{n\gamma_{a}}(t)},v^{i}=0$ where a barotropic equation of state $P=(\gamma_{a}-1)\rho$, $\gamma_{a}\in(1,2)$, is chosen and $C$ is a constant. One may study the global properties of the sufficiently small perturbations of the suitably re-scaled Einstein-Euler field equations about the background solutions described above. The re-scaled Einstein-Euler equations are generically a non-autonomous dynamical system since the scale factor $a(t)$ appears explicitly in the field equations. This is a vital difference from the pure vacuum case where after a suitable re-scaling by powers of the mean extrinsic curvature $\tau$, one may reduce the field equations to an autonomous form. This lack of autonomous character does not necessarily cause a problem provided that the scale factor satisfies a suitable integrability condition (which holds trivially for the current model).