# A dynamical systems formulation for inhomogeneous LRS-II spacetimes

Saikat Chakraborty  0000-0002-5472-304X Peter K.S. Dunsby  0000-0002-6271-9585
Rituparno Goswami  0000-0002-5355-3288 Amare Abebe  0000-0001-5475-2919

###### Abstract

We present a dynamical system formulation for inhomogeneous LRS-II spacetimes
using the covariant 1+1+2 decomposition approach. Our approach describes the
LRS-II dynamics from the point of view of a comoving observer. Promoting the
covariant radial derivatives of the covariant dynamical quantities to new
dynamical variables and utilizing the commutation relation between the
covariant temporal and radial derivatives, we were able to construct an
autonomous system of first-order ordinary differential equations along with
some purely algebraic constraints. Using our dynamical system formulation we
found several interesting features in the LRS-II phase space with dust, one of
them being that the homogeneous solutions constitute an invariant submanifold.
For the particular case of LTB, we were also able to recover the previously
known result that an expanding LTB tends to Milne in the absence of a
cosmological constant, providing a potential validation of our formalism.

##  1 Introduction

One of the key pillars of the standard cosmological model is the assumption of
the cosmological principle (sometimes called the Copernican principle).
According to this principle, averaged over a cosmologically relevant scale,
say a few percent of the particle horizon, it is a good enough approximation
to assume that matter in the universe is distributed in a smooth spatially
homogeneous and isotropic manner. It is further assumed that the spacetime
geometrical features trace that of the matter distribution. Together, these
two assumptions give rise to the modelling of the cosmological spacetime with
the Friedmann-Lemaítre-Robertson-Walker (FLRW) metric. The standard model of
cosmology, the so-called $\Lambda$CDM model is based on the FLRW spacetime
geometry. Small-scale inhomogeneities are treated as cosmological
perturbations [1].

Nevertheless, as is well known, the $\Lambda$CDM model suffers from the
cosmological constant problem [2], $H_{0}$-problem [3] and $\sigma_{8}$
problem [4]. These problems fuelled the search for alternative models to
describe the observable universe. Most attempt to find alternative models to
describe a dynamical dark energy, as opposed to a “constant” dark energy given
by the cosmological constant. These alternative models are constructed by
adding one or more additional dynamical degrees of freedom to the known
theory, either by adding additional matter components (canonical or non-
canonical fields, non-ideal fluids or combinations thereof [5, 6]) or altering
the theory of gravity altogether (scalar-tensor theories dominating this
direction [7, 8]). Of course, there are clear problems associated with both
approaches. If one adds additional matter components, one needs to worry about
their explanation from the particle physics side. On the other hand, modified
gravity theories are usually highly constrained by local gravity tests and
many a time are plagued by various stability issues. For an extensive review
of the dynamical analysis of such alternative models of dark energy, the
reader is referred to the beautiful review by Bahamonde et. al. [9].

An interesting alternative line of explaining the apparent acceleration of the
universe without the introduction of any new dynamical degree of freedom to
the theory is based on relaxing the cosmological principle, in particular the
assumption of homogeneity. Since a full non-perturbative treatment of generic
inhomogeneous spacetime dynamics is very difficult, the toy model of
spherically symmetric inhomogeneous spacetimes, namely the Lemaítre-Tolman-
Bondi (LTB) model was considered. One class of works have studied LTB models
with a purely formal notion of acceleration [10, 11, 12, 13]. Another class of
works focused on whether it is possible to explain without assuming any notion
of an acceleration the observed relationship between the luminosity distance
of redshift [14, 15, 16, 17], the very relationship which gave rise to the
notion of an accelerated universe under the assumption of homogeneity. Of
course, even these models come with their fair share of observational
constraints [18, 19, 20]. Nonetheless, two explicit theorems were proved by
Mustapha, Hellaby and Ellis in [21] which make it clear that homogeneity
cannot be proven without either a fully determinate theory of source evolution
or availability of distance measures that are independent of the source
evolution. This theorem makes it worthwhile to further study inhomogeneous
cosmology. Inhomogeneous LTB models are also the starting point of studies on
realistic gravitational collapse [22, 23, 24].

Dynamical system analysis has proved to be an invaluable tool for studying the
dynamics of homogeneous spacetimes. Since in standard cosmology we mostly deal
with homogeneous spacetimes, the dynamical system approach in cosmology is an
invaluable companion to theoretical cosmologists [25, 26]. For homogeneous
dynamical spacetimes, this formalism recasts the gravitational field equations
as a system of autonomous first-order ordinary differential equations along
with some algebraic constraints. The gravitational dynamics are then recast as
a flow in a phase space spanned by suitably defined dimensionless dynamical
variables (which are usually normalized by the expansion parameter of a
timeline congruence). When a self-gravitating system is inhomogeneous, as in
an LTB cosmology or LTB collapse model, the techniques of the dynamical system
cannot be applied in a straightforward manner. The reason is that now the
governing equations of the system are partial differential equations.

That is not to say that such attempts have not been made. One such approach
that could be found in the literature utilizes the so-called orthonormal frame
formalism [27]. The issue with this approach is that both the evolution and
constraint equations are partial differential equations, involving temporal
and spatial derivatives respectively. Although some techniques from dynamical
systems theory could be applied, the analysis becomes mathematically intricate
[28, 29, 30]. Another approach that is taken by Sussman introduces the so-
called quasi-local variables [31, 32]. Sussman’s approach does indeed give
rise to an effective system of ordinary differential equations with algebraic
constraints, but the quasilocal variables by themselves do not directly
correspond to the actual covariant geometric and thermodynamic quantities
relevant for a comoving observer.

In this work, we present a new dynamical system formulation for a class of
inhomogeneous spacetimes with a local rotational symmetry, namely the LRS-II.
The motivation to consider this class is that it contains the spherically
symmetric solutions e.g. LTB that are of physical interest. Our formulation is
based on the 1+1+2 covariant decomposition approach [33, 34], which itself is
built up on the 1+3 covariant decomposition approach [35, 36] 111That 1+1+2
covariant approach can potentially be helpful in constructing a dynamical
system formulation for LRS-II spacetimes is hinted in [37, 38]. However, the
authors considered only LRS-II vacuum or electrovacuum solutions for which
case it reduces to either static or homogeneous spacetime by virtue of the
Birkhoff’s theorem.. To proceed with this approach one needs to first specify
a timeline unit vector $u^{\mu}$ and a spacelike unit vector $e^{\mu}$, which
we choose to be along the worldline of a comoving observer and the direction
of the axis of local rotational symmetry respectively. In the covariant
decomposition approaches, all the dynamical quantities have a very clear and
coordinate-independent physical interpretation. Physically, our approach
involves going into the frame of a comoving observer and following the
dynamics from his/her point of view. The effect of spatial inhomogeneity is
that not all the comoving observers would experience the same evolutionary
phase (characterized by a fixed point), simultaneously, a point we also try to
elaborate on in the main text. The novelty of our approach is that we could
formulate an autonomous system with only time-dependent evolution equations
and algebraic constraints and our dynamical variables are also directly
related to covariant geometrical and thermodynamic quantities relevant to a
comoving observer.

The paper is organized as follows. In section 2 we present an exposition to
LRS-II spacetimes using the 1+1+2 covariant formalism. We dedicate a complete
subsection discussing various intricacies presented by the inhomogeneity of
the system. In section 3 we construct an autonomous dynamical system for LRS-
II using expansion normalized dimensionless dynamical variables. We dedicate
an entire subsection to discussing the interpretation of the dynamical
structure, something which we deem important here. Sections 4, 5 and 6 present
the invariant submanifolds, finite fixed points and asymptotic fixed points in
the LRS-II phase space. In section 7 we specialize to the spherically
symmetric LTB model and discuss its phase space. We conclude in section 8
summarizing the important and novel aspects of our work with an outlook of
possible further research along this line. Finally, appendix A is devoted to
exploring a possible different route to constructing an autonomous system for
cases with positive 3-curvatures, which proves to be not very fruitful.

##  2 Dynamics of LRS-II spacetime with a perfect fluid

We consider the LRS-II class of spacetimes that admits spherically symmetric
solutions and is rotation-free. The most general LRS-II metric can be written
as [39]

$ds^{2}=-\frac{1}{A^{2}(t,r)}dt^{2}+B^{2}(t,r)dr^{2}+C^{2}(t,r)[dy^{2}+D_{k}^{2}(y)dz^{2}]\,,$
(2.1)

where $t,\,r$ are affine parameters along $u^{\mu}$ and $e^{\mu}$ and
$k=(-1,0,1)$ denotes open, flat and closed geometry of the 2-sheets
respectively. The function $D_{k}(y)$ is

$D_{k}(y)=\begin{cases}\sin{y},\qquad k=+1\\\ y,\qquad k=0\\\ \sinh{y},\qquad
k=-1\end{cases}.$ (2.2)

Over the last few years, a new approach, originally introduced by Clarkson
[33] has been used extensively in the study of LRS-II spacetimes, both in
General Relativity and in $f(R)$ gravity [40, 41, 42, 43] .

The 1+1+2 covariant approach builds on the well-known 1+3 approach to
cosmology pioneered by Ehlers [35] and Ellis [36] by further splitting the
three spatial degrees of freedom into a part that lies along a spatial vector
$e^{a}$ and one perpendicular to $e^{a}$ (referred to as the “sheet”). In this
way any spatial tensor can be split into a part along $e^{a}$ and one
orthogonal to $e^{a}$.

In general, there is a set of evolution equations containing only
‘dot’-derivatives, a set of propagation equations containing only
‘hat’-derivatives and a set of mixed evolution-propagation equations
containing both hat and dot derivatives. The amazing simplicity offered by
confining our attention only to perfect fluids is that there we can get rid of
the mixed equations. The evolution and the propagation equations can be
completely decoupled.

The complete set of 1+1+2 covariant equations for LRS-II spacetimes in the
presence of a generic fluid can be found in [34]. Here we only give the
definitions of the key quantities related to this problem, together with their
evolution and constraint equations.

The set of non-vanishing covariant geometric and thermodynamic quantities
characterizing an LRS-II spacetime with a perfect fluid at the background
level are listed below along with their physical interpretations.

  * •

Geometric quantities:

    * –

$\theta\equiv D_{\mu}u^{\mu}$; the expansion of the congruence with the
timelike vector $u^{\mu}$.

    * –

$\phi\equiv d_{\mu}e^{\mu}$; the expansion of the 2-sheet orthogonal to
$u^{\mu},\,e^{\mu}$ along the preferred spatial direction specified by the
local rotational symmetry.

    * –

$\Sigma\equiv
e^{\mu}e^{\nu}\sigma_{\mu\nu}=e^{\mu}e^{\nu}D_{\langle\mu}u_{\nu\rangle}=e^{\mu}e^{\nu}\left(h^{\alpha}_{\,\,\,(\mu}h^{\beta}_{\,\,\,\nu)}-\frac{1}{3}h_{\mu\nu}h^{\alpha\beta}\right)D_{\alpha}u_{\beta}$;
projection of the shear tensor $\sigma_{\mu\nu}$ along the spacelike direction
$e^{\mu}$.

    * –

$\mathcal{E}\equiv
E_{\mu\nu}e^{\mu}e^{\nu}=C_{\mu\alpha\nu\beta}e^{\mu}e^{\nu}u^{\alpha}u^{\beta}$;
projection of the electric part $E_{\mu\nu}$ of the Weyl curvature tensor
$C_{\alpha\beta\gamma\delta}$ along the spacelike direction $e^{\mu}$.

    * –

${}^{3}R$; the 3-curvature of the 3-space orthogonal to $u^{\mu}$.

  * •

Thermodynamic quantities:

  * •

$\mu\equiv T_{\mu\nu}u^{\mu}u^{\nu}$; energy density of the perfect fluid as
measured locally by an observer with the 4-velocity $u^{\mu}$, $T_{\mu\nu}$
being the energy-momentum tensor.

  * •

$P\equiv\frac{1}{3}T_{\mu\nu}h^{\mu\nu}$; pressure of the perfect fluid.

For elaborate discussions on how these quantities arise in a
1+1+2-decomposition set-up, the reader is referred to the references [33, 34].
When the acceleration vector $a^{\mu}=u^{\nu}\nabla_{\nu}u^{\mu}$ for the
timelike congruence vanishes, as in the present case under consideration, the
4-vector $u^{\mu}$ satisfies the geodesic equation. Since the time coordinate
$t$ in the metric (2.1) is chosen as the affine parameter along $u^{\mu}$, we
can identify the timelike congruence with the congruence of comoving
geodesics. 222The explicit form of $u^{\mu},\,e^{\mu}$ for the metric (2.1)
are $u^{\mu}=(A(t,r),0,0,0),\,\,e^{\mu}=\left(0,\frac{1}{B(t,r)},0,0\right)$.

The evolution and propagation equations for LRS-II spacetimes with a perfect
barotropic fluid ($P=P(\mu)$) are as follows

  * •

Evolution equations:

$\displaystyle\dot{\phi}=-\frac{1}{2}\phi\left(\frac{2}{3}\theta-\Sigma\right),$
(2.3a)
$\displaystyle\dot{\theta}=-\frac{1}{3}\theta^{2}-\frac{3}{2}\Sigma^{2}-\frac{1}{2}(\mu+3P(\mu)),$
(2.3b) $\displaystyle\dot{\mu}=-\theta(\mu+P(\mu)),$ (2.3c)
$\displaystyle\dot{\Sigma}=-\frac{1}{2}\Sigma^{2}-\frac{2}{3}\theta\Sigma-\mathcal{E},$
(2.3d)
$\displaystyle\dot{\mathcal{E}}=-\frac{3}{2}\mathcal{E}\left(\frac{2}{3}\theta-\Sigma\right)-\frac{1}{2}(\mu+P(\mu))\Sigma.$
(2.3e)

  * •

Propagation equations:

$\displaystyle\hat{\phi}=-\frac{1}{2}\phi^{2}+\left(\frac{1}{3}\theta+\Sigma\right)\left(\frac{2}{3}\theta-\Sigma\right)-\frac{2}{3}\mu-\mathcal{E},$
(2.4a)
$\displaystyle\hat{\Sigma}-\frac{2}{3}\hat{\theta}=-\frac{3}{2}\phi\Sigma,$
(2.4b)
$\displaystyle\hat{\mathcal{E}}-\frac{1}{3}\hat{\mu}=-\frac{3}{2}\phi\mathcal{E},$
(2.4c) $\displaystyle\hat{P}=P,_{\mu}\hat{\mu}=0.$ (2.4d)

  * •

Friedmann constraint:

$\frac{\theta^{2}}{9}+\frac{{}^{3}R}{6}=\frac{\mu}{3}+\frac{\Sigma^{2}}{4}.$
(2.5)

The Gaussian curvature $K$ of the 2D sheet is given by
${}^{2}R_{\mu\nu}=KN_{\mu\nu}$. The following relation holds [40, 44] between
the Gaussian curvature and the 3-curvature

${}^{3}R=-2\left(\frac{1}{2}\hat{\phi}+\frac{3}{4}\phi^{2}-K\right)\,.$ (2.6)

Using the $\hat{\phi}$-equation in (2.4), one can get the following relation

$K=\frac{\mu}{3}-\mathcal{E}+\frac{1}{4}\phi^{2}-\left(\frac{\theta}{3}-\frac{\Sigma}{2}\right)^{2}\,.$
(2.7)

Its evolution and propagation equations are

$\displaystyle\dot{K}=-\frac{2}{3}\left(\frac{2}{3}\theta-\Sigma\right)K\,,$
(2.8a) $\displaystyle\hat{K}=-\phi K\,.$ (2.8b)

The covariant time and radial derivatives usually do not commute. For a scalar
quantity $\psi$, the following commutation relation holds [45]

$\hat{\dot{\psi}}-\dot{\hat{\psi}}=\left(\frac{1}{3}\theta+\Sigma\right)\hat{\psi}\,.$
(2.9)

An important observation is that, in the absence of the fluid momentum density
or anisotropic pressure, the last of the propagation equations implies either
a homogeneous matter distribution at all times ($\hat{\mu}=0$) or a fluid with
constant pressure. Henceforth we will consider a dust fluid for which the
pressure identically vanishes and therefore the last propagation equation is
identically satisfied.

The propagation equations (2.4) are spacelike constraints that must hold at
all the constant time hypersurfaces. In other words, these constraints must be
conserved in time. To show this explicitly, let us define the following
combinations:

$\displaystyle\mathcal{C}_{1}\equiv\hat{\phi}+\frac{1}{2}\phi^{2}-\left(\frac{1}{3}\theta+\Sigma\right)\left(\frac{2}{3}\theta-\Sigma\right)+\frac{2}{3}\mu+\mathcal{E},$
(2.10a)
$\displaystyle\mathcal{C}_{2}\equiv\hat{\Sigma}-\frac{2}{3}\hat{\theta}+\frac{3}{2}\phi\Sigma,$
(2.10b)
$\displaystyle\mathcal{C}_{3}\equiv\hat{\mathcal{E}}-\frac{1}{3}\hat{\mu}+\frac{3}{2}\phi\mathcal{E}\;.$
(2.10c)

The propagation equations can be represented as

$\\{\mathcal{C}_{1},\mathcal{C}_{2},\mathcal{C}_{3}\\}=\\{0,0,0\\}.$ (2.11)

Using the commutation relation (2.9) we can derive the following equations

$\displaystyle\dot{\mathcal{C}}_{1}=-\left(\frac{2}{3}\theta+\frac{1}{2}\Sigma\right)\mathcal{C}_{1}+\frac{1}{2}\phi\mathcal{C}_{2},$
(2.12a)
$\displaystyle\dot{\mathcal{C}}_{2}\equiv-\theta\mathcal{C}_{2}-\mathcal{C}_{3},$
(2.12b)
$\displaystyle\dot{\mathcal{C}}_{3}=-\frac{1}{2}(\mu+P-3\mathcal{E})\mathcal{C}_{2}-\left(\frac{4}{3}\theta-\frac{1}{2}\Sigma\right)\mathcal{C}_{3}\;.$
(2.12c)

Eqs.(2.12) essentially dictate the evolution of the spacelike constraints
given by Eq.(2.11). We see that the constraint (2.11) is conserved in time,
i.e. if it is satisfied initially, it is satisfied at all times.

The spacelike constraints in Eq.(2.4) are not purely algebraic as they involve
radial derivatives. However, using the commutation relation in Eq.(2.9) it is
possible to recast the governing equations as an autonomous system of first-
order differential equations with the constraints now being interpreted as
purely algebraic ones. The idea is to promote the hat derivatives to separate
dynamical variables and calculate their time evolution using the commutation
relation Eq.(2.9). The propagation equations in Eq.(2.4) can then be
interpreted as purely algebraic constraints. Hence we get an extended
dynamical system (considering dust fluid) given by the following dynamical
equations.

  * •

$\displaystyle\dot{\phi}=-\frac{1}{2}\phi\left(\frac{2}{3}\theta-\Sigma\right),$
(2.13a)
$\displaystyle\dot{\theta}=-\frac{1}{3}\theta^{2}-\frac{3}{2}\Sigma^{2}-\frac{1}{2}\mu,$
(2.13b) $\displaystyle\dot{\mu}=-\theta\mu,$ (2.13c)
$\displaystyle\dot{\Sigma}=-\frac{1}{2}\Sigma^{2}-\frac{2}{3}\theta\Sigma-\mathcal{E},$
(2.13d)
$\displaystyle\dot{\mathcal{E}}=-\frac{3}{2}\mathcal{E}\left(\frac{2}{3}\theta-\Sigma\right)-\frac{1}{2}\mu\Sigma,$
(2.13e)
$\displaystyle\dot{\hat{\phi}}=-\left(\frac{2}{3}\theta+\frac{1}{2}\Sigma\right)\hat{\phi}-\frac{1}{2}\phi\left(\frac{2}{3}\hat{\theta}-\hat{\Sigma}\right),$
(2.13f)
$\displaystyle\dot{\hat{\theta}}=-\theta\hat{\theta}-3\Sigma\hat{\Sigma}-\frac{1}{2}\hat{\mu}-\Sigma\hat{\theta},$
(2.13g)
$\displaystyle\dot{\hat{\mu}}=-\left(\frac{4}{3}\theta+\Sigma\right)\hat{\mu}-\mu\hat{\theta},$
(2.13h)
$\displaystyle\dot{\hat{\Sigma}}=-(\theta+2\Sigma)\hat{\Sigma}-\frac{2}{3}\Sigma\hat{\theta}-\hat{\mathcal{E}},$
(2.13i)
$\displaystyle\dot{\hat{\mathcal{E}}}=-\left(\frac{4}{3}\theta-\frac{1}{2}\Sigma\right)\hat{\mathcal{E}}-\frac{3}{2}\mathcal{E}\left(\frac{2}{3}\hat{\theta}-\hat{\Sigma}\right)-\frac{1}{2}\hat{\mu}\Sigma-\frac{1}{2}\mu\hat{\Sigma}\;,$
(2.13j)

along with the set of algebraic constraints given by (2.4). There are three
such algebraic constraints in the case of dust matter, allowing us to
eliminate $\hat{\phi},\,\hat{\Sigma},\,\hat{\mathcal{E}}$. The reduced
dynamical system is:

$\displaystyle\dot{\phi}=-\frac{1}{2}\phi\left(\frac{2}{3}\theta-\Sigma\right),$
(2.14a)
$\displaystyle\dot{\theta}=-\frac{1}{3}\theta^{2}-\frac{3}{2}\Sigma^{2}-\frac{1}{2}\mu,$
(2.14b) $\displaystyle\dot{\mu}=-\theta\mu,$ (2.14c)
$\displaystyle\dot{\Sigma}=-\frac{1}{2}\Sigma^{2}-\frac{2}{3}\theta\Sigma-\mathcal{E},$
(2.14d)
$\displaystyle\dot{\mathcal{E}}=-\frac{3}{2}\mathcal{E}\left(\frac{2}{3}\theta-\Sigma\right)-\frac{1}{2}\mu\Sigma,$
(2.14e)
$\displaystyle\dot{\hat{\theta}}=-\theta\hat{\theta}-3\Sigma\hat{\theta}+\frac{9}{2}\Sigma^{2}\phi-\frac{1}{2}\hat{\mu},$
(2.14f)
$\displaystyle\dot{\hat{\mu}}=-\left(\frac{4}{3}\theta+\Sigma\right)\hat{\mu}-\mu\hat{\theta}.$
(2.14g)

###  2.1 Some discussion about the structure of the system

It is worthwhile here to ponder upon the structure of the system of equations
(2.14). The system (2.14) by itself is a closed set, but one can notice a nice
decoupling here. The
$\\{\dot{\theta},\dot{\mu},\dot{\Sigma},\dot{\mathcal{E}}\\}$-equations are
decoupled from the other equations. $\theta,\,\mu,\,\Sigma,\,\mathcal{E}$ by
themselves form a closed set of equations since they do not contain the other
covariant quantities. This system is well-posed [46]. Once one knows
$\theta,\,\mu,\,\Sigma,\,\mathcal{E}$, either in the form of fixed points or
as numerical solutions, one can plug them into the
$\\{\dot{\phi},\dot{\hat{\theta}},\dot{\hat{\mu}}\\}$-equations and solve for
$\phi(t),\,\hat{\theta}(t),\,\hat{\mu}(t)$. The propagation equations (2.4)
then allow us to find out $\hat{\phi}$, $\hat{\Sigma}$ and
$\hat{\mathcal{E}}$.

We are considering here a congruence of comoving geodesics. Let us follow the
geodesic of a particular comoving observer and consider a covariant quantity
$Q$. If we go to the tangent space of a comoving observer, because of the
inhomogeneity of the system, the observer has at his disposal two derivatives;
$\dot{Q}$ and $\hat{Q}$. The $\dot{Q}$ gives the change in the value of $Q$ at
his tangent space, i.e. along his particular geodesic, while $\hat{Q}$ can be
interpreted as a “radial perturbation” of $Q$ on the same constant time
hypersurface, where time is measured along the unique timelike 4-vector field
$u^{\mu}$ 333This is not a generic spatial perturbation, though, but a
specific kind of perturbations respecting the local rotational symmetry.. This
can be seen as follows. Consider a comoving observer at $r=r_{0}$, who, at a
particular time $t_{0}$, has the value $Q=Q_{0}$. Now consider a neighboring
comoving observer at $r=r_{0}+\delta r$. Since both the observers are comoving
observers, their geodesics belong to the same congruence. The deviation vector
between their geodesics is given by $\xi^{\mu}=\delta re^{\mu}$. The value of
at the time $t_{0}$ for this neighbouring observer is given by

$Q=Q_{0}+\delta_{\xi}Q=Q_{0}+\xi^{\mu}\nabla_{\mu}Q=Q_{0}+\delta
re^{\mu}(-\dot{Q}u_{\mu}+\hat{Q}e_{\mu}+N^{\,\,\,\nu}_{\mu}D_{\nu}Q)=Q_{0}+\hat{Q}\delta
r\,,$ (2.15)

where in the last step we utilized the orthogonality of $e^{\mu}$ with
$u^{\mu}$ and $N_{\mu\nu}$. With respect to the comoving observer at
$r=r_{0}$, the non-vanishing $\delta_{\xi}Q=\hat{Q}\delta r$ is precisely what
is stopping the other comoving observers from having the same $Q$-value
$Q_{0}$, making the spatial distribution of $Q$ inhomogeneous. This justifies
the interpretation of $\hat{Q}$ as a perturbation. The
$\dot{\hat{Q}}$-equation gives the dynamics of this radial perturbation.

$\theta,\,\Sigma,\,\mathcal{E},\,\mu$ are the expansion, shear, electric part
of the Weyl tensor and the energy density locally measured by a comoving
observer. On the other hand, the covariant quantity $\phi$ is related to the
Gaussian curvature of a shell (Eq.(2.7)), which characterizes the “shape” of
the shells. The sign of the Gaussian curvature distinguishes between
spherically, cylindrically and hyperbolically symmetric spatial geometry. For
spherically and cylindrically symmetric situations, different comoving shells
will in general have different values of the Gaussian curvature. When this is
the case, $\phi$ can be thought of as labels for different shells.

The decoupling of the
$\\{\dot{\theta},\,\dot{\Sigma},\,\dot{\mathcal{E}},\,\dot{\mu}\\}$-equations
from the other equations imply that the qualitative nature of the evolution of
expansion, shear, electric part of the Weyl tensor and the energy density is
the same for all the comoving observers and is independent of how the radial
perturbations in the respective quantities evolve. However, the existence of
such radial perturbations implies that the evolution of different shells will
fall “out of phase”. Different shells will not experience the same
evolutionary phase simultaneously. This is the key contrast with a homogeneous
geometry e.g. FLRW, where _all_ the shells will experience the same
evolutionary phase simultaneously.

In general, when the fluid is not perfect, there will be mixed equations and
the radial perturbations do affect the shell dynamics. In such cases,
$\\{\dot{\theta},\dot{\mu},\dot{\Sigma},\dot{\mathcal{E}}\\}$-equations cannot
be decoupled from the other equations.

It should be clear from the above discussion that if $Q$ is a covariant
quantity, then in our approach $Q=0$ does _not_ automatically imply
$\hat{Q}=0$. Each comoving observer has a 4D tangent space. Consider locally
the 2D projection of the tangent space of a comoving observer spanned by the
4-vectors $u^{\mu}$ and $e^{\mu}$ (the local $u$-$e$ plane). In the local
$u$-$e$ plane of the comoving observer, $Q=0$ corresponds to a curve. The
normal vector to this curve is given by

$\nabla_{\mu}Q=-\dot{Q}u_{\mu}+\hat{Q}e_{\mu}+N^{\nu}_{\mu}D_{\nu}Q\,.$ (2.16)

The last term, $N^{\nu}_{\,\,\,\mu}D_{\nu}Q$ is the gradient of $Q$ on the
2-sheet orthogonal to $u^{\mu}$ and $e^{\mu}$. Since in the 1+1+2
decomposition of LRS spacetimes the spacelike vector field $e^{\mu}$ at each
spatial point is taken to be along the axis of local rotational symmetry, the
last term identically vanishes. Let the tangent to the curve $Q=0$ in the
local $u$-$e$ plane be given by

$Q^{\mu}=\alpha u^{\mu}+\beta e^{\mu}.$ (2.17)

The condition $Q^{\mu}\nabla_{\mu}Q=0$ helps to determine the slope to this
curve in the local $u$-$e$ plane.

$-\frac{\alpha}{\beta}=\frac{\hat{Q}}{\dot{Q}}.$ (2.18)

If this slope is positive, then $Q=0$ is attained by a comoving observer at a
larger radius at a later time and by a comoving observer at a smaller radius
at a smaller time. The opposite is true when the slope is negative. In special
cases when $\hat{Q}$ also vanishes along the $Q=0$ curve, the slope is zero
and all comoving observers attain $Q=0$ at the same time. In this case, the
$Q$ distribution is homogeneous.444Which does not mean that the spacetime
geometry as a whole should be homogeneous. The hat of other covariant
quantities can be non-vanishing.

##  3 A dynamical system formulation for LRS-II spacetime with dust

Throughout the rest of the paper, we specialize to that case of pressureless
dust. The first step towards a dynamical system analysis is to define
dimensionless dynamical variables:

$\displaystyle x_{1}\equiv\frac{\phi}{\theta},\quad
x_{2}\equiv\frac{3\mu}{\theta^{2}},\quad
x_{3}\equiv\frac{3}{2}\frac{\Sigma}{\theta},\quad
x_{4}\equiv\frac{\mathcal{E}}{\theta^{2}},$ (3.1) $\displaystyle
y_{1}\equiv\frac{\hat{\phi}}{\theta^{2}},\quad
y_{2}\equiv\frac{3\hat{\mu}}{\theta^{3}},\quad
y_{3}\equiv\frac{3}{2}\frac{\hat{\Sigma}}{\theta^{2}},\quad
y_{4}\equiv\frac{\hat{\mathcal{E}}}{\theta^{3}},$ $\displaystyle
z\equiv\frac{\hat{\theta}}{\theta^{2}}.$

In terms of these variables the Friedmann constraint (2.5) is

$x_{2}+x_{3}^{2}=1+\frac{3}{2}\frac{{}^{3}R}{\theta^{2}}.$ (3.2)

Next we introduce a dimensionless time variable $\tau$ as
$d\tau=\epsilon\theta dt$, where $\epsilon=+1$ for an expanding LRS-II
spacetime ($\theta>0$) and $\epsilon=-1$ for a contracting LRS-II spacetime
($\theta<0$). We imported the $\epsilon$ so that $\tau$ is always a
monotonically increasing function of time, and therefore can be justifiably
taken as a time variable on the phase space. So the dynamical system
formulation we develop here can be applied to study both inhomogeneous dark
energy models and the phenomenon of gravitational collapse.

An important relationship that is worth writing here is

$\frac{\dot{\theta}}{\theta^{2}}=-\frac{1}{6}(2+x_{2}+4x_{3}^{2}).$ (3.3)

Considering a dust fluid ($P(\mu)=0$) and using the equations (2.13), we get
the following dynamical system

$\displaystyle\frac{dx_{1}}{\epsilon
d\tau}=\frac{1}{6}x_{1}(2x_{3}+x_{2}+4x_{3}^{2}),$ (3.4a)
$\displaystyle\frac{dx_{2}}{\epsilon
d\tau}=-\frac{1}{3}x_{2}(1-x_{2}-4x_{3}^{2}),$ (3.4b)
$\displaystyle\frac{dx_{3}}{\epsilon
d\tau}=-\frac{1}{6}x_{3}(2+2x_{3}-x_{2}-4x_{3}^{2})-\frac{3}{2}x_{4},$ (3.4c)
$\displaystyle\frac{dx_{4}}{\epsilon
d\tau}=-\frac{1}{3}x_{4}(1-3x_{3}-x_{2}-4x_{3}^{2})-\frac{1}{9}x_{2}x_{3},$
(3.4d) $\displaystyle\frac{dy_{1}}{\epsilon
d\tau}=\frac{1}{3}y_{1}(x_{2}-x_{3}+4x_{3}^{2})+\frac{1}{3}x_{1}(y_{3}-z),$
(3.4e) $\displaystyle\frac{dy_{2}}{\epsilon
d\tau}=-\frac{1}{3}y_{2}+\frac{1}{6}y_{2}(3x_{2}-4x_{3})-x_{2}z+2x_{3}^{2}y_{2},$
(3.4f) $\displaystyle\frac{dy_{3}}{\epsilon
d\tau}=-\frac{1}{3}y_{3}-\frac{3}{2}y_{4}+\frac{1}{3}y_{3}(x_{2}-4x_{3})-\frac{2}{3}x_{3}z+\frac{4}{3}x_{3}^{2}y_{3},$
(3.4g) $\displaystyle\frac{dy_{4}}{\epsilon
d\tau}=-\frac{1}{3}y_{4}-\frac{1}{18}x_{2}(2y_{3}-9y_{4})-\frac{1}{9}x_{3}(y_{2}-3y_{4})+x_{4}(y_{3}-z)+2x_{3}^{2}y_{4},$
(3.4h) $\displaystyle\frac{dz}{\epsilon
d\tau}=-\frac{1}{6}y_{2}-\frac{1}{3}z+\frac{1}{3}x_{2}z-\frac{2}{3}x_{3}(2y_{3}+z)+\frac{4}{3}x_{3}^{2}z.$
(3.4i)

The propagation equations (2.4) can now be written explicitly as algebraic
constraints

$\displaystyle
y_{1}=-\frac{1}{2}x_{1}^{2}+\frac{2}{9}(1-x_{3})(1+2x_{3})-\frac{2}{9}x_{2}-x_{4},$
(3.5a) $\displaystyle y_{3}-z=-\frac{3}{2}x_{1}x_{3},$ (3.5b) $\displaystyle
y_{4}-\frac{1}{9}y_{2}=-\frac{3}{2}x_{1}x_{4}$ (3.5c)

Using the three algebraic constraints in Eq.(3.5) to eliminate
$y_{1},\,y_{3},\,y_{4}$ one gets a reduced dynamical system which corresponds
to the reduced system of equations (2.14)

$\displaystyle\frac{dx_{1}}{\epsilon
d\tau}=\frac{1}{6}x_{1}(2x_{3}+x_{2}+4x_{3}^{2}),$ (3.6a)
$\displaystyle\frac{dx_{2}}{\epsilon
d\tau}=-\frac{1}{3}x_{2}(1-x_{2}-4x_{3}^{2}),$ (3.6b)
$\displaystyle\frac{dx_{3}}{\epsilon
d\tau}=-\frac{1}{6}x_{3}(2+2x_{3}-x_{2}-4x_{3}^{2})-\frac{3}{2}x_{4},$ (3.6c)
$\displaystyle\frac{dx_{4}}{\epsilon
d\tau}=-\frac{1}{3}x_{4}(1-3x_{3}-x_{2}-4x_{3}^{2})-\frac{1}{9}x_{2}x_{3},$
(3.6d) $\displaystyle\frac{dy_{2}}{\epsilon
d\tau}=-\frac{1}{3}y_{2}+\frac{1}{6}y_{2}(3x_{2}-4x_{3})-x_{2}z+2x_{3}^{2}y_{2},$
(3.6e) $\displaystyle\frac{dz}{\epsilon
d\tau}=-\frac{1}{6}y_{2}-\frac{1}{3}z+\frac{1}{3}x_{2}z-2x_{3}z-2x_{1}x_{3}^{2}+\frac{4}{3}x_{3}^{2}z.$
(3.6f)

###  3.1 Interpretation of the dynamical structures

The approach that we are going to take here is the following. Since, in a
generically inhomogeneous situation, the dynamics with respect to observers at
different comoving shells are not equivalent, we will go to the frame of a
particular comoving observer and interpret the dynamical structures, namely
fixed points and invariant submanifolds, from that observer’s point of view.
As explained before, a particular comoving observer has, at his disposal the
covariant geometrical and thermodynamic quantities as well as their radial
gradients as given by the hat derivatives. Both these sets have to be taken
into account to get a full picture of the dynamics. The dynamics are described
in the 6D phase space $\\{x_{1},x_{2},x_{3},x_{4},y_{2},z\\}$ by the phase
flow given by the autonomous system (3.6).

In general, if the comoving observer at $r=r_{0}$ is undergoing an
evolutionary phase characterized by a particular fixed point, a neighbouring
observer at $r=r_{0}+\delta r$ need not necessarily undergo the same
evolutionary phase, because of generally non-vanishing radial perturbations.
This can be seen more clearly as follows. If
$\\{x_{1}^{*},x_{2}^{*},x_{3}^{*},x_{4}^{*},y_{2}^{*},z^{*}\\}$ are the
coordinates of a fixed point, then the fixed point is characterized by the
vanishing of the following covariant quantities:

$\\{\phi-x_{1}^{*}\theta,3\mu-
x_{2}^{*}\theta^{2},3\Sigma-2x_{3}^{*}\theta,\mathcal{E}-x_{4}^{*}\theta^{2},3\hat{\mu}-y_{2}^{*}\theta^{3},\hat{\theta}-z^{*}\theta^{2}\\}=\\{0,0,0,0,0,0\\}\,.$

According to the discussion at the end of section 2.1, the vanishing of a
covariant quantity $Q$ does not automatically imply the vanishing of its
radial perturbation $\hat{Q}$. This gives rise to a non-vanishing slope of the
$Q=0$ curve in the local $u$-$e$ plane of a comoving observer, which means
that a neighbouring comoving observer at $r=r_{0}+\delta r$ will attain the
same fixed point at a slightly earlier or a slightly later time. The dynamics
of the radial perturbations with respect to a comoving observer, in a
covariant manner, are given by the $y$-equations.

For a certain fixed point, if
$\\{y_{1},y_{2},y_{3},y_{4},z\\}=\\{0,0,0,0,0\\}$, then that implies
$\\{\hat{\phi},\hat{\theta},\hat{\Sigma},\hat{\mathcal{E}},\hat{\mu}\\}=\\{0,0,0,0,0\\}$;
the radial perturbations in all the covariant quantities of interest vanishes.
In such a scenario, observers at different comoving shells measure the same
matter density and experience the same spacetime geometry. Such a fixed point
necessarily corresponds to a completely homogeneous scenario. An example of
such a geometry is the Kantowski-Sachs spacetime. However, the converse is not
necessarily true. For a fixed point to represent a homogeneous geometry and a
homogeneous matter distribution, not all the hat derivatives should
necessarily vanish. The prime example for such a case is the FLRW geometry,
for which $\hat{\phi}\neq 0$ but the hat derivatives of all the other
covariant geometric and thermodynamic quantities vanish.

The above discussion regarding the fixed points also extends to points on
invariant submanifolds, as they are also characterized by the vanishing of
certain covariant quantities.

##  4 Invariant submanifolds

In this section, we present the invariant submanifolds. Each of them divides
the entire phase space into two disjoint sectors. The interpretations of these
invariant submanifolds are given assuming that the expansion $\theta$ remains
finite. For a singular situation when $\theta\rightarrow\pm\infty$, the below
interpretations need not necessarily hold true.

The invariant submanifolds containing inhomogeneous solutions are as follows:

  1. 1.

$\mathcal{S}_{\phi}$: The invariant submanifold $x_{1}=0$ corresponds to
$\phi=0$. Recalling that $\phi$ measures the expansion of the 2D sheet along
the spacelike direction $e^{\mu}$, this invariant submanifold contains all the
solutions where a comoving observer measures zero expansion of the 2D sheet
along that direction. For the metric (2.1), it can be calculated that

$\phi=\frac{1}{B(t,r)}\left(2\frac{C_{,r}(t,r)}{C(t,r)}-\frac{A_{,r}(t,r)}{A(t,r)}\right)$
(4.1)

This might include solutions at spatial infinity $B(t,r)\rightarrow\infty$
(since at spatial infinity one can consider the radial lines to be parallel to
each other) and also Kantowski-Sachs solutions
$A=1,\,B=B(t),\,C=C(t),\,D(\theta)=\sin\theta$.

The stability of this invariant submanifold depends on the sign of the
quantity

$2x_{3}+x_{2}+4x_{3}^{2}=\frac{3\Sigma}{\theta}\left(1+\frac{3\Sigma}{\theta}\right)+\frac{3\mu}{\theta^{2}}.$
(4.2)

  2. 2.

$\mathcal{S}_{K}$: This is the submanifold of vanishing Gaussian curvature of
the 2D sheet. Let us define

$x_{K}\equiv\frac{9K}{\theta^{2}}=x_{2}-9x_{4}+\frac{9}{4}x_{1}^{2}-(1-x_{3})^{2}\,,$
(4.3)

where for the second equality we have utilized the expression (2.7). It is
easy to find out the evolution equation for $x_{K}$ using the evolution
equation for Gaussian curvature in (2.8):

$\frac{dx_{K}}{\epsilon
d\tau}=\frac{1}{9}x_{K}(2+3x_{2}+4x_{3}+12x_{3}^{2})\,.$ (4.4)

For the metric (2.1), it can be calculated that

$K=\frac{1}{C^{2}(t,r)}\frac{D^{\prime\prime}(y)}{D(y)}=\begin{cases}-\frac{1}{C^{2}(t,r)},\qquad
k=+1\\\ 0,\qquad k=0\\\ \frac{1}{C^{2}(t,r)},\qquad k=-1\end{cases}\,.$ (4.5)

This invariant submanifold, which contains cylindrically symmetric solutions
$D(y)=y$, divides the entire phase space into two disjoint subsections,
containing solutions of positive Gaussian curvature (i.e. spherically
symmetric solutions) and negative Gaussian curvature (i.e. hyperbolically
symmetric solutions). Apart from cylindrically symmetric solutions, it is also
possible for $\mathcal{S}_{K}$ to contain spherically or hyperbolically
symmetric solutions for a comoving observer at an infinite spatial distance
away from the origin, which is taken to be the symmetry centre555For
homogeneous spacetimes, there is no unique symmetry centres. The origin can be
placed anywhere.. This is because, even for spherically and hyperbolically
symmetric solutions, the Gaussian curvature of the shells will decrease as
they are further away from the origin. Existence of $\mathcal{S}_{K}$ means
that the Gaussian curvature of the 2D sheets cannot change sign.

The stability of this invariant submanifold depends on the quantity

$2+3x_{2}+4x_{3}+12x_{3}^{2}=2+6\frac{\Sigma}{\theta}+\frac{9\mu}{\theta^{2}}+\frac{27\Sigma^{2}}{\theta^{2}}\,.$
(4.6)

  3. 3.

$\mathcal{S}_{\mu}$: The invariant submanifold $x_{2}=0$ corresponds to the
vacuum scenario666In presence of matter there was a unique choice for the
timelike unit 4-vector $u^{\mu}$ that was the basis of the 1+3-decompositiion,
namely, the unit vector tangent to the worldline of the matter particles. In
vacuum, we do not have that choice. For cosmology, we can choose the “frame”,
i.e. the unit 4-vector $u^{\mu}$ such that the CMB dipole moment
vanishes.[36]. More discussion on this important issue is presented at the
beginning of the next section.. On this submanifold the Friedmann equation
(3.2) becomes

$x_{3}^{2}=1+\frac{3}{2}\frac{{}^{3}R}{\theta^{2}}.$ (4.7)

As expected, the vacuum case represents an invariant submanifold because there
is no matter production mechanism in classical physics. The stability of this
invariant submanifold depends on the sign of the quantity

$1-4x_{3}^{2}=1-\frac{9\Sigma^{2}}{\theta^{2}}=-3\left(1+2\frac{{}^{3}R}{\theta^{2}}\right).$
(4.8)

Following are some of the generic conclusions regarding the stability of this
invariant submanifold that we can reach from the above expression:

     * •

When shear is negligible ($|\Sigma/\theta|\ll 1$), $\mathcal{S}_{\mu}$ is of
attracting nature for and expanding LRS-II ($\epsilon=+1$) and of repelling
nature for a contracting LRS-II ($\epsilon=-1$). This is consistent with the
intuitive conclusion that the any pre-existing small amount of matter density
must dilute away with an expansion whereas it must go up with a contraction.
Note that this is in general not true in presence of non-negligible shear.

     * •

For spatially flat or positively curved LRS-II (${}^{3}R\geq 0$),
$\mathcal{S}_{\mu}$ is always of repelling nature for an expanding LRS-II
($\epsilon=+1$) and of attracting nature for a contracting LRS-II
($\epsilon=-1$). One may be tempted to interpret this result directly as any
pre-existing small amount of matter density is actually going up with an
expansion and diluting with contraction. This sounds counter-intuitive and one
can clearly see from the $\dot{\mu}$-equation in the system (2.14), matter
density definitely does go down with an expansion and go up with a
contraction. Indeed, it would be a wrong interpretation. Instead, the correct
interpretation should be that any pre-existing small amount of matter density
becomes dominant over the shear with expansion and becomes subdominant with
respect to the shear with contraction. This behaviour is similar to what one
expects for the case of homogeneous and anisotropic Bianchi spacetime in
presence of an isotropic fluid, where the matter density goes as
$\sim\frac{1}{a^{3}}$ and metric anisotropy goes as $\sim\frac{1}{a^{6}}$.
This same behaviour is preserved even when inhomogeneity is introduced in both
the metric and the matter sectors.

Note that one cannot conclude the same for spatially negatively curved LRS-II.

     * •

The above two conclusions are not in conflict as the isotropic vacuum limit
can only be supported with negative spatial curvature. The isotropic vacuum
limit with negative spatial curvature can be interpreted as an inhomogeneous
version of the well-known Milne solutions.

As explained at the end of the last section, a particular comoving observer
being on the invariant submanifold $\mathcal{S}_{\mu}$ does not mean that a
neighbouring comoving observer will also be on the same invariant submanifold.
In other words, even if a particular comoving observer experiences vacuum
locally, a neighbouring observer can experience locally a non-vanishing energy
density. This is because of a generally non-vanishing $\hat{\mu}$, which is
interpreted as some radial perturbations in $\mu$. The evolution equation of
this radial perturbation can be obtained from the last two equations of (2.14)

$\ddot{\hat{\mu}}+\left(\frac{4}{3}\theta+\Sigma\right)\dot{\hat{\mu}}-\left(\frac{4}{9}\theta^{2}+\frac{5}{2}\Sigma^{2}+\frac{2}{3}\theta\Sigma+\mathcal{E}\right)\hat{\mu}=0\,.$
(4.9)

However, it is more illuminating to write the evolution equation for the
dimensionless quantity $y_{2}$ here,

$\frac{dy_{2}}{\epsilon d\tau}=-\frac{1}{3}y_{2}(1+2x_{3}-6x_{3}^{2})\,.$
(4.10)

Whether the radial perturbations decay ($y_{2}\rightarrow 0$) or not during an
expanding or a contracting LRS-II, depends on the sign of the quantity
$(1+2x_{3}-6x_{3}^{2})$. If the radial perturbation decays, the situation
asymptotically tends to be vacuum for _all_ the observers.

  4. 4.

$\mathcal{S}_{\Sigma\mathcal{E}}$: The invariant submanifold
$\\{x_{3},x_{4}\\}=\\{0,0\\}$ corresponds to
$\\{\Sigma,\mathcal{E}\\}=\\{0,0\\}$. On this submanifold the Friedmann
equation (3.2) becomes

$x_{2}=1+\frac{3}{2}\frac{{}^{3}R}{\theta^{2}}.$ (4.11)

The stability of this invariant submanifold can be investigated from the
Jacobian

$\left(\begin{smallmatrix}\frac{\partial}{\partial
x_{3}}\left(\frac{dx_{3}}{d\tau}\right)&\frac{\partial}{\partial
x_{4}}\left(\frac{dx_{3}}{d\tau}\right)\\\ \frac{\partial}{\partial
x_{3}}\left(\frac{dx_{4}}{d\tau}\right)&\frac{\partial}{\partial
x_{4}}\left(\frac{dx_{4}}{d\tau}\right)\end{smallmatrix}\right)_{x_{3}=0,x_{4}=0},$

whose eigenvalues are
$\frac{\epsilon}{12}\left(3x_{2}-4\pm\sqrt{x_{2}^{2}+24x_{2}}\right)$.
Following are some of the generic conclusions that can be drawn:

     * •

In a vacuum ($x_{2}=0$), $\mathcal{S}_{\Sigma\mathcal{E}}$ is always of
attracting (repelling) nature for an expanding (contracting) LRS-II.

     * •

In the presence of matter, the stability of $\mathcal{S}_{\Sigma\mathcal{E}}$
depends entirely on the matter density $\mu$ (or, equivalently, the
3-curvature ${}^{3}R$).

       * –

For an expanding LRS-II spacetime, $\mathcal{S}_{\Sigma\mathcal{E}}$ is
attracting when $0\leq x_{2}<3-\sqrt{7}$.

       * –

For a contracting LRS-II spacetime, $\mathcal{S}_{\Sigma\mathcal{E}}$ is
attracting when $x_{2}>3+\sqrt{7}$.

       * –

For a spatially flat (${}^{3}R=0$) expanding or contracting LRS-II, there will
always be one direction in the phase space along which this invariant
submanifold is attracting, and one direction along which it is repelling.

For the special case of a spherically symmetric spacetime (i.e. the Gaussian
curvature $K>0$), $\\{\Sigma,\mathcal{E}\\}=\\{0,0\\}$ implies a locally FLRW
spacetime. A comoving observer passing through an evolutionary phase with
$\\{\Sigma,\mathcal{E}\\}=\\{0,0\\}$ (i.e. on the invariant submanifold
$\\{x_{3},x_{4}\\}=\\{0,0\\}$) experiences locally an FLRW spacetime. However,
the observer will experience, in general, non-vanishing radial perturbations
in the spacetime geometry and matter density, as given by the quantities
$\\{\hat{\theta},\hat{\Sigma},\hat{\mathcal{E}},\hat{\mu}\\}$. The evolution
equations for $\\{\hat{\theta},\hat{\mu}\\}$, can be obtained from the last
two equations of (2.14)

$\displaystyle\ddot{\hat{\theta}}+\theta\dot{\hat{\theta}}-\left(\frac{1}{3}\theta^{2}+\mu\right)\hat{\theta}=\frac{2}{3}\theta\hat{\mu}\,,$
(4.12a)
$\displaystyle\ddot{\hat{\mu}}+\frac{4}{3}\theta\dot{\hat{\mu}}-\left(\frac{4}{9}\theta^{2}+\frac{7}{6}\mu\right)\hat{\mu}=2\mu\theta\hat{\theta}\,.$
(4.12b)

Subsequently, $\\{\hat{\Sigma},\hat{\mathcal{E}}\\}$ can be found from
Eq.(2.4). These non-vanishing radial perturbations will stop the geometry from
being totally homogeneous. If a comoving observer at $r=r_{*}$ experiences
locally an FLRW geometry, a neighbouring comoving observer $r_{1}=r_{*}+\delta
r$ will not, in general, experience an FLRW geometry. Similar arguments apply
for the cases of cylindrically and hyperbolically symmetric spacetimes.

The points on the invariant submanifolds
$\mathcal{S}_{\phi},\,\mathcal{S}_{K},\,\mathcal{S}_{\mu},\,\mathcal{S}_{\Sigma\mathcal{E}}$
correspond, in general, to inhomogeneous solutions. Take, for example, the
invariant submanifold $\mathcal{S}_{\mu}$. If a shell is initially on
$\mathcal{S}_{\mu}$, then it will always stay on it. Physically, if the matter
density $\mu$ vanishes on a particular shell, then it will always remain
vanishing on that shell forever. At the same time, the other shells will, in
general, have a nonzero matter density, and their matter densities will also
change with time. Similar interpretation goes for the invariant submanifolds
$\mathcal{S}_{\phi}$ and $\mathcal{S}_{\Sigma\mathcal{E}}$.

Below we present a special invariant submanifold, on which all the points
correspond to homogeneous solutions.

  1. 1.

$\mathcal{S}_{\rm hom}$: A more restricted invariant submanifold is given by
$\\{x_{3},x_{4},y_{2},z\\}=\\{0,0,0,0\\}$, which corresponds to homogeneous
solutions. The points on this invariant submanifold that lie within the region
$x_{K}>0$, which correspond to spherically symmetric solutions, represent the
homogeneous FLRW geometry. The points on this submanifold that lie within the
region $x_{K}\leq 0$ represent homogeneous geometries that are hyperbolically
or cylindrically symmetric counterparts of FLRW geometry.

As one can see from the dynamical system (3.6), inhomogeneity creeps into the
dynamics through the dynamical variables $y_{2}$ and $z$. Note that, by
itself, $\\{y_{2},z\\}=\\{0,0\\}$ is _not_ an invariant submanifold. Even if
we start initially with a homogeneous geometry ($z=0$) and a homogeneous
matter distribution ($y_{2}=0$), non-vanishing sheet expansion and shear
($x_{1}x_{3}\neq 0$) can induce inhomogeneity in the geometry ($z$), which, in
turn, induces inhomogeneity in the matter distribution ($y_{2}$) by virtue of
the $x_{2}z$ term. However, if we start with the invariant sub-manifold
$\\{x_{3},x_{4},y_{2},z\\}=\\{0,0,0,0\\}$, then there is no way for
inhomogeneity to creep into the system. This implies that a homogeneous
universe stays homogeneous in the absence of any radial perturbations, which
is consistent with $\mathcal{S}_{\rm hom}$ being an invariant submanifold.

Note the difference between $\mathcal{S}_{\Sigma\mathcal{E}}$ and
$\mathcal{S}_{\rm hom}$. $\mathcal{S}_{\Sigma\mathcal{E}}$ had vanishing shear
$\Sigma$ and electric Weyl tensor $\mathcal{E}$
($\\{x_{3},x_{4}\\}=\\{0,0\\}$), but that did not automatically imply
vanishing $\hat{\Sigma}$ and $\hat{\mathcal{E}}$
($\\{y_{3},y_{4}\\}\neq\\{0,0\\}$). Not all the comoving shells were
experiencing an evolutionary phase with vanishing shear and electric Weyl
tensor simultaneously. Therefore, the solutions in
$\mathcal{S}_{\Sigma\mathcal{E}}$ were inhomogeneous solutions. On the other
hand, $\mathcal{S}_{\rm hom}$ also has vanishing shear $\Sigma$ and electric
Weyl tensor $\mathcal{E}$ ($\\{x_{3},x_{4}\\}=\\{0,0\\}$), and by virtue of
also having $\\{y_{2},z\\}=\\{0,0\\}$, the constraint equations (3.5)
guarantee vanishing $\\{\hat{\Sigma}$ and $\hat{\mathcal{E}}\\}$
($\\{y_{3},y_{4}\\}\neq\\{0,0\\}$). In this case, _all_ the comoving observers
experience an evolutionary phase with vanishing shear and electric Weyl
tensor, which is precisely the case of a homogeneous FLRW spacetime geometry.
The invariant submanifold $\mathcal{S}_{\rm hom}$ is actually a subset of
$\mathcal{S}_{\Sigma\mathcal{E}}$: $\mathcal{S}_{\rm
hom}\subset\mathcal{S}_{\Sigma\mathcal{E}}$, that is singled out by the
vanishing of the radial perturbations.

For a spherically symmetric solution on the invariant submanifold
$\mathcal{S}_{\Sigma\mathcal{E}}$, the local spacetime geometry with respect
to the comoving observer at $r=r_{*}$ is FLRW. But because of the non-
vanishing radial perturbations, a neighbouring comoving observer at
$r_{1}=r_{*}+\delta r$ will not, in general, experience an FLRW geometry. For
the invariant submanifold $\mathcal{S}_{\rm hom}$, there are no such radial
perturbations. If the comoving observer at $r_{0}=r_{*}$ experiences an FLRW
geometry, a neighbouring comoving observer at $r_{1}=r_{*}+\delta r$ will also
experience the same FLRW geometry. This argument is inductive. One can then
focus on this observer at $r_{1}=r_{*}+\delta r$ and apply the same argument.
By the principle of induction, this argument can be extended to _all_ the
comoving observers, which implies that all the comoving observers experience
the same local FLRW geometry. Similar argument applies to cylindrically and
hyperbolically symmetric spacetimes.

The stability of this invariant submanifold can be investigated from the
Jacobian

$\left(\begin{smallmatrix}\frac{\partial}{\partial
x_{3}}\left(\frac{dx_{3}}{d\tau}\right)&\frac{\partial}{\partial
x_{4}}\left(\frac{dx_{3}}{d\tau}\right)&\frac{\partial}{\partial
y_{2}}\left(\frac{dx_{3}}{d\tau}\right)&\frac{\partial}{\partial
z}\left(\frac{dx_{3}}{d\tau}\right)\\\ \frac{\partial}{\partial
x_{3}}\left(\frac{dx_{4}}{d\tau}\right)&\frac{\partial}{\partial
x_{4}}\left(\frac{dx_{4}}{d\tau}\right)&\frac{\partial}{\partial
y_{2}}\left(\frac{dx_{4}}{d\tau}\right)&\frac{\partial}{\partial
z}\left(\frac{dx_{4}}{d\tau}\right)\\\ \frac{\partial}{\partial
x_{3}}\left(\frac{dy_{2}}{d\tau}\right)&\frac{\partial}{\partial
x_{4}}\left(\frac{dy_{2}}{d\tau}\right)&\frac{\partial}{\partial
y_{2}}\left(\frac{dy_{2}}{d\tau}\right)&\frac{\partial}{\partial
z}\left(\frac{dy_{2}}{d\tau}\right)\\\ \frac{\partial}{\partial
x_{3}}\left(\frac{dz}{d\tau}\right)&\frac{\partial}{\partial
x_{4}}\left(\frac{dz}{d\tau}\right)&\frac{\partial}{\partial
y_{2}}\left(\frac{dz}{d\tau}\right)&\frac{\partial}{\partial
z}\left(\frac{dz}{d\tau}\right)\\\
\end{smallmatrix}\right)_{x_{3}=0,x_{4}=0,y_{2}=0,z=0},$

whose eigenvalues are

$\frac{\epsilon}{12}\left(3x_{2}-4\pm\sqrt{x_{2}^{2}+24x_{2}}\right),\,\,\frac{\epsilon}{12}\left(5x_{2}-4\pm\sqrt{x_{2}^{2}+24x_{2}}\right)$

The first pair of eigenvalues are exactly the ones we obtained for
$\mathcal{S}_{\Sigma\mathcal{E}}$, whereas it can be checked that the
eigenvectors corresponding to the last pair of eigenvalues lie entirely on
$\mathcal{S}_{\Sigma\mathcal{E}}$. Since $\mathcal{S}_{\rm
hom}\subset\mathcal{S}_{\Sigma\mathcal{E}}$, we conclude that a part of the
stability of $\mathcal{S}_{\rm hom}$ comes from the flow in the vicinity of
$\mathcal{S}_{\Sigma\mathcal{E}}$, which is given by the first pair of
eigenvalues. The other part of the stability of $\mathcal{S}_{\rm hom}$ comes
actually from the flow within $\mathcal{S}_{\Sigma\mathcal{E}}$, which is
given by the last pair of eigenvalues. This last pair of eigenvalues are
essentially related to the dynamics of the radial perturbations on
$\mathcal{S}_{\Sigma\mathcal{E}}$. When both of them are negative, one can say
that the radial perturbations tend to die out, driving the geometry towards a
globally homogeneous one.

The stability of this invariant submanifold dictates the evolution of some
LRS-II geometries that are almost homogeneous. Following are some of the
generic conclusions that can be drawn:

     * •

In a vacuum ($x_{2}=0$), $\mathcal{S}_{\rm hom}$ is always of attracting
(repelling) nature for an expanding (contracting) LRS-II.

     * •

In presence of matter, the stability of $\mathcal{S}_{\rm hom}$ depends
entirely on the matter density $\mu$ (or, equivalently, the 3-curvature
${}^{3}R$).

       * –

For an expanding LRS-II spacetime, $\mathcal{S}_{\rm hom}$ is attracting when
$0\leq x_{2}<\frac{4}{3}-\frac{\sqrt{10}}{3}$.

       * –

For a contracting LRS-II spacetime, $\mathcal{S}_{\rm hom}$ is attracting when
$x_{2}>3+\sqrt{7}$.

       * –

For a spatially flat (${}^{3}R=0$) expanding or contracting LRS-II, there will
always be two directions in the phase space along which this invariant
submanifold is attracting, and two directions along which it is repelling.

Before ending this section, some important discussions are in order:

  * •

Firstly, it is worth mentioning that for LRS-II, spatially flat solutions do
_not_ in general represent an invariant submanifold. If we define
$x_{5}\equiv\frac{3}{2}\frac{{}^{3}R}{\theta^{2}}$, then, using the Friedmann
equation (3.2), one can derive

$\frac{dx_{5}}{\epsilon
d\tau}=\frac{1}{3}x_{5}(x_{2}+2x_{3}+4x_{3}^{2})+\frac{3}{2}x_{3}(x_{1}^{2}+2y_{1}),$
(4.13)

where the second term, which depends explicitly on the shear, prevents
$x_{5}=0$ from being an invariant submanifold. The 3-curvature locally
experienced by a comoving observer can change between being positive and
negative. Only when the shear vanishes, spatially flat solutions represent an
invariant submanifold.

  * •

For the points at the intersection between the invariant submanifolds
$\mathcal{S}_{\Sigma\mathcal{E}}$ and $S_{\mu}$, the evolution asymptotically
homogenizes with expansion whereas inhomogeneities grow with contraction. It
is easy to find out the behaviour of the radial perturbations for points
belonging to the intersection of these two invariant submanifolds, as putting
$\\{x_{2},x_{3},x_{4}\\}=\\{0,0,0\\}$ leads to a linear system for $y_{2}$ and
$z$:

$\displaystyle\frac{dy_{2}}{\epsilon d\tau}=-\frac{1}{3}y_{2}\,,$ (4.14a)
$\displaystyle\frac{dz}{\epsilon d\tau}=-\frac{1}{6}y_{2}-\frac{1}{3}z\,,$
(4.14b)

whose solution is
$\\{y_{2}(\tau)\,,z(\tau)\\}=\\{c_{1}e^{-\epsilon\tau/3}\,,c_{1}e^{-\epsilon\tau/3}\left(c_{2}-\frac{1}{6}\epsilon\tau\right)\\}$.
Using the propagation equations (3.5), we can deduce the complete behaviour of
the radial perturbations

$y_{2}(\tau)=c_{1}e^{-\epsilon\tau/3}=9y_{4}(\tau)\,,\quad
z(\tau)=c_{1}e^{-\epsilon\tau/3}\left(c_{2}-\frac{1}{6}\epsilon\tau\right)=y_{3}(\tau)$
(4.15)

Evidently, the phase flow is towards $\mathcal{S}_{\rm hom}$ for $\epsilon=+1$
and away from $\mathcal{S}_{\rm hom}$ for $\epsilon=-1$.

##  5 Finite fixed points

We can obtain five isolated fixed points and two lines of fixed points for the
system (3.6), which are listed in Table 1.

Fixed point | $(x_{1},x_{2},x_{3},x_{4},y_{2},z)$ | Stability |  | Actual
---
spacetime geometry
$\mathcal{L}_{1}$ | $(\lambda,0,0,0,0,0)$ |  | Attractor for expanding LRS-II
---
Repeller for contracting LRS-II
| Minkowski
---
(Milne w.r.t. comoving observer)
$\mathcal{L}_{2}$ | $\left(\lambda,0,-\frac{1}{2},0,0,\frac{1}{2}\lambda\right)$ | Saddle | Minkowski
$\mathcal{P}_{1}$ | $\left(0,0,-1,-\frac{4}{9},0,0\right)$ |  | Repeller for expanding LRS-II
---
Attractor for contracting LRS-II
| Homogeneous
---
Spatially flat
Vacuum
$\mathcal{P}_{2}$ | $(0,1,0,0,0,0)$ | Saddle |  | Homogeneous
---
Spatially flat
Matter dominated
$\mathcal{P}_{3}$ | $\left(0,0,\frac{1}{4},-\frac{1}{16},0,0\right)$ | Saddle |  | Homogeneous
---
Spatially negatively curved
Vacuum
$\mathcal{P}_{4}$ | $(0,0,1,0,0,0)$ | Saddle | Minkowski

Table 1: Finite fixed points that can be obtained from the system (3.6). The
stability can be analyzed by considering the Jacobian eigenvalues.
$\mathcal{L}_{1}$ and $\mathcal{L}_{2}$ are lines of fixed points, with
$\lambda$ denoting the parameter along the line.

Let us discuss below what kind of solutions these fixed points represent. In
interpreting the spacetime geometry corresponding to these fixed points, we
must keep in mind that, when a matter is present, we have a preferred choice
of the timelike 4-vector $u^{\mu}$ at each spacetime point. The appropriate
foliation is the family of hypersurfaces of constant energy density. We can,
then, choose the 4-vector $u^{\mu}$ to be normal to this hypersurface of
constant energy density. This choice is unique and has an invariant physical
meaning. The interpretations of the spacetime geometry for fixed points that
correspond to nonvacuum solutions are made with respect to the a comoving
observer having this 4-velocity.

When there is no matter, there is no unique choice of foliation, as any
spacelike hypersurface is a hypersurface of uniform energy density equal to
zero. Consequently, there is no unique choice for the timelike 4-vector
$u^{\mu}$. The interpretations of the spacetime geometry for fixed points that
correspond to a vacuum solution are, therefore ambiguous; one has the freedom
to switch to a different $u^{\mu}$, which amounts to a different choice of
congruence. Therefore, for those vacuum fixed points whose interpretations are
not straightforward with respect to a comoving observer (which are, in fact,
the case for all the vacuum fixed points except for $\mathcal{L}_{1}$), we
will adopt alternative approaches to find out the _actual_ spacetime geometry
corresponding to the fixed points. These actual spacetime geometries need not
necessarily correspond to the spacetime geometry seen by the comoving
observers of the LRS-II geometry.

  1. 1.

$\mathcal{L}_{1}$: The line of fixed points $\mathcal{L}_{1}$ lies at the
intersections of the invariant submanifold $\mathcal{S}_{\mu}$ and
$\mathcal{S}_{\rm hom}$, which implies this is a vacuum and homogeneous
solution with respect to the comoving observers. The constraint equation (3.2)
shows that such a fixed point is possible only when the 3-curvature is
negative ${}^{3}R<0$, which identifies this solution as Milne with respect to
the comoving observer.

On the other hand, since it is a vacuum case, there is no preferred choice for
congruence of timelike geodesics. It is known that one can redefine the
coordinates in such a way, i.e. switch to a different congruence of geodesics
such that the Milne geometry is recast in a Minkowski form. Another way of
arguing that the geometry is effectively Minkowski is that since both matter
density (and hence, by virtue of the Einstein field equation, the Ricci
tensor) vanishes and electric Weyl tensor vanishes
($x_{2}=x_{4}=0\Rightarrow\mu=\mathcal{E}=0$), the Riemann curvature tensor
must vanish. The spacetime geometry is actually Minkowski, although with
respect to the LRS-II comoving observers ($u^{\mu}=(1,0,0,0)$), locally it is
Milne.

Therefore, homogeneous Milne arises as an attractor in an expanding LRS-II and
as a repelled in a contracting LRS-II. The attractive nature of the Milne
solution in GR was found in Refs.[46, 47] in the study of LTB dynamics. We
show that the same conclusion holds for the generic case of LRS-II.

  2. 2.

$\mathcal{L}_{2}$: The line of fixed points $\mathcal{L}_{2}$ falls on the
invariant submanifold $\mathcal{S}_{\mu}$. If a comoving observer is going
through an evolutionary phase characterized by this fixed point, the matter
density locally with respect to the observer vanishes. Moreover, since
$y_{2}=0$, i.e. $\hat{\mu}=0$, the comoving observer also sees no radial
perturbation in matter density. Also, the line is characterized by the
relation

$z=\frac{1}{2}x_{1}\,\,\Rightarrow\hat{\theta}=\frac{1}{2}\phi\theta.$ (5.1)

Also, $x_{3}=-\frac{1}{2}$ implies the following relation between the local
expansion and shear

$\theta+3\Sigma=0.$ (5.2)

$x_{3}=-\frac{1}{2}$ also implies, via the constraint equation (3.2), that
locally the 3-curvature is negative ${}^{3}R<0$. Since $x_{4}=0$, the electric
Weyl tensor locally vanishes. Moreover, from the propagation equation, (3.5)
one can see that $y_{2}=0$ and $x_{4}=0$ also makes $y_{4}=0$, i.e. the
comoving observer does not experience any radial perturbation of the electric
Weyl tensor. Unlike $\mathcal{L}_{1}$, it is, however, not straightforward to
interpret the spacetime geometry with respect to a comoving LRS-II observer.

On the other hand, since it is a vacuum case, the Ricci tensor vanishes by
virtue of the Einstein field equation. Since the Weyl curvature tensor also
vanishes, this definitely means the vanishing of the Riemann curvature tensor.
Therefore, the spacetime geometry corresponding to this solution is Minkowski.
A Minkowski spacetime satisfies the relations (5.1),(5.2), but is apparently
inconsistent with the fact that we are getting ${}^{3}R<0$. However, since it
is a vacuum solution, there is no preferred choice for congruence of timelike
geodesics. We can switch to a different congruence, i.e. redefine the
coordinates in such a way that the Minkowski geometry is recast in a Milne
form, reconciling the fact that we are getting a negative 3-curvature.

  3. 3.

$\mathcal{P}_{1}$: The isolated fixed point $\mathcal{P}_{1}$ lies at the
intersection of the invariant submanifolds $\mathcal{S}_{\phi}$,
$\mathcal{S}_{\mu}$ and $\mathcal{S}_{K}$. It represents a vacuum solution. It
can be found from the constraint equation (3.2) that ${}^{3}R=0$, i.e. the
solution is spatially flat. Since $y_{2}=0$ ($\hat{\mu}=0$) and $z=0$
($\hat{\theta}=0$), the comoving observer does not experience any radial
perturbation in matter density and the rate of expansion is the same for all
the comoving observers. From the propagation equations (3.5) one can also
calculate that $\\{y_{1},y_{3},y_{4}\\}=\\{0,0,0\\}$, i.e.
$\\{\hat{\phi},\hat{\Sigma},\hat{\mathcal{E}}\\}=\\{0,0,0\\}$. This implies
that the comoving observer experiences no radial perturbation in _any_ of the
dynamical quantities. All the LRS-II comoving observers see the same spacetime
geometry, i.e. the situation is homogeneous. With respect to the LRS-II
comoving observers, the spacetime geometry is spatially flat, homogeneous and
anisotropic.

There is another way of arriving at the conclusion that the spacetime geometry
must be homogeneous. This fixed point is a vacuum solution of GR. An extension
of Birkhoff’s theorem holds for LRS-II spacetimes in GR [44, 48], stating that
the solution possesses an additional symmetry which makes the spacetime either
spatially homogeneous with $\phi=0$ or static with $\theta=\Sigma=0$. The
condition $\theta=\Sigma=0$ is clearly not satisfied for this fixed point
whereas the condition $\phi=0$ is satisfied, suggesting that the spacetime
must be homogeneous.

  4. 4.

$\mathcal{P}_{2}$: The isolated fixed point $\mathcal{P}_{2}$ lies at the
intersection of the invariant submanifolds $\mathcal{S}_{\phi}$,
$\mathcal{S}_{\rm hom}$ and $\mathcal{S}_{K}$. From the Friedmann constraint
(3.2), one can see that the 3-curvature corresponding to this fixed point
solution is zero ${}^{3}R=0$. Therefore, it represents a spatially flat
homogeneous solution. A comoving observer describes the spacetime geometry to
be a spatially flat and homogeneous. Note that, since this is not a vacuum
solution, there is a preferred choice of the congruence of timelike geodesics,
namely, those that are along the worldline of the comoving observers.

The time-evolution corresponding to this solution can be found out from
eq.(3.3), which when written in terms of the Hubble parameter
$H\equiv\frac{\theta}{3}\equiv\frac{\dot{a}}{a}$, gives

$\frac{\dot{H}}{H^{2}}=-\frac{3}{2}\,.$ (5.3)

This corresponds to the time evolution

$a(t)\sim\begin{cases}(+t)^{2/3}\,,\qquad(t>0)\qquad\text{for an expanding
universe}\\\ (-t)^{2/3}\,,\qquad(t<0)\qquad\text{for a contracting
universe}\end{cases}$ (5.4)

Since ${}^{3}R=\Sigma=0$, the original Friedmann equation (2.5) gives, in
terms of the Hubble parameter,

$3H^{2}=\mu\,.$ (5.5)

We, therefore, get the matter-dominated epoch as a saddle point, i.e.
intermediate state of LRS-II evolution with dust.

  5. 5.

$\mathcal{P}_{3}$: The isolated fixed point $\mathcal{P}_{3}$ lies at the
intersection of the invariant submanifolds $\mathcal{S}_{\phi}$,
$\mathcal{S}_{\mu}$ and $\mathcal{S}_{K}$. From the Friedmann constraint
(3.2), one can see that the 3-curvature corresponding to this fixed point
solution is negative ${}^{3}R<0$. It is a vacuum solution. Since $y_{2}=z=0$,
a comoving observer does not experience any radial perturbation in matter
distribution and expansion rate ($\hat{\mu}=\hat{\theta}=0$). From the
propagation equations (3.5) one can also calculate that
$\\{y_{3},y_{4}\\}=\\{0,0\\}$, i.e.
$\\{\hat{\Sigma},\hat{\mathcal{E}}\\}=\\{0,0\\}$ but $y_{1}\neq 0$,
i.e.$\hat{\phi}\neq 0$.

However, just like the point $\mathcal{P}_{1}$, the fixed point
$\mathcal{P}_{3}$ is also a vacuum LRS-II solution of GR and therefore
possesses an additional symmetry which makes the spacetime either spatially
homogeneous or static. Since $\phi=0$ for this fixed point, this suggests that
the spacetime must be homogeneous [44, 48].

  6. 6.

$\mathcal{P}_{4}$: The isolated fixed point $\mathcal{P}_{4}$ also lies at the
intersection of the invariant submanifolds $\mathcal{S}_{\phi}$,
$\mathcal{S}_{\mu}$ and $\mathcal{S}_{K}$. From (3.2), the spatial 3-curvature
vanishes ${}^{3}R=0$. Therefore, this is a spatially flat vacuum solution.
Since this fixed point corresponds to a vacuum situation and the electric part
of the Weyl tensor also vanishes ($x_{4}=0$), the Riemann curvature tensor
vanishes, which specifies the actual spacetime geometry to be Minkowski.

Before moving on to the next section about the fixed points at the infinity of
the phase space, some discussions are in order regarding the finite fixed
points.

  * •

In table 1, we have listed the _actual_ spacetime geometry corresponding to
the fixed points representing vacuum solutions. we emphasize again that these
actual spacetime geometries are not always the geometry seen by a comoving
observer. With respect to a different congruence, these actual geometries
should be apparent. The idea is similar to when one tries to rewrite the
static Schwarschild or Minkowski metric with respect to the coordinates
attached to a comoving observer, one gets back a non-static metric, but the
actual geometry still remains Schwarzschild or Minkowski. In the absence of
any matter, the observer always has this freedom. No such ambiguity appears
for fixed points that represent nonvacuum solutions.

  * •

One notices that all the isolated finite fixed points in the phase space of
LRS-II in the presence of dust fall at the intersection of the invariant
submanifolds $\mathcal{S}_{\phi}$ and $\mathcal{S}_{K}$. The corresponding
evolutionary phases are achieved by a comoving observer only when (s)he
experiences zero expansion of the 2-sheet along the preferred radial direction
and zero Gaussian curvature of the shell. We have given the expressions of
$\phi$ and $K$ in terms of the metric functions for the LRS-II metric (2.1) in
Eqs.(4.1) and (4.5) respectively. From those expressions, it is not easy to
interpret the conditions $\phi=0$ and $K=0$ from the coordinate point of view
in its full generality. The interpretation becomes easier when we specialize
to a particular metric, e.g. LTB, as we will see in section 7.

  * •

$\theta\rightarrow-\infty$ represents a _shell-focusing_ singularity, when the
congruence of comoving geodesics runs into a caustic. On the other hand,
$\theta\rightarrow+\infty$ can be interpreted as a _shell-diverging_
singularity, when the comoving geodesics are infinitely away from each It is
possible for all the finite fixed points listed in the table 1 to represent
singular situations. Take, for example, the fixed point $\mathcal{P}_{2}$. It
can represent a situation where $\theta$ diverges faster than
$\Sigma,\,\mathcal{E}$ and at the same rate as $\sqrt{3\mu}$. This is
precisely the big-bang or big-crunch singularity ingrained in an FLRW
evolution. In fact, if an expansion normalized dynamical variable $x_{i}$
vanishes for a fixed point, one possible interpretation of the fixed point can
be a singular situation where $\theta$ diverges faster than whatever covariant
quantity is there in the numerator of $x_{i}$.

  * •

Perhaps the most remarkable observation is, even if we have started with an
inhomogeneous spacetime in the presence of an inhomogeneous matter (dust)
distribution, all the finite isolated fixed points that we could find actually
represent a homogeneous solution. As we will see in the next section, this is
not the case with fixed points at the infinity of the phase space. It is
possible for some of those asymptotic fixed points (belonging to the family
$\mathcal{B}^{(4)}$; see section 6) to represent inhomogeneous geometries.

##  6 Fixed points at infinity

To investigate the global behaviour of the phase space, it is necessary to
find fixed points that may reside at the infinity of the phase space. To find
the fixed points at the infinity of the phase space, we compactify the 6D
phase space spanned by the dynamical variables
$x_{1},\,x_{2},\,x_{3},\,x_{4},\,y_{2},\,z$ using the usual Poincare
compactification prescription. The procedure is as follows. First, we define
spherical polar coordinates as follows:

$\displaystyle x_{1}$
$\displaystyle=r\cos{\alpha}\cos{\beta}\cos{\gamma}\cos{\delta}\cos{\theta}\,,$
(6.1) $\displaystyle x_{2}$
$\displaystyle=r\sin{\alpha}\cos{\beta}\cos{\gamma}\cos{\delta}\cos{\theta}\,,$
$\displaystyle x_{3}$ $\displaystyle=r\sin{\theta}\,,$ $\displaystyle x_{4}$
$\displaystyle=r\sin{\beta}\cos{\gamma}\cos{\delta}\cos{\theta}\,,$
$\displaystyle y_{2}$ $\displaystyle=r\sin{\gamma}\cos{\delta}\cos{\theta}\,,$
$\displaystyle z$ $\displaystyle=r\cos{\delta}\cos{\theta}\,.$

Now, we have a 6D phase space spanned by the radial variable $r$ and five
angular variables $\alpha,\,\beta,\,\gamma,\,\delta,\,\theta$. Next, we
compactify the radial direction as

$r\rightarrow R=\frac{r}{\sqrt{1+r^{2}}}\,.$ (6.2)

Now that the phase space is compactified, we go back to Cartesian coordinates
again by defining

$\displaystyle X_{1}$
$\displaystyle=R\cos{\alpha}\cos{\beta}\cos{\gamma}\cos{\delta}\cos{\theta}\,,$
(6.3) $\displaystyle X_{2}$
$\displaystyle=R\sin{\alpha}\cos{\beta}\cos{\gamma}\cos{\delta}\cos{\theta}\,,$
$\displaystyle X_{3}$ $\displaystyle=R\sin{\theta}\,,$ $\displaystyle X_{4}$
$\displaystyle=R\sin{\beta}\cos{\gamma}\cos{\delta}\cos{\theta}\,,$
$\displaystyle Y_{2}$ $\displaystyle=R\sin{\gamma}\cos{\delta}\cos{\theta}\,,$
$\displaystyle Z$ $\displaystyle=R\cos{\delta}\cos{\theta}\,.$

The dynamical variables $X_{1},\,X_{2},\,X_{3},\,X_{4},\,Y_{2},\,Z$ are now
compact, limited within the 6-sphere given by

$0\leq X_{1}^{2}+X_{2}^{2}+X_{3}^{2}+X_{4}^{2}+Y_{2}^{2}+Z^{2}\leq 1$ (6.4)

The direct transformation between the compact and noncompact dynamical
variables are

$\displaystyle X_{1}=\frac{x_{1}}{\sqrt{1+r^{2}}}\,,$ $\displaystyle
X_{2}=\frac{x_{2}}{\sqrt{1+r^{2}}}\,,\quad
X_{3}=\frac{x_{3}}{\sqrt{1+r^{2}}}\,,\quad
X_{4}=\frac{x_{4}}{\sqrt{1+r^{2}}}\,,$ (6.5) $\displaystyle
Y_{2}=\frac{y_{2}}{\sqrt{1+r^{2}}}\,,\quad Z=\frac{z}{\sqrt{1+r^{2}}}\,.$

and the inverse transformation is

$\displaystyle x_{1}=\frac{X_{1}}{\sqrt{1-R^{2}}}\,,$ $\displaystyle
x_{2}=\frac{X_{2}}{\sqrt{1-R^{2}}}\,,\quad
x_{3}=\frac{X_{3}}{\sqrt{1-R^{2}}}\,,\quad
x_{4}=\frac{X_{4}}{\sqrt{1-R^{2}}}\,,$ (6.6) $\displaystyle
y_{2}=\frac{Y_{2}}{\sqrt{1-R^{2}}}\,,\quad z=\frac{Z}{\sqrt{1-R^{2}}}\,.$

The above transformations can be used to write the dynamical system (3.6) in
terms of the compact dynamical variables (6.5). The expressions are huge and
we do not find it necessary to write their explicit forms here. One encounters
a singularity of the dynamical system at the boundary of the phase space
$R\rightarrow 1$. This can be regularized by redefining the phase space time
variable as

$\tau\rightarrow\bar{\tau}:d\bar{\tau}=\frac{d\tau}{1-R^{2}}\,.$ (6.7)

In terms of $\bar{\tau}$, one can get a regular dynamical system in terms of
the compact dynamical variables. As we already investigated the finite fixed
points in the previous section, in this section we specifically list the
asymptotic fixed points, i.e. fixed points with $R=1$.

There are three different families of asymptotic fixed points. They are as
follows:

  1. 1.

A pair of isolated fixed points given by
$\mathcal{O}_{\pm}\equiv(X_{1},X_{2},X_{3},X_{4},Y_{2},Z)=(0,0,\pm 1,0,0,0)$,
which corresponds to $x_{3}=\frac{3\Sigma}{2\theta}\rightarrow\pm\infty$. Both
these points are saddles. Both these asymptotic fixed points lie on the
invariant submanifolds $\mathcal{S}_{\phi}$ and $\mathcal{S}_{\mu}$. They are
vacuum solutions. From Eq.(4.3), we find that these fixed points lie on the
region $x_{K}<0$, i.e. the corresponding spacetime geometry, as seen by a
comoving observer, must be cylindrically symmetric. Since this is a vacuum
case, the Ricci tensor vanishes by virtue of the Einstein field equation.
$X_{4}=0$ imply that the electric Weyl tensor also vanishes. Together, they
imply the vanishing of the Riemann tensor. The actual spacetime geometry is
therefore Minkowski. A comoving observer, of course, does _not_ see the
geometry as Minkowski, as evident from the fact that
$x_{3}\equiv\frac{3\Sigma}{2\theta}=\pm\infty$. Nonetheless, with respect to a
different timelike congruence, the Minkowski geometry should be apparent.

  2. 2.

A 1-parameter family of fixed points $\mathcal{A}^{(1)}$ given by
$X_{2}=X_{4}=Y_{2}=0,\,Z=3X_{1},\,10X_{1}^{2}+X_{3}^{2}=1$. this corresponds
to the quantity $(10x_{1}^{2}+x_{3}^{2})$ diverging, $z=3x_{1}$ and all other
variables vanishing. All the fixed points in this family are saddles. This
entire family resides on the invariant submanifold $\mathcal{S}_{\mu}$. Just
like $\mathcal{O}_{\pm}$, the actual spacetime geometry corresponding to all
the fixed points within this family is Minkowski.

  3. 3.

A 4-parameter family of fixed points $\mathcal{B}^{(4)}$ given by
$X_{3}=0,\,X_{1}^{2}+X_{2}^{2}+X_{4}^{2}+Y_{2}^{2}+Z^{2}=1$ which correspond
to $x_{3}=0$ and the quantity
$(x_{1}^{2}+x_{2}^{2}+x_{4}^{2}+y_{2}^{2}+z^{2})$ diverging. All the fixed
points belonging to this family are also saddles. Different fixed points of
this family represent different spacetime geometries, but they are all
isotropic ($x_{3}=0$).

The fixed points at infinity of the phase space cannot represent singular
situations. Such fixed points are characterized by some expansion normalized
dynamical variables $x_{i}$ going to infinity. Therefore, either
$|\theta|\rightarrow 0$ or the covariant quantity at the numerator of $x_{i}$
diverges faster than $\theta$.

##  7 Lemaítre-Tolman-Bondi (LTB) spacetime

The spherically symmetric Lemaítre-Tolman-Bondi (LTB) spacetime, which belongs
to LRS-II class, is of special interest because of its relevance in the study
of gravitational collapse [49, 50, 51, 52, 53, 54] and in attempts to explain
late-time cosmological observations without invoking any dark energy component
[16, 17, 20, 19, 18]. Different aspects of LTB spacetime has been studied
before using the 1+1+2 formalism [55, 56]. This section is devoted to
describing the LTB spacetime geometry in light of the dynamical system
formulation developed in the previous sections.

In terms of 3D spherical polar coordinates LTB metric is written as

$ds^{2}=-dt^{2}+\frac{R^{\prime
2}(t,r)}{1-\mathcal{K}(r)}dr^{2}+R^{2}(t,r)d\Omega^{2}.$ (7.1)

For this metric the normalized timelike vector $u^{\mu}$ and the normalized
spacelike vector $e^{\mu}$ are, respectively,

$u^{\mu}=(1,0,0,0),\qquad
e^{\mu}=\left(0,\frac{\sqrt{1-\mathcal{K}(r)}}{R^{\prime}},0,0\right).$ (7.2)

For this particular choice of vectors, the covariantly defined time and radial
derivatives for a scalar quantity $\psi$ become

$\dot{\psi}=\partial_{t}\psi,\qquad\hat{\psi}=\frac{\sqrt{1-\mathcal{K}(r)}}{R^{\prime}}\partial_{r}\psi.$
(7.3)

The coordinate dependent expressions of the relevant non-vanishing geometric
scalars are as follows 777For the metric functions $\dot{(...)}$ means
derivative w.r.t. the time coordinate $t$, not to be confused with the
covariant time derivatives of the covariant quantities. $(...)^{\prime}$ means
derivative w.r.t. $r$.

$\displaystyle\phi\equiv d_{\mu}e^{\mu}=\frac{2\sqrt{1-\mathcal{K}(r)}}{R}\,,$
(7.4a) $\displaystyle\theta\equiv
D_{\mu}u^{\mu}=\frac{\dot{R}^{\prime}}{R^{\prime}}+2\frac{\dot{R}}{R}\,,$
(7.4b)
$\displaystyle\Sigma\equiv\sigma_{\mu\nu}e^{\mu}e^{\nu}=D_{\langle\mu}u_{\nu\rangle}e^{\mu}e^{\nu}=e^{\mu}e^{\nu}\left(h^{\alpha}_{\,\,\,(\mu}h^{\beta}_{\,\,\,\nu)}-\frac{1}{3}h_{\mu\nu}h^{\alpha\beta}\right)D_{\alpha}u_{\beta}=\frac{2}{3}\left(\frac{\dot{R}^{\prime}}{R^{\prime}}-\frac{\dot{R}}{R}\right)\,,$
(7.4c)
$\displaystyle\mathcal{E}\equiv\frac{1}{6R}\left(\frac{\mathcal{K}^{\prime}}{R^{\prime}}-2\frac{\mathcal{K}}{R}\right)+\frac{1}{3}\left(\frac{\ddot{R}}{R}-\frac{\ddot{R^{\prime}}}{R^{\prime}}\right)-\frac{1}{3}\frac{\dot{R}}{R}\left(\frac{\dot{R}}{R}-\frac{\dot{R^{\prime}}}{R^{\prime}}\right)\,,$
(7.4d) ${}^{3}R=2\frac{(\mathcal{K}R)^{\prime}}{R^{2}R^{\prime}}\,,$ (7.4e)
$\displaystyle K=\frac{1}{R^{2}}\,.$ (7.4f)

The above expressions specifically allow us to interpret the solutions on the
invariant submanifold $\mathcal{S}_{\phi}$ and $\mathcal{S}_{K}$ in a much
clearer manner in the context of LTB. Note that one must have $K(r)<1$ in an
LTB spacetime. Then, both the conditions $\phi=0$ and $K=0$ represent
solutions at spatial infinity $R(t,r)\rightarrow\infty$. Consequently, for an
LTB spacetime, the evolutionary phases represented by the isolated finite
fixed points
$\mathcal{P}_{1},\,\mathcal{P}_{2},\,\mathcal{P}_{3},\,\mathcal{P}_{4}$ can be
achieved by a comoving observer only when the observer is at spatial infinity.

Note that, since LTB is spherically symmetric, its phase space is bounded by

$x_{K}>0\Leftrightarrow x_{2}-9x_{4}+\frac{9}{4}x_{1}^{2}-(1-x_{3})^{2}>0\,.$
(7.5)

This implies that the asymptotic fixed point $\mathcal{O}^{\pm}$ is outside
the range of the LTB phase space. The other families of asymptotic fixed
points $\mathcal{A}^{(1)}$ and $\mathcal{B}^{(4)}$ can be achieved within the
LTB phase space.

The existence of the line of attractors $\mathcal{L}_{1}$ suggests that an
expanding LTB cosmology will asymptotically approach Milne. This conclusion is
in agreement with those available in the literature [46, 47].

The original Birkhoff’s theorem states that a spherically symmetric vacuum
solution in GR is necessarily either a part of the extended Schwarzschild
manifold or Minkowski. When the Weyl curvature tensor is non-vanishing, it is
necessarily a part of the Schwarzschild spacetime. However, it can be either
the static Schwarzschild exterior spacetime ($\theta=\Sigma=0$) or the
spatially homogeneous Schwarzschild interior spacetime ($\phi=0$). In light of
this, we can conclude that the isolated fixed points
$\mathcal{P}_{1},\,\mathcal{P}_{3}$ both represent Schwarzschild interior
solutions.

For the sake of completeness, we list in table 2 the finite and asymptotic
fixed points of the LTB phase space.

Fixed point | $(x_{1},x_{2},x_{3},x_{4},y_{2},z)$ | Stability |  | Actual
---
spacetime geometry
$\mathcal{L}_{1}$ | $(\lambda,0,0,0,0,0)$ |  | Attractor for expanding LTB
---
Repeller for contracting LTB
| Minkowski
---
(Milne w.r.t. comoving observer)
$\mathcal{L}_{2}$ | $\left(\lambda,0,-\frac{1}{2},0,0,\frac{1}{2}\lambda\right)$ | Saddle | Minkowski
$\mathcal{P}_{1}$ | $\left(0,0,-1,-\frac{4}{9},0,0\right)$ |  | Repeller for expanding LTB
---
Attractor for contracting LTB
Scharzschild interior
$\mathcal{P}_{2}$ | $(0,1,0,0,0,0)$ | Saddle |  | Spatially flat
---
matter dominated FLRW
$\mathcal{P}_{3}$ | $\left(0,0,\frac{1}{4},-\frac{1}{16},0,0\right)$ | Saddle | Scharzschild interior
$\mathcal{P}_{4}$ | $(0,0,1,0,0,0)$ | Saddle | Minkowski
$\mathcal{A}^{(1)}$ |  |  $(10x_{1}^{2}+x_{3}^{2})$ diverging, $z=3x_{1}$,
---
$x_{2}=x_{4}=y_{2}=0$
Saddle | Minkowski
$\mathcal{B}^{(4)}$ |  |  $(x_{1}^{2}+x_{2}^{2}+x_{4}^{2}+y_{2}^{2}+z^{2})$ diverging,
---
$x_{3}=0$
Saddle | Minkowski

Table 2: Finite and asymptotic fixed points of the LTB phase space.

The conclusions drawn regarding the stability of different invariant
submanifolds of the LRS-II phase space, which we listed in section 4, were
independent of the Gaussian curvature. Therefore they hold in exactly the same
form also for LTB. For LTB, the invariant submanifold $\mathcal{S}_{\rm hom}$
contains FLRW solutions and the stability of this submanifold determines the
evolution of almost-FLRW LTB geometries. From the generic conclusions that we
could draw regarding the stability of $\mathcal{S}_{\rm hom}$ in section 4, we
can conclude that the homogenization of an almost-FLRW LTB depends crucially
on the matter density. An almost-FLRW expanding LTB will isotropize provided
the matter density is in the range
$0\leq\frac{\mu}{\theta^{2}}<\frac{4}{9}-\frac{\sqrt{10}}{9}$, whereas an
almost-FLRW contracting LTB will isotropize provided
$\frac{\mu}{\theta^{2}}>1+\frac{\sqrt{7}}{3}$. Qualitatively one can conclude
that a large matter density hinders the homogenization of an almost-FLRW
expanding LTB whereas assists in the homogenization of an almost-FLRW
contracting LTB.

##  8 Conclusion

In this work, we present a new dynamical system formulation for LRS-II
spacetimes utilizing the 1+1+2 covariant decomposition approach [33, 34]. Our
approach describes the dynamics of an LRS-II spacetime from the point of view
of a comoving observer. A crucial point that we have emphasized throughout the
text is that for a comoving observer, the relevant dynamical quantities are
not only the covariant geometric and thermodynamic quantities but also their
covariantly defined “radial derivatives”, which we interpret as a specific
kind of perturbation that preserves the local rotational symmetry but stops
the spacetime from being globally homogeneous. This realization helps us
obtain an autonomous dynamical system containing a set of first-order ordinary
differential equations involving only time derivatives and a set of purely
algebraic constraints. Another nice feature of our formulation is that the
dimensionless dynamical variables that span the LRS-II phase space are
expansion normalized and directly related to actual covariant geometrical and
thermodynamic quantities, making it easier to interpret a fixed point or an
invariant submanifold with respect to a comoving observer. Due to the
inhomogeneity of the system, or in other words due to the existence of the
non-vanishing radial perturbations, not all the comoving observers will go
through the same evolutionary phase simultaneously.

We identify several invariant submanifolds and fixed points in sections 4, 5
and 6. We found that homogeneous solutions contain an invariant submanifold
whose stability depends only on the quantity $\frac{\mu}{\theta^{2}}$ (or
equivalently $\frac{{}^{3}R}{\theta^{2}}$). For LTB, one can qualitatively say
that a large matter density hinders the homogenization of an almost-FLRW
expanding LTB whereas assists in the homogenization of an almost-FLRW
contracting LTB. Interestingly, we find that spatially flat solutions in
general do not constitute an invariant submanifold in the LRS-II phase space
because of a non-vanishing shear. Finding a solution corresponding to vacuum
fixed points becomes ambiguous due to not having a preferred choice for
$u^{\mu}$ (i.e. a preferred frame), something which is a precondition for the
applicability of the 1+3 formalism. Special care is taken to interpret the
solutions corresponding to the vacuum fixed points. We find a spatially flat
matter-dominated FLRW solution as a saddle fixed point, i.e. an intermediate
epoch of evolution, as one might expect in the presence of nonrelativistic
matter. We also recovered a future attractor which a comoving observer would
identify with Milne, consistent with earlier results that expanding LTB
solutions without the cosmological constant tends to Milne [46, 47].

Let us finish with an outlook of possible further research in this line. The
logical next step is to generalize the present analysis including a
cosmological constant term (e.g. as in the $\Lambda$LTB model [57]) and/or a
non-vanishing pressure. In this work, we got a decoupled set of evolution and
propagation equations because we chose to confine our attention to perfect
fluids. Also, for perfect fluids the worldlines of comoving observers are
geodesics, since the acceleration vector $a^{\mu}$ vanishes. However, neither
the aforementioned decoupling nor the vanishing of acceleration is a necessary
condition to obtain an autonomous dynamical system in our approach; they just
merely make the calculations simpler. Even if one considers a non-perfect
fluid one could still proceed following our approach and construct an
autonomous dynamical system that describes the dynamics from the point of view
of a comoving observer. It is interesting to explore a formal comparison
between the 1+1+2-based dynamical system formulation that we developed here
with the dynamical system based on the orthonormal frame approach [28, 29, 30]
and Sussman’s approach in terms of quasi-local variables [31, 32]; whether or
not there exists a one-to-one correspondence between fixed points and
invariant submanifolds. It might be interesting to revisit the question
whether the cosmic no-hair theorem (an inflating spacetime isotropizes
asymptotically) holds in the presence of inhomogeneity, a question which has
been addressed previously using the orthonormal frame formalism [30]. Apart
from cosmology, we expect that this formulation should also find vast
application in the study of gravitational collapse [22, 23, 24]. In this
context, it is interesting to see if we can somehow incorporate the
information regarding the apparent horizon in the phase space of a collapsing
LTB and subsequently address the cosmic censorship conjecture 888For a study
of cosmic censorship in the 1+1+2 language, see [58, 59].. All these
interesting lines of research are reserved as potential future projects.

## Acknowledgement

SC acknowledges funding support from the NSRF via the Program Management Unit
for Human Resources and Institutional Development, Research and Innovation
(Thailand) [grant number B01F650006]. PKSD acknowledges the First Rand Bank
(SA) for financial support.

##  Appendix A An alternate approach to phase space analysis for ${}^{3}R>0$

For non-negative spatial curvature ${}^{3}R>0$, it is possible to formulate a
dynamical system with an alternative choice of dynamical variables. The
formulation we present here has the merit that for some homogeneous spacetimes
it can lead to a compact phase space, allowing for a global phase space
analysis [60]. For LRS-II, as we will see below, they do not lead to a compact
phase space. Nonetheless, it is still interesting to explore what this
alternative formulation gives us and how it compares with the usual
formulation in terms of the expansion normalized dynamical variables that we
have considered in the main body of the paper.

Inspired by the compact phase space formalism used in [60], we define the
following quantities

$\displaystyle
D\equiv\sqrt{\frac{\theta^{2}}{9}+\frac{{}^{3}R}{6}}=\sqrt{\frac{\mu}{3}+\frac{\Sigma^{2}}{4}}\,,$
(A.1) $\displaystyle U_{1}\equiv\frac{\phi}{D},\quad
U_{2}\equiv\frac{\mu}{3D^{2}},\quad U_{3}\equiv\frac{\Sigma}{2D},\quad
U_{4}\equiv\frac{\mathcal{E}}{D^{2}},\quad U_{5}\equiv\frac{\theta}{3D},\quad
U_{6}\equiv\sqrt{\frac{{}^{3}R}{6}}\frac{1}{D}\,.$ (A.2)

The Friedmann equation implies following two constraints

$U_{5}^{2}+U_{6}^{2}=U_{2}+U_{3}^{2}=1,$ (A.3)

using which we choose to eliminate $U_{2}$ and $U_{6}$. At the tangent space
of the observer, two sets of dynamical quantities are relevant: the original
covariant quantities as well as their covariant radial variations, given by
the hat derivatives. Therefore, we also define the following variables

$V_{1}\equiv\frac{\hat{\phi}}{D^{2}},\quad
V_{2}\equiv\frac{\hat{\mu}}{3D^{3}},\quad
V_{3}\equiv\frac{\hat{\Sigma}}{2D^{2}},\quad
V_{4}\equiv\frac{\hat{\mathcal{E}}}{D^{3}},\quad
V_{5}\equiv\frac{\hat{\theta}}{3D^{2}}\,.$ (A.4)

Together, Eqs.(A.2) and (A.4) provide the full set of dynamical variables with
respect to which we can write a closed set of autonomous equations describing
the LRS-II dynamics.

Clearly, except for $U_{3}$ and $U_{5}$, no other dimensions in the phase
space is compactified; $-1<U_{3},U_{5}<1$. However, $U_{5}$ being bounded
makes it easier to characterize the singularity; a big bang
($\theta\rightarrow+\infty$) corresponds to $U_{5}\rightarrow 1$ and a big
crunch ($\theta\rightarrow-\infty$) corresponds to $U_{5}\rightarrow-1$. It is
to be kept in mind that although singularities are represented by specific
points in the phase space, in the local $u$-$e$ plane they are represented by
curves (singularity curve).

The propagation equations (2.4) imply the following constraint equations

$\displaystyle V_{1}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}U_{1}^{2}+2(U_{5}+2U_{3})(U_{5}-U_{3})-2U_{2}-U_{4}\,,$
(A.5a) $\displaystyle V_{3}-V_{5}$ $\displaystyle=$
$\displaystyle-\frac{3}{2}U_{1}U_{3}\,,$ (A.5b) $\displaystyle V_{4}-V_{2}$
$\displaystyle=$ $\displaystyle-\frac{3}{2}U_{1}U_{4}\,.$ (A.5c)

using which we choose to eliminate $V_{1},\,V_{3},\,V_{4}$.

The quantity $D$ satisfies the following relations

$\frac{\dot{D}}{D^{2}}=-\frac{3}{2}U_{5}-\frac{1}{2}U_{3}(2U_{3}^{2}+U_{5}U_{3}+U_{4}),$
(A.6)

We define the following time variable on the phase space $\tau_{c}$ as

$d\tau_{c}=Ddt\,.$ (A.7)

We are now in a position to write the dynamical system in terms of the
dynamical variables defined in Eqs.(A.2) and (A.4). Using the relation (A.6),
we can write the system as

$\displaystyle\frac{dU_{1}}{d\tau_{c}}=U_{1}\left[\frac{1}{2}U_{5}+U_{3}+\frac{1}{2}U_{3}(2U_{3}^{2}+U_{5}U_{3}+U_{4})\right]\,,$
(A.8a)
$\displaystyle\frac{dU_{3}}{d\tau_{c}}=-\frac{1}{2}U_{4}-U_{3}(U_{3}+2U_{5})+\frac{1}{2}U_{3}^{2}(2U_{3}^{2}+U_{5}U_{3}+U_{4})\,,$
(A.8b)
$\displaystyle\frac{dU_{4}}{d\tau_{c}}=-3U_{3}(1-U_{4}-U_{3}^{2})+U_{3}X_{4}(2U_{3}^{2}+U_{5}U_{3}+U_{4})\,,$
(A.8c)
$\displaystyle\frac{dU_{5}}{d\tau_{c}}=-\frac{1}{2}(1-U_{5}^{2})-\frac{3}{2}U_{3}^{2}+\frac{1}{2}U_{5}U_{3}(2U_{3}^{2}+U_{5}U_{3}+U_{4})\,,$
(A.8d)
$\displaystyle\frac{dV_{2}}{d\tau_{c}}=\frac{1}{2}V_{2}(U_{5}-4U_{3})-3V_{5}(1-U_{3}^{2})+\frac{3}{2}U_{3}V_{2}(2U_{3}^{2}+U_{5}U_{3}+U_{4})\,,$
(A.8e)
$\displaystyle\frac{dV_{5}}{d\tau_{c}}=6U_{1}U_{3}^{2}-\frac{1}{2}V_{2}-6U_{3}V_{5}+U_{3}V_{5}(2U_{3}^{2}+U_{5}U_{3}+U_{4})\,.$
(A.8f)

The above system presents only two fixed points:
$P_{\pm}\equiv\\{U_{1},U_{3},U_{4},U_{5},V_{2},V_{5}\\}=\\{0,0,0,\pm
1,0,0\\}$, which represents big-bang and big-crunch singularities. Both of
them are saddle fixed points. Both these fixed points have $U_{2}=1$ and
$U_{6}=0$, i.e. matter dominates the spacetime dynamics on the approach to the
singularities whereas the contributions from the 3-curvature asymptotically
vanishes. $P_{\pm}$ resides on the invariant submanifold $\mathcal{S}_{\rm
hom}$, which now is given by $\\{U_{3},U_{4},V_{2},V_{5}\\}=\\{0,0,0,0\\}$,
implying that these are isotropic singularities.

The mapping between the sets of dynamical variables $\\{U_{i},V_{i}\\}$ and
$\\{x_{i},y_{i}\\}$ are as follows

$\displaystyle x_{1}=\frac{U_{1}}{3U_{5}},\quad
x_{2}=\frac{U_{2}}{U_{5}^{2}},\quad x_{3}=\frac{U_{3}}{U_{5}},\quad
x_{4}\equiv\frac{U_{4}}{9U_{5}^{2}},$ (A.9) $\displaystyle
y_{1}=\frac{V_{1}}{9U_{5}^{2}},\quad y_{2}=\frac{V_{2}}{3U_{5}^{3}},\quad
y_{3}=\frac{V_{3}}{3U_{5}^{2}},\quad y_{4}=\frac{V_{4}}{27U_{5}^{3}},$
$\displaystyle z=\frac{V_{5}}{3U_{5}^{2}}.$

Substituting back the coordinates $\\{U_{i},V_{i}\\}$’s of $P_{\pm}$ in
Eq.(A.9), we get the coordinates
$\\{x_{1},x_{2},x_{3},x_{4},y_{2},z\\}=\\{0,1,0,0,0,0\\}$, which corresponds
to the spatially flat matter-dominated homogeneous fixed point
$\mathcal{P}_{2}$ tn Tab.1. In the LTB phase space this point represents the
spatially flat matter dominated FLRW fixed point (Tab.2). The result is
consistent since both $P_{\pm}$ and $\mathcal{P}_{2}$ are saddles. Therefore,
the fixed points $P_{\pm}$ represent the big bang and big crunch singularities
already ingrained in the matter dominated homogeneous solutions. The singular
nature becomes apparent in the phase space formulation presented in this
section. However, this phase space formulation fails to give us any of the
other fixed points.

## References

  * [1] Christos G. Tsagas, Anthony Challinor, and Roy Maartens.  Relativistic cosmology and large-scale structure.  Phys. Rept., 465:61–147, 2008.
  * [2] Steven Weinberg.  The Cosmological Constant Problem.  Rev. Mod. Phys., 61:1–23, 1989.
  * [3] Eleonora Di Valentino et al.  Snowmass2021 - Letter of interest cosmology intertwined II: The hubble constant tension.  Astropart. Phys., 131:102605, 2021.
  * [4] Eleonora Di Valentino et al.  Cosmology Intertwined III: $f\sigma_{8}$ and $S_{8}$.  Astropart. Phys., 131:102604, 2021.
  * [5] Edmund J. Copeland, M. Sami, and Shinji Tsujikawa.  Dynamics of dark energy.  Int. J. Mod. Phys. D, 15:1753–1936, 2006.
  * [6] Luca Amendola and Shinji Tsujikawa.  Dark Energy: Theory and Observations.  Cambridge University Press, 2010.
  * [7] E. Papantonopoulos.  Modifications of Einstein’s Theory of Gravity at Large Distances.  Lecture Notes in Physics. Springer International Publishing, 2014.
  * [8] Yashar Akrami et al.  Modified Gravity and Cosmology: An Update by the CANTATA Network.  Springer, 2021.
  * [9] Sebastian Bahamonde, Christian G. Böhmer, Sante Carloni, Edmund J. Copeland, Wei Fang, and Nicola Tamanini.  Dynamical systems applied to cosmology: dark energy and modified gravity.  Phys. Rept., 775-777:1–122, 2018.
  * [10] Yasusada Nambu and Masayuki Tanimoto.  Accelerating universe via spatial averaging.  7 2005.
  * [11] John W. Moffat.  Late-time inhomogeneity and acceleration without dark energy.  JCAP, 05:001, 2006.
  * [12] Pantelis S. Apostolopoulos, Nikolaos Brouzakis, Nikolaos Tetradis, and Eleftheria Tzavara.  Cosmological acceleration and gravitational collapse.  JCAP, 06:009, 2006.
  * [13] Reza Mansouri.  Structured frw universe leads to acceleration: a non-perturbative approach.  12 2005.
  * [14] R. Ali Vanderveld, Eanna E. Flanagan, and Ira Wasserman.  Mimicking dark energy with Lemaitre-Tolman-Bondi models: Weak central singularities and critical points.  Phys. Rev. D, 74:023506, 2006.
  * [15] Havard Alnes, Morad Amarzguioui, and Oyvind Gron.  An inhomogeneous alternative to dark energy?  Phys. Rev. D, 73:083519, 2006.
  * [16] David Garfinkle.  Inhomogeneous spacetimes as a dark energy model.  Class. Quant. Grav., 23:4811–4818, 2006.
  * [17] Hideo Iguchi, Takashi Nakamura, and Ken-ichi Nakao.  Is dark energy the only solution to the apparent acceleration of the present universe?  Prog. Theor. Phys., 108:809–818, 2002.
  * [18] Miguel Zumalacarregui, Juan Garcia-Bellido, and Pilar Ruiz-Lapuente.  Tension in the Void: Cosmic Rulers Strain Inhomogeneous Cosmologies.  JCAP, 10:009, 2012.
  * [19] Valerio Marra and Alessio Notari.  Observational constraints on inhomogeneous cosmological models without dark energy.  Class. Quant. Grav., 28:164004, 2011.
  * [20] Leandros Perivolaropoulos.  Large Scale Cosmological Anomalies and Inhomogeneous Dark Energy.  Galaxies, 2:22–61, 2014.
  * [21] Nazeem Mustapha, Charles Hellaby, and G. F. R. Ellis.  Large scale inhomogeneity versus source evolution: Can we distinguish them observationally?  Mon. Not. Roy. Astron. Soc., 292:817–830, 1997.
  * [22] Pankaj S. Joshi.  Gravitational collapse: The Story so far.  Pramana, 55:529–544, 2000.
  * [23] Pankaj S. Joshi, Daniele Malafarina, and Ramesh Narayan.  Equilibrium configurations from gravitational collapse.  Class. Quant. Grav., 28:235018, 2011.
  * [24] Daniele Malafarina.  A Brief Review of Relativistic Gravitational Collapse.  Astrophys. Space Sci. Libr., 440:169–197, 2016.
  * [25] Dynamical Systems in Cosmology.  Cambridge University Press, 1997.
  * [26] Alan A. Coley.  Dynamical systems in cosmology.  In Spanish Relativity Meeting (ERE 99), 9 1999.
  * [27] Henk van Elst and Claes Uggla.  General relativistic (1+3) orthonormal frame approach revisited.  Class. Quant. Grav., 14:2673–2695, 1997.
  * [28] Henk van Elst, Claes Uggla, and John Wainwright.  Dynamical systems approach to G(2) cosmology.  Class. Quant. Grav., 19:51–82, 2002.
  * [29] Claes Uggla, Henk van Elst, John Wainwright, and George F. R. Ellis.  The Past attractor in inhomogeneous cosmology.  Phys. Rev. D, 68:103502, 2003.
  * [30] Woei Chet Lim, Henk van Elst, Claes Uggla, and John Wainwright.  Asymptotic isotropization in inhomogeneous cosmology.  Phys. Rev. D, 69:103507, 2004.
  * [31] Roberto A. Sussman.  A Dynamical system approach to inhomogeneous dust solutions.  Class. Quant. Grav., 25:015012, 2008.
  * [32] Roberto A. Sussman and German Izquierdo.  A Dynamical systems study of the inhomogeneous Lambda-CDM model.  Class. Quant. Grav., 28:045006, 2011.
  * [33] Christopher A. Clarkson and Richard K. Barrett.  Covariant perturbations of Schwarzschild black holes.  Class. Quant. Grav., 20:3855–3884, 2003.
  * [34] Chris Clarkson.  A Covariant approach for perturbations of rotationally symmetric spacetimes.  Phys. Rev. D, 76:104034, 2007.
  * [35] J. Ehlers.  Contributions to the relativistic mechanics of continuous media.  Abh. Akad. Wiss. Lit. Mainz. Nat. Kl., 11:793–837, 1961.
  * [36] George F. R. Ellis and Henk van Elst.  Cosmological models: Cargese lectures 1998.  NATO Sci. Ser. C, 541:1–116, 1999.
  * [37] Apratim Ganguly, Radouane Gannouji, Rituparno Goswami, and Subharthi Ray.  Global structure of black holes via the dynamical system.  Class. Quant. Grav., 32(10):105006, 2015.
  * [38] Miguel Cruz, Apratim Ganguly, Radouane Gannouji, Genly Leon, and Emmanuel N. Saridakis.  Global structure of static spherically symmetric solutions surrounded by quintessence.  Class. Quant. Grav., 34(12):125014, 2017.
  * [39] J. M. Stewart and G. F. R. Ellis.  Solutions of Einstein’s equations for a fluid which exhibit local rotational symmetry.  J. Math. Phys., 9:1072–1082, 1968.
  * [40] Gerold Betschart and Chris A. Clarkson.  Scalar and electromagnetic perturbations on LRS class II space-times.  Class. Quant. Grav., 21:5587–5608, 2004.
  * [41] Anne Marie Nzioki, Rituparno Goswami, and Peter K. S. Dunsby.  Covariant perturbations of Schwarzschild black holes in f(R) gravity.  Int. J. Mod. Phys. D, 26(06):1750048, 2016.
  * [42] Anne Marie Nzioki, Peter K. S. Dunsby, Rituparno Goswami, and Sante Carloni.  A Geometrical Approach to Strong Gravitational Lensing in f(R) Gravity.  Phys. Rev. D, 83:024030, 2011.
  * [43] Anne Marie Nzioki, Sante Carloni, Rituparno Goswami, and Peter K. S. Dunsby.  A New framework for studying spherically symmetric static solutions in f(R) gravity.  Phys. Rev. D, 81:084028, 2010.
  * [44] Rituparno Goswami and George F R Ellis.  Almost Birkhoff Theorem in General Relativity.  Gen. Rel. Grav., 43:2157–2170, 2011.
  * [45] Sayuri Singh, George F. R. Ellis, Rituparno Goswami, and Sunil D. Maharaj.  New class of LRS spacetimes with simultaneous rotation and spatial twist.  Phys. Rev. D, 94(10):104040, 2016.
  * [46] A. A. Coley, W. C. Lim, and G. Leon.  Spherically symmetric cosmology: Resource paper.  3 2008.
  * [47] John Wainwright and Shawn Andrews.  The Dynamics of Lemaitre-Tolman cosmologies.  Class. Quant. Grav., 26:085017, 2009.
  * [48] Rituparno Goswami and George F. R. Ellis.  Birkhoff Theorem and Matter.  Gen. Rel. Grav., 44:2037–2050, 2012.
  * [49] P. S. Joshi and I. H. Dwivedi.  Naked singularities in spherically symmetric inhomogeneous Tolman-Bondi dust cloud collapse.  Phys. Rev. D, 47:5357–5369, 1993.
  * [50] Pankaj S. Joshi and Daniele Malafarina.  All black holes in Lemaître–Tolman–Bondi inhomogeneous dust collapse.  Class. Quant. Grav., 32(14):145004, 2015.
  * [51] Pankaj S. Joshi.  On the genericity of spacetime singularities.  Pramana, 69:119–136, 2007.
  * [52] Pankaj S. Joshi, Naresh Dadhich, and Roy Maartens.  Why do naked singularities form in gravitational collapse?  Phys. Rev. D, 65:101501, 2002.
  * [53] Pankaj S. Joshi, Rituparno Goswami, and Naresh Dadhich.  Why do naked singularities form in gravitational collapse? 2.  Phys. Rev. D, 70:087502, 2004.
  * [54] Pankaj S. Joshi and Rituparno Goswami.  Role of initial data in spherical collapse.  Phys. Rev. D, 69:064027, 2004.
  * [55] Peter Dunsby, Naureen Goheer, Bob Osano, and Jean-Philippe Uzan.  How close can an Inhomogeneous Universe mimic the Concordance Model?  JCAP, 06:017, 2010.
  * [56] James P. Zibin.  Scalar Perturbations on Lemaitre-Tolman-Bondi Spacetimes.  Phys. Rev. D, 78:043504, 2008.
  * [57] Peter Sundell and Iiro Vilja.  Inhomogeneity of the $\Lambda$LTB models.  1 2016.
  * [58] Aymen I. M. Hamid, Rituparno Goswami, and Sunil D. Maharaj.  Cosmic Censorship Conjecture revisited: Covariantly.  Class. Quant. Grav., 31:135010, 2014.
  * [59] Aymen I. M. Hamid, Rituparno Goswami, and Sunil D. Maharaj.  Cosmic censorship conjecture revisited: Covariantly.  In 14th Marcel Grossmann Meeting on Recent Developments in Theoretical and Experimental General Relativity, Astrophysics, and Relativistic Field Theories, volume 3, pages 2615–2620, 2017.
  * [60] Deon M. Solomons, Peter Dunsby, and George Ellis.  Bounce behaviour in Kantowski-Sachs and Bianchi cosmologies.  Class. Quant. Grav., 23:6585–6597, 2006.