# Bouncing cosmology from nonlinear dark energy with two cosmological
constants

Molly<EMAIL_ADDRESS>and Marco
<EMAIL_ADDRESS>aInstitute of Cosmology & Gravitation,
University of Portsmouth, Dennis Sciama Building, Burnaby Road, Portsmouth,
PO1 3FX, United Kingdom
bINFN Sezione di Trieste, Via Valerio 2, 34127 Trieste, Italy

###### Abstract

We explore the dynamics of FLRW cosmologies which consist of dark matter,
radiation and dark energy with a quadratic equation of state. Standard
cosmological singularities arise due to energy conditions which are violated
by dark energy, therefore we focus our analysis on non-singular bouncing and
cyclic cosmologies, in particular focusing on the possibility of closed FLRW
models always having a bounce for any initial conditions. We analyse the range
of dynamical behaviour admitted by the system, and find a class of closed
models that admit a non-singular bounce, with early- and late-time accelerated
expansion connected by a decelerating phase. In all cases, we find the
bouncing models are only relevant when dark matter and radiation appear at a
certain energy scale, and so require a period such as reheating. We then
investigate imposing an upper bound on the dark matter and radiation, such
that their energy densities cannot become infinite. In this case, we find that
all closed models bounce, and a class of models exist with early- and late-
time acceleration, connected by a decelerating phase. We also consider
parameter values for the dark energy component, such that the discrepancy
between the observed value of $\Lambda$ and the theoretical estimates of the
contributions to the effective cosmological constant expected from quantum
field theory would be explained. However, we find that the class of models
left does not allow for an early- and late-time accelerated expansion,
connected by a decelerating period where large-scale structure could form.
Nonetheless, our qualitative analysis serves as a basis for the construction
of more realistic models with realistic quantitative behaviour.

##  I Introduction

The $\Lambda$-cold-dark-matter ($\Lambda$CDM) model has provided a successful
framework to describe the history of the Universe, and is consistent with a
wealth of observational data Aghanim et al. (2018a, b). Despite its
robustness, there are problems which require addressing. The issue of the
inevitability of a singularity as the origin of our Universe has been debated
for many years. Assuming the strong energy condition holds, singularities
appear to be unavoidable Penrose (1965); Hawking (1965, 1966a, 1966b, 1967);
Hawking and Penrose (1970); Hawking and Ellis (1973); Wald (1984), however
their current interpretation is that they represent points where General
Relativity breaks down Joshi (1993); Hawking and Israel (1979); Gibbons et al.
(2003); Ashtekar et al. (2015). Broadly speaking, there are two ways to tackle
this problem: solve the singularity problem by developing a modified theory of
gravity to unify General Relativity with this high energy regime, or avoid the
singularity by breaking the standard energy conditions in an alternative
origin story of the Universe in General Relativity. For the latter,
nonsingular bouncing models have been proposed as a way to evade a singularity
Gasperini and Veneziano (1993); Khoury et al. (2001); Kallosh et al. (2001);
Steinhardt and Turok (2002); Novello and Bergliaffa (2008); Wands (2008);
Maier et al. (2012, 2013); Battefeld and Peter (2014); Cai et al. (2017);
Brandenberger and Peter (2017); Cai and Piao (2017); Ijjas and Steinhardt
(2018, 2019); Ganguly and Bruni (2019); Ilyas et al. (2020); Battista (2021);
Zhu et al. (2021).

Another problem facing $\Lambda$CDM is that of $\Lambda$ itself. Observations
have provided strong evidence for the accelerated expansion of our Universe
Riess et al. (1998); Perlmutter et al. (1998), which requires a component that
violates the strong energy condition Hawking and Ellis (1973), known as dark
energy. The cosmological constant $\Lambda$ is the simplest explanation we
have for dark energy, however the observed value of $\Lambda$ is at odds by up
to 120 orders of magnitude with theoretical estimates of the contributions to
the effective cosmological constant expected from quantum field theory (QFT)
Weinberg (1989); Straumann (1999); Ellis et al. (2012).

Taking into account a dynamic dark energy component beyond $\Lambda$, it is
worth considering whether a bounce can be produced. A non-interacting vacuum
is equivalent to $\Lambda$ in General Relativity Lemaître (1934), so a simple
extension is to consider a vacuum that interacts with other components,
thereby becoming dynamical Wands et al. (2012). It has been shown that non-
singular bounces can arise for both a linear and non-linear interaction of a
vacuum with a perfect fluid in Friedmann-Lemaître-Robertson-Walker (FLRW)
cosmologies Bruni et al. (2022).

For non-interacting models, dark energy represented by a single fluid with a
quadratic equation of state (EoS) in FLRW can admit nonsingular bounces with
the right combination of parameters Ananda and Bruni (2006a). In addition, a
quadratic EoS finds motivations from e.g. brane-world models Shiromizu et al.
(2000); Bridgman et al. (2002); Gergely and Maartens (2002); Langlois (2003);
Maartens and Koyama (2010), loop quantum cosmology Vandersloot (2005); Parisi
et al. (2007) and k-essence Scherrer (2004); Giannakis and Hu (2005) (see also
Ananda and Bruni (2006a, b) and references therein). In general, a nonlinear
EoS can allow for the existence of effective cosmological constants, appearing
as asymptotic states of the dynamical dark energy. Although simple, analysing
a quadratic EoS serves well to illustrate the general qualitative behaviour of
a system that can admit two effective cosmological constants, a high energy
one acting as a repellor and a low energy one acting as a future attractor for
the expanding models. Therefore, it is useful to understand the global
dynamics as it is only dependent on the existence of the effective
cosmological constants and not on the specific model used to produce them.

In this paper, we analyse a subset of the quadratic EoS from Ananda and Bruni
(2006a, b) which admit two effective cosmological constants and produce a
bounce for models with positively curved spatial sections, i.e. closed FLRW
models. This dark energy only violates the strong energy condition, and not
the null energy condition, hence a bounce can only occur when the curvature is
positive. Here, we extend the scenario in Ananda and Bruni (2006a, b) to
include non-interacting dark matter and radiation in order to understand the
effect these components have on the bounce. In particular, we investigate
whether all closed models will bounce when matter and radiation are included,
thereby providing a scenario where the bounce is generic, i.e. does not depend
on initial conditions. Standard radiation and dark matter diverge more quickly
at high-enough energies than the curvature and dark energy considered here. We
find that we obtain bouncing models for certain values of parameters, but at
high energies, if dark matter and radiation become dominant, some models
become singular. We therefore consider models where dark matter and radiation
are produced at a certain energy scale after the bounce, by imposing an upper
bound on their energy densities. In this case, we find that all closed models
bounce.

We would also like to understand whether such a model could alleviate the
cosmological constant problem, and unite theoretical estimates from QFT with
observational constraints on $\Lambda$. We restrict our analysis such that the
dark energy has an effective EoS parameter satisfying $w>-1$, so that it
decays during the expansion between the high energy and low energy
cosmological constants, even allowing a decelerated phase in general, and do
not consider phantom dark energy ($w<-1$). The analysis for this paper is
available through GitHub, provided the reader has a Mathematica license
Burkmar (2023).

The paper is organized as follows. In Section II we present the system of
equations we are analysing in terms of dimensionless variables as well as in
compactified form. In Section III we show the sub-manifolds of the system,
where the dimensions of the dynamics is reduced and we explore the stability
of the fixed points. The sub-manifolds of the system include cases analogous
to $\Lambda$CDM, and ones in which only dark energy is present which highlight
the dynamics before including dark matter and radiation. In Section IV we
study the dynamics of the full system, where both dark matter and radiation
are present in standard form, i.e. with EoS parameters $w_{m}=0$ and
$w_{r}=1/3$. We show that a class of closed models exist where bounces are
obtained with early- and late-time accelerated expansion phases, connected by
a decelerating phase. The caveat with this full system is that not all closed
models bounce, and above a certain energy scale matter and radiation become
dominant, and, depending on the initial conditions, a subset of closed models
evolve between two singularities. Therefore, these models require a post-
bounce process such as reheating, when matter and radiation are created. In
this spirit, in Section V we re-consider the equations for dark matter and
radiation, introducing scale-dependent EoS’s that limit their energy densities
to an upper bound, i.e. an energy scale at which they appear, such that at
lower energies we still have $w_{m}\rightarrow 0$ and $w_{r}\rightarrow 1/3$.
We present the phase spaces for the new system of equations, showing that in
general all closed models bounce. We find that it is possible for these models
to evolve with early- and late-time accelerated expansion separated by a
deceleration period, however the cosmological constant problem cannot be
solved simultaneously. We conclude in Section VI. In this paper, we employ
natural units such that $8\pi G=c=1$.

##  II Cosmological Dynamical System

###  II.1 Cosmology with nonlinear EoS dark energy

FLRW models with no cosmological constant term in Einstein’s field equations
evolve according to a dynamical system, consisting of the energy conservation
equations for each component and the Raychaudhuri equation. We consider the
dynamics of a universe filled with dark energy, pressureless dark matter and
radiation, which evolve according to

$\dot{\rho}_{x}=-3H(\rho_{x}+P_{x})$ (1)

$\dot{\rho}_{m}=-3H\rho_{m}$ (2)

$\dot{\rho}_{r}=-4H\rho_{r}\,,$ (3)

where $\rho_{x}$ and $P_{x}$ are the energy density and isotropic pressure of
dark energy, $\rho_{m}$ is the energy density of dark matter and $\rho_{r}$ is
the energy density of radiation; overdots are derivatives with respect to time
$t$. $H$ is the Hubble expansion function, and the Raychaudhuri equation
describing it’s evolution is

$\dot{H}=-H^{2}-\frac{1}{6}(\rho_{x}+3P_{x}+\rho_{m}+2\rho_{r})\,.$ (4)

We assume that the dark energy EoS is barotropic, $P_{x}=P_{x}(\rho_{x})$,
with $P_{x}(\rho_{x})\geq-\rho_{x}$, such that Eq. (1) has two fixed points
($\dot{\rho}_{x}=0$), corresponding to two effective cosmological constants:

$P_{x}(\rho_{\Lambda})=-\rho_{\Lambda},~{}~{}~{}P_{x}(\rho_{*})=-\rho_{*}.$
(5)

It follows from the assumption $P_{x}(\rho_{x})\geq-\rho_{x}$ that the dark
energy density decreases/increases during a period of expansion/contraction.
Then, during expansion ($H>0$), we take $\rho_{\Lambda}$ to represent the
attractor at low energy close to the dark energy density we observe today, and
$\rho_{*}$ the repellor at high energy; their roles are inverted during
contraction. We assume $\rho_{*}$ is positive and is an energy scale between
the Planck scale and that typical of inflation, to ensure the evolution is
always classical, but does not interfere with particle production at lower
energies.

In this paper we assume the same quadratic EoS used in Ananda and Bruni
(2006a, b), which is the simplest non-linear case to study and, as we said in
the introduction (see also Ananda and Bruni (2006a, b)), has other motivations
from various scenarios Shiromizu et al. (2000); Bridgman et al. (2002);
Gergely and Maartens (2002); Langlois (2003); Maartens and Koyama (2010);
Vandersloot (2005); Parisi et al. (2007); Scherrer (2004); Giannakis and Hu
(2005). With respect to Ananda and Bruni (2006a, b), we restrict parameters
such that the linear term is always positive, $\alpha>0$, the quadratic term
is always negative, $\beta=-1/\rho_{*}$, and the constant pressure term is
negative, $P_{0}=-\rho_{\Lambda}$. With this choice of parameters our EoS is

$P_{x}=-\rho_{\Lambda}+\frac{\rho_{\Lambda}}{\rho_{*}}\rho_{x}-\frac{\rho_{x}^{2}}{\rho_{*}}\,.$
(6)

This secures the two effective cosmological constants in Eq. (5), as now the
energy conservation equation of the dark energy (1) can be written as

$\dot{\rho_{x}}=-3H(\rho_{x}-\rho_{\Lambda})\left(1-\frac{\rho_{x}}{\rho_{*}}\right).$
(7)

We note however that the qualitative dynamics that follows from the the
specific EoS (6) is going to be representative of the dynamics for any EoS
satisfying the condition $P_{x}(\rho_{x})+\rho_{x}>0$ for
$\rho_{\Lambda}<\rho_{x}<\rho_{*}$, i.e. between the two fixed points (5)
satisfying $P_{x}(\rho_{x})=-\rho_{x}$, as these conditions are enough to
describe any monotonically decreasing $\rho_{x}$ between $\rho_{*}$ and
$\rho_{\Lambda}$ during expansion, see e.g. the EoS in Hogg and Bruni (2021).

The effective EoS parameter $w_{x}=P_{x}/\rho_{x}$ for the dark energy is as

$w_{x}=-\frac{\rho_{\Lambda}}{\rho_{x}}+\frac{\rho_{\Lambda}}{\rho_{*}}-\frac{\rho_{x}}{\rho_{*}}\,,$
(8)

and the EoS parameters for dark matter and radiation are $w_{m}=0$ and
$w_{r}=1/3$, respectively. This system admits a first integral, the Friedmann
equation, which is written in terms of each energy density as

$H^{2}=\frac{\rho_{m}}{3}+\frac{\rho_{r}}{3}+\frac{\rho_{x}}{3}-\frac{k}{a^{2}}\,,$
(9)

where $k$ is the curvature and $a$ is the cosmic scale factor, connected to
the Hubble expansion function through the expression $H=\dot{a}/a$. We set
$a_{0}=1$ today, therefore $k$ is an arbitrary constant which is positive,
negative or zero for closed, open and flat models, respectively.

###  II.2 Dimensionless Variables

Examining the equations above, a dimensional analysis suggests that the
dynamics really only depends on a single parameter, the dimensionless ratio of
the two effective cosmological constants, if we use dimensionless variables.
Following Ananda and Bruni (2006a, b), we define these variables as:

$x=\frac{\rho_{x}}{\rho_{*}}~{}~{}y=\frac{H}{\sqrt{\rho_{*}}}~{}~{}z=\frac{\rho_{m}}{\rho_{*}}~{}~{}r=\frac{\rho_{r}}{\rho_{*}}~{}~{}\mathcal{R}=\frac{\rho_{\Lambda}}{\rho_{*}}~{}~{}\eta=\sqrt{\rho_{*}}t.$
(10)

The variable $\mathcal{R}$ is the ratio of the low energy effective
cosmological constant $\rho_{\Lambda}$ to the high energy effective
cosmological constant $\rho_{*}$ and takes values in the range $(0,1)$, $x$ is
the normalised dark energy density, varying between the two dimensionless
effective cosmological constants in the range $[\mathcal{R},1]$. The variable
$y$ is the normalised Hubble parameter, $z$ the normalised dark matter energy
density, $r$ the normalised radiation energy density and $\eta$ the normalised
time variable. We consider the region of phase space where the energy
densities for matter and radiation are always positive, i.e. $z,r>0$.
Equations (2), (3) and (7) then become

$z^{\prime}=-3yz$ (11)

$r^{\prime}=-4yr$ (12)

$x^{\prime}=-3y(x-\mathcal{R})(1-x)\,,$ (13)

where the primes indicate differentiation with respect to $\eta$. The
Raychaudhuri equation (4) can now be expressed as

$y^{\prime}=-y^{2}-\frac{1}{6}[z+2r-3\mathcal{R}+(1+3\mathcal{R})x-3x^{2}]\,,$
(14)

and the Friedmann equation (9) becomes

$y^{2}=\frac{x}{3}+\frac{z}{3}+\frac{r}{3}-\frac{k}{\rho_{*}a^{2}}.$ (15)

The effective EoS parameter for the dark energy (8) in dimensionless variables
becomes

$w_{x}=-\frac{\mathcal{R}}{x}+\mathcal{R}-x,$ (16)

therefore, the necessary condition for acceleration is

$-\frac{\mathcal{R}}{x}+\mathcal{R}-x<-1/3.$ (17)

We can also express the evolution of the cosmic scale factor $a$ using our
dimensionless variables

$a^{\prime}=ay\,.$ (18)

The Friedmann equation (15) can then be rearranged for a kinetic term
$a^{\prime 2}/2$, a potential $U$ and total energy $E$

$\frac{a^{\prime 2}}{2}+U=E\,,$ (19)

where $U$ is given by

$U=-\frac{a^{2}}{6}(x+z+r)$ (20)

and $E$ by

$E=-\frac{k}{2\rho_{*}}\,.$ (21)

$E$ is zero, positive and negative for flat, open and closed models,
respectively.

From Eq. (18) and Eq.s (11) and (12), we get the standard behaviour for matter
and radiation, $z(a)\sim a^{-3}$ and $r(a)\sim a^{-4}$, respectively. For $x$
(13), we obtain

$x(a)=\frac{a^{-3(1-\mathcal{R})}+c_{a}\mathcal{R}}{a^{-3(1-\mathcal{R})}+c_{a}}\,,$
(22)

where the constant of integration is

$c_{a}=\frac{1-x_{0}}{x_{0}-\mathcal{R}}\,.$ (23)

In the following, however, we will express $z$, $r$ and $a$ as functions of
$x$.

###  II.3 Compactified Variables

In our system, $y$ takes values in the range $(-\infty,\infty)$, and $z$ and
$r$ in the range $[0,\infty)$. $x$ is already limited to the range
$[\mathcal{R},1]$. In order to analyse the full dynamics, it is desirable to
deal with a compact phase space. Therefore, we compactify the $y$, $z$ and $r$
variables in the following way:

$Y=\frac{y}{\sqrt{1+y^{2}}}$ (24)

$Z=\frac{z}{1+z}$ (25)

$R=\frac{r}{1+r}\,,$ (26)

so that $Y$ takes values in the range $[-1,1]$ and $R$ and $Z$ in the range
$[0,1]$. Using these variables, Eq. (11) - (14) become

$Z^{\prime}=\frac{-3YZ(1-Z)}{\sqrt{1-Y^{2}}}$ (27)

$R^{\prime}=\frac{-4YR(1-R)}{\sqrt{1-Y^{2}}}$ (28)

$x^{\prime}=\frac{-3Y(x-\mathcal{R})(1-x)}{\sqrt{1-Y^{2}}}$ (29)

$Y^{\prime}=-Y^{2}\sqrt{1-Y^{2}}-\frac{(1-Y^{2})^{\frac{3}{2}}}{6}\\\
\times\left[\frac{Z}{1-Z}+\frac{2R}{1-R}+x(1+3\mathcal{R})-3\mathcal{R}-3x^{2}\right]\,,$
(30)

and the Friedmann equation (15) can be expressed as

$\frac{Y^{2}}{1-Y^{2}}=\frac{x}{3}+\frac{Z}{3(1-Z)}+\frac{R}{3(1-R)}-\frac{k}{\rho_{*}a^{2}}\,.$
(31)

Now, setting $a_{0}=1$, we find

$\frac{k}{\rho_{*}}=\frac{x_{0}}{3}+\frac{Z_{0}}{3(1-Z_{0})}+\frac{R_{0}}{3(1-R_{0})}-\frac{Y_{0}^{2}}{1-Y_{0}^{2}}\,.$
(32)

Finally, the potential $U$ (20) in terms of compact variables is given by

$U=-\frac{a^{2}}{6}\left(x+\frac{Z}{1-Z}+\frac{R}{1-R}\right)\,.$ (33)

##  III Sub-manifolds of the System

To analyse our autonomous dynamical system, we first need to find the sub-
manifolds and the fixed points, and determine their linear stability character
Arrowsmith and Place (1982, 1992); Strogatz (2015). The fixed points
$u_{j}^{*}$ satisfy $f_{i}(u_{j}^{*})=0$, where $f_{i}$ are the first-order
derivatives of the independent variables $u_{j}$ with respect to time. We then
linearize about each fixed point, first finding the Jacobian matrix, which has
the form

$\textbf{M}_{ij}=\frac{\partial f_{i}}{\partial u_{j}}.$ (34)

Table 1:  The sub-manifolds of the system with their fixed points, eigenvalues and stability character. $E$ denotes an Einstein universe, $dS_{\pm}$ an expanding (+) or contracting (-) de-Sitter universe, and $S_{\pm}$ a singularity with infinite expansion (+) or contraction (-). The fixed points at $Y=\pm 1$ representing infinities do not admit a linearization, and therefore eigenvalues cannot be found (denoted by N/A in the Eigenvalues column). We provide details of how we find their stability character in the text. Sub-manifold | Fixed Points | Name | Eigenvalues | Stability Character
---|---|---|---|---
$x=1$, $Z=0$ |  $Y=\pm\frac{1}{2}$, $R=0$ | $dS_{1\pm}$ | $\begin{pmatrix}\mp\frac{4}{\sqrt{3}}\\\ \mp\frac{2}{\sqrt{3}}\\\ \end{pmatrix}$ | Attractor/Repellor
|  $Y=0$, $R=\frac{1}{2}$ | $E$ | $\begin{pmatrix}\frac{2}{\sqrt{3}}\\\ -\frac{2}{\sqrt{3}}\end{pmatrix}$ | Saddle
|  $Y=\pm 1$, $R=0$ | $dS_{2\pm}$ | N/A | Saddle
|  $Y=\pm 1$, $R=1$ | $S_{\pm}$ | N/A | Repellor/Attractor
$x=1$, $R=0$ |  $Y=\pm\frac{1}{2}$, $Z=0$ | $dS_{1\pm}$ | $\begin{pmatrix}\mp\sqrt{3}\\\ \mp\frac{2}{\sqrt{3}}\\\ \end{pmatrix}$ | Attractor/Repellor
|  $Y=0$, $Z=\frac{2}{3}$ | $E$ | $\begin{pmatrix}1\\\ -1\end{pmatrix}$ | Saddle
|  $Y=\pm 1$, $Z=0$ | $dS_{2\pm}$ | N/A | Saddle
|  $Y=\pm 1$, $Z=1$ | $S_{\pm}$ | N/A | Repellor/Attractor
$x=1$ |  $Y=\pm\frac{1}{2}$, $Z=0$, $R=0$ | $dS_{1\pm}$ | $\begin{pmatrix}\mp\frac{4}{\sqrt{3}}\\\ \mp\sqrt{3}\\\ \mp\frac{2}{\sqrt{3}}\\\ \end{pmatrix}$ | Attractor/Repellor
|  $Y=0$, $R=\frac{-2+3Z}{-4+5Z}$ | $E$ | $\begin{pmatrix}0\\\ \sqrt{\frac{8-9Z}{6(1-Z)}}\\\ -\sqrt{\frac{8-9Z}{6(1-Z)}}\end{pmatrix}$ | Saddle
|  $Y=\pm 1$, $Z=0$, $R=0$ | $dS_{2\pm}$ | N/A | Saddle
|  $Y=\pm 1$, $Z=1$, $R=1$ | $S_{\pm}$ | N/A | Repellor/Attractor
$x=\mathcal{R}$, $Z=0$ |  $Y=\pm\sqrt{\frac{\mathcal{R}}{3+\mathcal{R}}}$, $R=0$ | $dS_{1\pm}$ | $\begin{pmatrix}\mp\frac{4\sqrt{\mathcal{R}}}{\sqrt{3}}\\\ \mp\frac{2\sqrt{\mathcal{R}}}{\sqrt{3}}\\\ \end{pmatrix}$ | Attractor/Repellor
|  $Y=0$, $R=\frac{\mathcal{R}}{1+\mathcal{R}}$ | $E$ | $\begin{pmatrix}\frac{2\sqrt{\mathcal{R}}}{3}\\\ -\frac{2\sqrt{\mathcal{R}}}{3}\end{pmatrix}$ | Saddle
|  $Y=\pm 1$, $R=0$ | $dS_{2\pm}$ | N/A | Saddle
|  $Y=\pm 1$, $R=1$ | $S_{\pm}$ | N/A | Repellor/Attractor
$x=\mathcal{R}$, $R=0$ |  $Y=\pm\sqrt{\frac{\mathcal{R}}{3+\mathcal{R}}}$, $Z=0$ | $dS_{1\pm}$ | $\begin{pmatrix}\mp\sqrt{3\mathcal{R}}\\\ \mp\frac{2\sqrt{\mathcal{R}}}{\sqrt{3}}\end{pmatrix}$ | Attractor/Repellor
|  $Y=0$, $Z=\frac{2\mathcal{R}}{1+2\mathcal{R}}$ | $E$ | $\begin{pmatrix}\sqrt{\mathcal{R}}\\\ -\sqrt{\mathcal{R}}\end{pmatrix}$ | Saddle
|  $Y=\pm 1$, $Z=0$ | $dS_{2\pm}$ | N/A | Saddle
|  $Y=\pm 1$, $Z=1$ | $S_{\pm}$ | N/A | Repellor/Attractor
$x=\mathcal{R}$ |  $Y=\pm\sqrt{\frac{\mathcal{R}}{3+\mathcal{R}}}$, $Z=0$, $R=0$ | $dS_{1\pm}$ | $\begin{pmatrix}\mp\frac{4\sqrt{\mathcal{R}}}{\sqrt{3}}\\\ \mp\sqrt{3\mathcal{R}}\\\ \mp\frac{2\sqrt{\mathcal{R}}}{\sqrt{3}}\end{pmatrix}$ | Attractor/Repellor
|  $Y=0$, $R=\frac{Z-2\mathcal{R}+2Z\mathcal{R}}{-2+3Z-2\mathcal{R}+2Z\mathcal{R}}$ | $E$ | $\begin{pmatrix}0\\\ \sqrt{\frac{8\mathcal{R}-Z(1+8\mathcal{R})}{6(1-Z)}}\\\ -\sqrt{\frac{8\mathcal{R}-Z(1+8\mathcal{R})}{6(1-Z)}}\end{pmatrix}$ | Saddle
|  $Y=\pm 1$, $Z=0$, $R=0$ | $dS_{2\pm}$ | N/A | Saddle
|  $Y=\pm 1$, $Z=1$, $R=1$ | $S_{\pm}$ | N/A | Repellor/Attractor
$Z=0$, $R=0$ |  $Y=\pm\sqrt{\frac{\mathcal{R}}{3+\mathcal{R}}}$, $x=\mathcal{R}$ | $dS_{1\pm}$ | $\begin{pmatrix}\mp\frac{2\sqrt{\mathcal{R}}}{\sqrt{3}}\\\ \mp\sqrt{3\mathcal{R}}(1-\mathcal{R})\end{pmatrix}$ | Attractor/Repellor
|  $Y=\pm\frac{1}{2}$, $x=1$ | $dS_{2\pm}$ | $\begin{pmatrix}\mp\frac{2}{\sqrt{3}}\\\ \pm\sqrt{3}(1-\mathcal{R})\end{pmatrix}$ | Saddle
|  $Y=0$, $x=\frac{1}{6}(1+3\mathcal{R}\pm\sqrt{1-30\mathcal{R}+9\mathcal{R}^{2}})$ | $E$ | $\begin{pmatrix}\sqrt{\frac{-1+30\mathcal{R}-9\mathcal{R}^{2}\mp(1+3\mathcal{R})\sqrt{1-30\mathcal{R}+9\mathcal{R}^{2}}}{18}}\\\ -\sqrt{\frac{-1+30\mathcal{R}-9\mathcal{R}^{2}\mp(1+3\mathcal{R})\sqrt{1-30\mathcal{R}+9\mathcal{R}^{2}}}{18}}\end{pmatrix}$ | Saddle, Cusp or Centre
|  $Y=\pm 1$, $x=\mathcal{R}$ | $dS_{3\pm}$ | N/A | Saddle
|  $Y=\pm 1$, $x=1$ | $dS_{4\pm}$ | N/A | Repellor/Attractor
$Z=0$ |  $Y=\pm\sqrt{\frac{\mathcal{R}}{3+\mathcal{R}}}$, $x=\mathcal{R}$, $R=0$ | $dS_{1\pm}$ | $\begin{pmatrix}\mp\frac{4\sqrt{R}}{\sqrt{3}}\\\ \mp\frac{2\sqrt{\mathcal{R}}}{\sqrt{3}}\\\ \mp\sqrt{3\mathcal{R}}(1-\mathcal{R})\end{pmatrix}$ | Attractor/Repellor
|  $Y=0$, $R=\frac{-x+3x^{2}+3\mathcal{R}-3x\mathcal{R}}{2-x+3x^{2}+3\mathcal{R}-3x\mathcal{R}}$ | $E$ | $\begin{pmatrix}0\\\ \sqrt{\frac{x(9\mathcal{R}^{2}+18\mathcal{R}-1)-x^{2}(27\mathcal{R}+9)+18x^{3}+9\mathcal{R}(1-\mathcal{R})}{6}}\\\ -\sqrt{\frac{x(9\mathcal{R}^{2}+18\mathcal{R}-1)-x^{2}(27\mathcal{R}+9)+18x^{3}+9\mathcal{R}(1-\mathcal{R})}{6}}\end{pmatrix}$ | Saddle
|  $Y=\pm 1$, $x=\mathcal{R}$, $R=0$ | $dS_{2\pm}$ | N/A | Saddle
|  $Y=\pm 1$, $x=1$, $R=1$ | $S_{\pm}$ | N/A | Repellor/Attractor
$R=0$ |  $Y=\pm\sqrt{\frac{\mathcal{R}}{3+\mathcal{R}}}$, $x=\mathcal{R}$, $Z=0$ | $dS_{1\pm}$ | $\begin{pmatrix}\mp\sqrt{3R}\\\ \mp\frac{2\sqrt{\mathcal{R}}}{\sqrt{3}}\\\ \mp\sqrt{3\mathcal{R}}(1-\mathcal{R})\end{pmatrix}$ | Attractor/Repellor
|  $Y=0$, $Z=\frac{-x+3x^{2}+3\mathcal{R}-3x\mathcal{R}}{1-x+3x^{2}+3\mathcal{R}-3x\mathcal{R}}$ | $E$ | $\begin{pmatrix}0\\\ \sqrt{\frac{x(7\mathcal{R}+3\mathcal{R}^{2})-x^{2}(9\mathcal{R}+4)+6x^{3}+\mathcal{R}(2-3\mathcal{R})}{2}}\\\ -\sqrt{\frac{x(7\mathcal{R}+3\mathcal{R}^{2})-x^{2}(9\mathcal{R}+4)+6x^{3}+\mathcal{R}(2-3\mathcal{R})}{2}}\end{pmatrix}$ | Saddle
|  $Y=\pm 1$, $x=\mathcal{R}$, $Z=0$ | $dS_{2\pm}$ | N/A | Saddle
|  $Y=\pm 1$, $x=1$, $Z=1$ | $S_{\pm}$ | N/A | Repellor/Attractor
$k=0$ |  $Y=\pm\sqrt{\frac{\mathcal{R}}{3+\mathcal{R}}}$, $x=\mathcal{R}$, $R=0$ ($Z=0$) | $dS_{\pm}$ | $\begin{pmatrix}\mp\frac{4\sqrt{\mathcal{R}}}{\sqrt{3}}\\\ \mp\sqrt{3\mathcal{R}}\\\ \mp\sqrt{3\mathcal{R}}(1-\mathcal{R})\end{pmatrix}$ | Attractor/Repellor
|  $Y=0$, $R=\frac{3x^{2}-3\mathcal{R}x+3\mathcal{R}}{1+3x^{2}-3\mathcal{R}x+3\mathcal{R}}$ | $E$ | $\begin{pmatrix}0\\\ \sqrt{\frac{x(5\mathcal{R}+3\mathcal{R}^{2})-x^{2}(2+9\mathcal{R})+6x^{3}+\mathcal{R}(4-3\mathcal{R})}{2}}\\\ -\sqrt{\frac{x(5\mathcal{R}+3\mathcal{R}^{2})-x^{2}(2+9\mathcal{R})+6x^{3}+\mathcal{R}(4-3\mathcal{R})}{2}}\\\ \end{pmatrix}$ | Saddle
|  $Y=\pm 1$, $x=\mathcal{R}$, $R=0$ ($Z=1$) | $S_{1\pm}$ | N/A | Saddle
|  $Y=\pm 1$, $x=1$, $R=1$ ($Z=0$) | $S_{2\pm}$ | N/A | Repellor/Attractor

Evaluating this Jacobian matrix at each of the fixed points and finding the
eigenvalues then tells us the stability of the fixed points. If the
eigenvalues have non-zero real parts, the fixed point is said to be
hyperbolic. If all real parts of the eigenvalues are positive, then the fixed
point is a repellor, and if the eigenvalues have negative real parts, then the
fixed point is an attractor. If there are positive and negative real parts of
the eigenvalues, then the fixed point is a saddle point or a cusp. Finally, if
the eigenvalues are purely imaginary, and the real parts are therefore zero,
then the fixed point is a centre if the system is linear. For a non-linear
system this requires verification from numerical integration and plots of the
phase space. The case of complex eigenvalues is not relevant here.

We also classify the fixed points by the type of universe model they
represent. The de Sitter models correspond to a cosmological constant and
constant energy density of matter and radiation, which in our system
corresponds to $x^{\prime}=Z^{\prime}=R^{\prime}=0$. $Y$ can vary for
positively and negatively curved de Sitter models. For flat models, $Y$ is
constant, $Y^{\prime}=0$, giving rise to de Sitter fixed points. These occur
at the effective cosmological constants $x=\mathcal{R}$ and $x=1$, and when
$Z=R=0$. An Einstein universe is static, and in our system is represented by
fixed points at $Y^{\prime}=Y=0$. Finally, our system admits singularities,
which occur when the energy densities of matter and radiation become infinite
at $Z=R=1$.

In the following subsections we consider the sub-manifolds of the system,
which simplifies the dynamics by reducing the dimensions of the system Boehmer
and Chan (2014). This helps us to understand the conditions for the existence
of each fixed point, and understand their nature. A summary of the fixed
points of each sub-manifold, along with their eigenvalues and stability
character, can be found in Table LABEL:tab:Sub-manifolds. The fixed points
labelled $E$ represent a static Einstein universe and the fixed points
labelled $dS_{\pm}$ represent spatially flat expanding (+) and contracting (-)
de-Sitter models. $S_{\pm}$ denotes singularities with infinite expansion (+)
or contraction (-) and infinite energy density. Note we cannot compute the
eigenvalues of any fixed point at $Y=\pm 1$, as the Jacobian becomes singular
and we therefore cannot Taylor expand around them. However, from the
Raychaudhuri equation (30), we find that $Y^{\prime}$ is negative along the de
Sitter lines $x=\mathcal{R}$ and $x=1$ around the fixed points at $Y=\pm 1$.

###  III.1 $\Lambda$CDM Dynamics

We begin with the sub-manifolds with an effective cosmological constant and so
are analogous to $\Lambda$CDM, except here $x=1$ and $x=\mathcal{R}$ are
asymptotic values rather than true constants. In each of these cases,
$x^{\prime}=0$ (29).

####  III.1.1 $x=1$

We begin with the $x=1$ sub-manifold. The resulting dynamics is 3-dimensional
(3-D), given by the equations for dark matter (27), radiation (28), and the
Hubble variable (30), which reduces to

$Y^{\prime}=-Y^{2}\sqrt{1-Y^{2}}-\frac{(1-Y^{2})^{\frac{3}{2}}}{6}\\\
\times\left(\frac{Z}{1-Z}+\frac{2R}{1-R}-2\right)\,.$ (35)

In this case, the Friedmann equation (31) becomes

$\frac{Y^{2}}{1-Y^{2}}=\frac{1}{3}+\frac{Z}{3(1-Z)}+\frac{R}{3(1-R)}-\frac{k}{\rho_{*}a^{2}}\,,$
(36)

and the potential (33) is

$U=-\frac{a^{2}}{6}\left(1+\frac{Z}{1-Z}+\frac{R}{1-R}\right)\,.$ (37)

For the sub-manifolds where $x$ is a constant, we express $a$ in terms of $Z$.
Integrating $a^{\prime}$ (18) with respect to $Z$ (27), we find

$a=\left(\frac{1-Z}{Z}\frac{Z_{0}}{1-Z_{0}}\right)^{\frac{1}{3}}\,.$ (38)

We can then solve for the fixed points of the system. There are de Sitter
fixed points where $Z=R=0$, and $Y=\pm 1/2$ and $Y=\pm 1$. The fixed points at
$Y=\pm 1$ and at $Z=R=0$ are coordinate singularities of the de Sitter
spacetime when represented as an FLRW, see Hawking and Ellis (1973). At
$Z=R=1$, (27) and (28) have fixed points, and at $Y=\pm 1$ (35) has fixed
points. Together, the asymptotic points at $Z=R=1$ and $Y=\pm 1$ represent
singularities, with infinite expansion (+) or contraction (-) and infinite
energy density. At the Einstein fixed point, $Y=0$ and the condition
$Y^{\prime}=0$ reduces Eq. (35) to a constraint between $Z$ and $R$:

$\frac{Z_{E}}{1-Z_{E}}+\frac{2R_{E}}{1-R_{E}}-2=0\,,$ (39)

where the subscript $E$ refers to the values of variables at the Einstein
point. On the other hand, matter and radiation do not evolve independently,
and from (2) and (3), they can be related by

$\rho_{r}=\frac{1}{(3H_{0}^{2})^{1/3}}\frac{\Omega_{r}}{\Omega_{m}^{4/3}}\rho_{m}^{4/3},$
(40)

where $H_{0}$, $\Omega_{r}$ and $\Omega_{m}$ can be fixed using Planck values
Aghanim et al. (2018b). This results in (27) and (28) admitting a first
integral $c_{rz}$, which can be used to write a relation between $R$ and $Z$

$\frac{R}{1-R}=c_{rz}\left(\frac{Z}{1-Z}\right)^{\frac{4}{3}},$ (41)

where $c_{rz}$ depends on $\Omega_{m}$, $\Omega_{r}$ and $\rho_{*}$. Together,
(39) and (41) give the values of $Z$ and $R$ at the Einstein fixed point for a
given value of $c_{rz}$. Once the first integral $c_{rz}$ is fixed, a surface
in phase space is defined and the dynamics is effectively 2-dimensional (2-D).
In order to fix $c_{rz}$, we first have to give a value to
$\rho_{*}=\rho_{\Lambda}/\mathcal{R}$, where the future asymptotic value
$\rho_{\Lambda}$ is by definition some fraction of our current observed dark
energy density,

$\rho_{\Lambda}=\alpha\rho_{x,0}\,,$ (42)

and $0<\alpha<1$. For the purpose of our qualitative analysis, we fix
$\alpha=0.5$ as currently we are not far from this asymptotic value as our
Universe is already accelerating. For a reasonable model to solve the
cosmological constant problem, we would fix $10^{-120}<\mathcal{R}<10^{-60}$
so the dark energy evolved between the high estimate of the contributions to
the effective cosmological constant from QFT and the low observed value.
However, the dynamics cannot be depicted in plots of the phase space for these
small values of $\mathcal{R}$. We therefore choose $\mathcal{R}=0.05$ in order
to illustrate the evolution of our variables, given that the value of
$\mathcal{R}$ will not qualitatively change the dynamics. We therefore find

$c_{rz}=\frac{\Omega_{r}}{(\Omega_{m})^{\frac{4}{3}}}\left(\frac{\alpha\Omega_{\Lambda}}{\mathcal{R}}\right)^{\frac{1}{3}}=\Omega_{r}\left(\frac{\Omega_{\Lambda}}{\Omega_{m}^{4}}\right)^{\frac{1}{3}}\left(\frac{\alpha}{\mathcal{R}}\right)^{\frac{1}{3}}\\\
\simeq 0.00083\,,$ (43)

where we have used $\Omega_{\Lambda}=0.6889$ and $\Omega_{m}=0.3111$, and
calculated $\Omega_{r}=9.1824\times 10^{-5}$ using $\Omega_{m}$ and the
redshift of matter-radiation equality $z=3387$ Aghanim et al. (2018b).

The phase space for the $x=1$ sub-manifold is shown in Fig. 1, where the
dynamics takes place on the first integral surface $c_{rz}$, shown in yellow.
The green outer most thick curve is the flat Friedmann separatrix (FFS), which
separates the closed models (in between the green curves) from the open
models. This separatrix on the 2-D surface is a sub-set of the general
separatrix hypersurface, which has the general form

$\frac{Y^{2}}{1-Y^{2}}=\frac{x}{3}+\frac{Z}{3(1-Z)}+\frac{R}{3(1-R)},$ (44)

where here $x=1$ and $R$ and $Z$ are related by (41). The inner most curve is
the closed Friedmann separatrix (CFS), which in general is the hypersurface
consisting of the Einstein fixed points present in the phase space, and
trajectories asymptotic to them. This shows the boundaries between the
different types of closed model. The CFS has the general form

$\frac{Y^{2}}{1-Y^{2}}=\frac{x}{3}+\frac{Z}{3(1-Z)}+\frac{R}{3(1-R)}\\\
-\frac{1}{a^{2}}\left[\frac{x_{E}}{3}+\frac{Z_{E}}{3(1-Z_{E})}+\frac{R_{E}}{3(1-R_{E})}\right]\,,$
(45)

where again here $x=x_{E}=1$ and the relations in Eq. (38) and Eq. (41) apply.

Figure 1: The phase space for the $x=1$ sub-manifold, with the dynamics shown
on the yellow $c_{rz}$ surface, given by (41). The green curves represent flat
models asymptotic to the de Sitter fixed points at $Z=0$ and singularities at
$Z=1$. The fixed point at $Y=0$ (a saddle) represents the Einstein model, and
the black separatrices asymptotic to it separate the different types of closed
models in the phase space.

The 2-D dynamics taking place on the first integral surface in the 3-D phase
space can also be projected onto the $Z$-$Y$ plane in order to make the
dynamics more visible, again using (41). The projected sub-manifold can be
seen in Fig. 2. The vertical black lines along $Z=0$ and $Z=1$ show the de
Sitter models in the phase space, where $x^{\prime}=Z^{\prime}=R^{\prime}=0$.
Expanding (contracting) open and flat models evolve from (to) a singularity to
(from) a flat de Sitter fixed point. Closed models outside the CFS evolve in
the same way. Within the CFS there are bouncing models, which contract until
they reach a maximum energy density at a minimum $a$ where they bounce, and
then expand again. The bounce occurs when the dark energy component is
dominant over the matter and radiation, and the magnitude of the curvature
term is equal to the sum of all other contributions to the Friedmann equation
(36), giving $Y=0$. There are also turn-around models, which expand until they
reach a minimum energy density at a maximum $a$ before contracting again, that
evolve between the two singularities. The turn-around occurs at $Y=0$, when
the dark matter and radiation are dominant over the dark energy, and the
magnitude of the curvature term becomes large enough to equal the sum of all
other contributions to the Friedmann equation (36). Note that for all
trajectories in the $Z$-$Y$ plane, the expansion (contraction) is accelerating
($a^{\prime\prime}>0$) to the left of the Einstein point. Therefore in a
$\Lambda$CDM cosmology, such as the one we are illustrating here, bouncing
models are always accelerating, and re-collapsing models are never
accelerating.

Figure 2: $x=1$ sub-manifold, top panel: the projection on the $Z$-$Y$ plane
of the 2-D dynamics taking place on the first integral surface in the 3-D
phase space in Fig. 1; the qualitative behaviour is typical of any
$\Lambda$CDM model, independently of the parameter values. Bottom panel: the
corresponding potential in Eq. (37), where trajectories of the same colour in
the two panels correspond to each other. Bounded trajectories are below the
maximum of the potential, and unbounded are above.

There are also two other 2-D submanifolds at $x=1$; one for $Z=0$ (radiation
only) and the other when $R=0$ (matter only). The dynamics of these two sub-
cases is qualitatively the same as in Fig. 2. In the $Z=0$ case the Einstein
point occurs at $R_{E}=1/2$, and in the $R=0$ case it occurs at $Z_{E}=2/3$.

####  III.1.2 $x=\mathcal{R}$

Next we consider the case where $x=\mathcal{R}$. The remaining 3-dimensional
dynamics is given by the equations for the dark matter (27), radiation (28),
and the Hubble function (30), which reduces to

$Y^{\prime}=-Y^{2}\sqrt{1-Y^{2}}-\frac{(1-Y^{2})^{\frac{3}{2}}}{6}\\\
\times\left(\frac{Z}{1-Z}+\frac{2R}{1-R}-2\mathcal{R}\right)\,.$ (46)

The Friedmann equation (31) when $x=\mathcal{R}$ becomes

$\frac{Y^{2}}{1-Y^{2}}=\frac{\mathcal{R}}{3}+\frac{Z}{3(1-Z)}+\frac{R}{3(1-R)}-\frac{k}{\rho_{*}a^{2}}\,,$
(47)

and the potential (33) is

$U=-\frac{a^{2}}{6}\left(\mathcal{R}+\frac{Z}{1-Z}+\frac{R}{1-R}\right)\,.$
(48)

The dynamics of this system is qualitatively the same as in Fig. 1. The
singularities are present at $Z=R=1$ and $Y=\pm 1$, and there are de Sitter
points at $Y=\sqrt{\mathcal{R}/(3+\mathcal{R})}$ and $Y=\pm 1$ along $Z=R=0$.
At the Einstein fixed point, the Hubble expansion variable (46) reduces to

$\frac{Z_{E}}{1-Z_{E}}+\frac{2R_{E}}{1-R_{E}}-2\mathcal{R}=0\,,$ (49)

which is solved numerically using (41) to find the values of $Z_{E}$ and
$R_{E}$ at this point.

As before, this sub-manifold can be reduced to 2-D at $Z=0$ and at $R=0$,
which are both qualitatively the same as in Fig. 2. The de Sitter points and
singularities are as outlined above. The Einstein point in the $Z=0$ case
occurs at $R_{E}=\mathcal{R}/(1+\mathcal{R})$, and in the $R=0$ case at
$Z_{E}=2\mathcal{R}/(1+2\mathcal{R})$.

###  III.2 Dynamic Dark Energy

We now consider the sub-manifolds with dynamic dark energy so that $x$ is not
fixed, and therefore $x^{\prime}\neq 0$ in Eq. (29).

####  III.2.1 $Z=0$ and $R=0$

We consider the sub-manifold where matter and radiation are not present, and
there is only dark energy. In this case $Z^{\prime}=0$ in Eq. (27) and
$R^{\prime}=0$ in Eq. (28). The remaining dynamics is 2-D and given by the
equations for the dark energy (29) and the Hubble expansion variable (30)
which becomes

$Y^{\prime}=-Y^{2}\sqrt{1-Y^{2}}-\frac{(1-Y^{2})^{\frac{3}{2}}}{6}\\\
\times[x(1+3\mathcal{R})-3\mathcal{R}-3x^{2}].$ (50)

When $Z=R=0$, the Friedmann equation (31) becomes

$\frac{Y^{2}}{1-Y^{2}}=\frac{x}{3}-\frac{k}{\rho_{*}a^{2}}\,,$ (51)

and the potential (33) is

$U=-\frac{xa^{2}}{6}\,.$ (52)

In order to plot the potential and the CFS in the phase space, we must express
the scale factor $a$ in terms of $x$. Integrating $a^{\prime}$ (18) with
respect to $x$ (13), we find

$a(x)=\left(c_{a}\frac{x-\mathcal{R}}{1-x}\right)^{\frac{-1}{3(1-\mathcal{R})}}\,,$
(53)

i.e. the inverse of Eq. (22), where $c_{a}$ is as in Eq. (23). There are de
Sitter points at $Y=\sqrt{\mathcal{R}/(3+\mathcal{R})}$ and $Y=\pm 1$ along
$x=\mathcal{R}$, and at $Y=\pm 1/2$ and $Y=\pm 1$ along $x=1$. The fixed
points at $Y=\pm 1$ and at $x=\mathcal{R}$ and $x=1$ are coordinate
singularities of the de Sitter spacetime when represented as an FLRW, see
Hawking and Ellis (1973). At the Einstein fixed point, $Y=0$, Eq. (50) reduces
to

$x_{E}(1+3\mathcal{R})-3\mathcal{R}-3x_{E}^{2}=0.$ (54)

Rearranging for $x_{E}$, we find

$x_{E\pm}=\frac{1+3\mathcal{R}}{6}\pm\sqrt{\frac{(1+3\mathcal{R})^{2}}{36}-\mathcal{R}}.$
(55)

Given the condition $0<\mathcal{R}<1$, both roots in Eq. (55) will be
positive, provided that the term under the square root is greater or equal to
zero. Thus, for Einstein points to exist, we find

$\mathcal{R}\leq\mathcal{R}_{M}\equiv\frac{5}{3}-\frac{2\sqrt{6}}{3}\,.$ (56)

Therefore, there are three cases for this sub-manifold: no Einstein points are
admitted when $\mathcal{R}>\mathcal{R}_{M}$, two Einstein points exist when
$\mathcal{R}<\mathcal{R}_{M}$, and there is a limiting case in between when
$\mathcal{R}=\mathcal{R}_{M}$ where one Einstein point is admitted. The open
and flat models evolve in the same way for each case. Expanding (contracting)
open models evolve between two de Sitter fixed points, from (to) an open
geometry to (from) a flat geometry. Expanding (contracting) flat models evolve
between two de Sitter fixed points along the FFS. In general, all closed
models admit a bounce, however the qualitative behaviour changes depending on
the number of Einstein points admitted.

The $\mathcal{R}>\mathcal{R}_{M}$ case is shown in Fig. 3, where the system
admits no Einstein fixed points. The two vertical black lines along
$x=\mathcal{R}$ and $x=1$ highlight the de Sitter models in the phase space.
All closed models bounce, evolving between two de Sitter fixed points. The
limiting $\mathcal{R}=\mathcal{R}_{M}$ case is shown in Fig. 4. Here, the
system admits one Einstein point, which is a cusp, so all closed models
bounce, evolving between two de Sitter fixed points.

Finally, the case where $\mathcal{R}<\mathcal{R}_{M}$ is shown in Fig. 5. Two
Einstein points are admitted: the fixed point at $x_{E}\simeq 0.07$ is a
saddle point, which is part of the CFS (the black loop in Fig. 5), and the
other at $x_{E}\simeq 0.28$ is a centre. Closed models within the CFS either
bounce once, evolving between two de Sitter fixed points, or are cyclic around
the centre fixed point, repeatedly contracting until they bounce, then
expanding until they turn-around. Closed models outside the CFS also bounce
and evolve between two de Sitter fixed points. These closed models outside the
CFS evolve with an early- and late-time acceleration, connected by a
decelerating period for $x$ between the two Einstein points. This is the only
case for this set of sub-manifolds to have a period of deceleration, as here
$\mathcal{R}$ is small enough for the effective EoS to be larger than $-1/3$
for $x$ values between the two Einstein points, which violates the condition
for acceleration in Eq. (17).

Figure 3: Top panel: the phase space of the $Z=0$, $R=0$ sub-manifold, with
$\mathcal{R}>\mathcal{R}_{M}$. Bottom panel: corresponding potential in Eq.
(52), where trajectories of the same colour in the two panels correspond to
each other. In this case, there are no Einstein points present, and so all
closed models bounce, however there is never a decelerating phase. Figure 4:
Top panel: the phase space of the $Z=0$, $R=0$ sub-manifold, with
$\mathcal{R}=\mathcal{R}_{M}$. Bottom panel: corresponding potential in Eq.
(52), where trajectories of the same colour in the two panels correspond to
each other. One Einstein point (cusp) is present, and all closed models
bounce, however there is never a decelerating period. Figure 5: Top panel: the
phase space of the $Z=0$, $R=0$ sub-manifold, with
$\mathcal{R}<\mathcal{R}_{M}$. Bottom panel: corresponding potential in Eq.
(52), where trajectories of the same colour in the two panels correspond to
each other. Two Einstein points at $x_{E}\simeq 0.07$ (a saddle) and
$x_{E}\simeq 0.28$ (a centre) are present, therefore all closed models bounce,
either once or in repeated cycles. Deceleration occurs for $x$ between the two
Einstein points. The closed models outside the CFS evolve with an initial and
final acceleration, with a decelerating period in between.

####  III.2.2 $Z=0$, $R=0$

A three-dimensional sub-manifold exists where $Z=0$, and another when $R=0$.
Adding only dark matter to the system has the same effect as adding only
radiation; the difference between the two are the values of the first integral
of $Z(x)$ and $R(x)$ needed to obtain each case. The dynamics is qualitatively
the same as the cases we present of the full system in section IV, where we
include both dark matter and radiation in our analysis, therefore we do not
include these sub-manifolds explicitly here.

####  III.2.3 Spatially flat models: $k=0$

Finally, we consider the sub-manifold where the models are all spatially flat,
shown in Fig. 6, where we have compactified the 4-dimensional (4-D) phase
space to 2-D on the $x$-$Y$ plane. In order to do this, we express $z$ and $r$
as functions of $x$. From Eq. (11) and Eq. (13), we find

$z(x)=c_{z}\left(\frac{x-\mathcal{R}}{1-x}\right)^{\frac{1}{1-\mathcal{R}}}\,,$
(57)

where the constant of integration is

$c_{z}=z_{0}\left(\frac{1-x_{0}}{x_{0}-\mathcal{R}}\right)^{\frac{1}{1-\mathcal{R}}}\,.$
(58)

Here, to express $c_{z}$, we set $z_{0}$ using the Friedmann equation (31)
when $k=0$:

$\frac{Y^{2}}{1-Y^{2}}=\frac{x}{3}+\frac{Z}{3(1-Z)}+\frac{R}{3(1-R)}\,.$ (59)

We can then express $z_{0}$ as

$z_{0}=\frac{Z_{0}}{(1-Z_{0})}=\frac{3Y_{0}^{2}}{1-Y_{0}^{2}}-x_{0}-\frac{R_{0}}{(1-R_{0})}\,.$
(60)

Then to express $r$ in terms of $x$, from Eq. (12) and Eq. (13) we find

$r(x)=c_{r}\left(\frac{x-\mathcal{R}}{1-x}\right)^{\frac{4}{3(1-\mathcal{R})}}\,,$
(61)

where the constant of integration is

$c_{r}=r_{0}\left(\frac{1-x_{0}}{x_{0}-\mathcal{R}}\right)^{\frac{4}{3(1-\mathcal{R})}}\,.$
(62)

We set $c_{r}$ in a similar way to $c_{rz}$, using the density parameters
Aghanim et al. (2018b), and our definitions of $\mathcal{R}$ (10) and
$\rho_{\Lambda}$ (42). Keeping our choice of $\mathcal{R}=0.05$, we then find

$c_{r}=\frac{\mathcal{R}}{\alpha}\frac{\Omega_{r}}{\Omega_{\Lambda}}\left(\frac{\alpha-\mathcal{R}}{\mathcal{R}(1-\alpha)}\right)^{\frac{4}{3(1-\mathcal{R})}}\simeq
0.00077\,.$ (63)

The potential is as in Eq. (33). The shape of the potential in this case
changes with the initial conditions, and as $k=0$ the total energy of the
system $E=0$ (21). Therefore the behaviour of each trajectory in the phase
space can be seen along $U=0$ in the bottom panel of Fig. 6. Here, the green
separatrix represents models where $Z=0$, so there is no dark matter
component, and the innermost black curve is the CFS. There are no de Sitter
lines along $x=\mathcal{R}$ or $x=1$ as in the flat case, $Z$ varies for
constant $x$. We do not consider models within the FFS as $Z<0$ in this
region. The models outside this separatrix always have positive dark matter
energy density, $Z>0$. Expanding (contracting) models evolve from (to) an
initial (a future) singularity to (from) a de Sitter fixed point. Since these
models do not exhibit any bouncing behaviour, they are not interesting with
respect to our investigation.

Figure 6: Top panel: the phase space of the $k=0$ sub-manifold. Bottom panel:
corresponding potential in Eq. (33), where trajectories of the same colour in
the two panels correspond to each other. Models within the green $Z=0$
separatrix always have negative dark matter energy density $Z<0$. The shape of
the potential changes with the initial conditions, and as these are all flat
models $k=0$ the total energy of the system $E=0$ (21). Therefore, the
behaviour of each trajectory can be seen along $U=0$ for each potential.

##  IV The Full System

To show the full range of dynamics, we project the full 4-D dynamics to 2-D on
the $x$-$Y$ plane, by expressing $z$ and $r$ in terms of $x$ as in Eq. (57)
and Eq. (61), respectively. We set $c_{r}$ as in Eq. (63), keeping
$\mathcal{R}=0.05$ in order for the dynamics to be visible, and to be able to
analyse the full range of the dynamics. We do not fix $c_{z}$ in the full
system, but instead investigate how the dynamics varies for different values
of this parameter. In order to produce plots of the phase space, we must also
express the scale factor $a$ in terms of $x$, using Eq. (53). The Friedmann
equation (31) can therefore be written just in terms of $x$ and $Y$, and the
potential (20) can be written in terms of $x$ only. The fixed points for the
full system are given in Table 2.

As before, $E$ represents a static Einstein universe, $dS_{\pm}$ represents
spatially flat expanding (+) and contracting (-) de-Sitter models, and
$S_{\pm}$ represents singularities with infinite expansion (+) or contraction
(-) and infinite energy density. Note that the $dS_{2\pm}$ fixed points are
coordinate singularities of the de Sitter spacetime when represented as an
FLRW, see Hawking and Ellis (1973).

Table 3 shows the eigenvalues of the fixed points for this system, and the
linear stability classification for each point is given in Table 4. Again, we
cannot include the eigenvalues of the fixed points at $Y=\pm 1$ in Table 3, as
the Jacobian becomes singular. Therefore, we find that the singularity at
$Y=+1$ is a repellor and the singularity at $Y=-1$ is an attractor, and the de
Sitter points at $Y=\pm 1$ are saddle points, with trajectories moving away
from the $Y=+1$ point and towards the $Y=-1$ point along $x=\mathcal{R}$. Note
that the linear stability analysis of the Einstein points needs to be done for
each individual case, as these can be a centre, a cusp or a saddle point,
whereas the linear stability of the de-Sitter points and singularities are
general for all cases considered here.

Name | $x$ | $Y$ | $Z$ | $R$
---|---|---|---|---
$E$ | $x$ | $0$ | $Z$ | $\frac{-3x^{2}(1-Z)+Z(1+3\mathcal{R})-3\mathcal{R}+x(1-Z)(1+3\mathcal{R})}{-2-3x^{2}(1-Z)-3\mathcal{R}+3Z(1+\mathcal{R})+x(1-Z)(1+3\mathcal{R})}$
dS1± | $\mathcal{R}$ | $\pm\sqrt{\frac{\mathcal{R}}{3+\mathcal{R}}}$ | $0$ | $0$
dS2± | $\mathcal{R}$ | $\pm 1$ | $0$ | $0$
S± | $1$ | $\pm 1$ | $1$ | $1$
Table 2:  The fixed points of the full system. $E$ denotes an Einstein universe, $dS_{\pm}$ an expanding (+) or contracting (-) de-Sitter universe, and $S_{\pm}$ a singularity with infinite expansion (+) or contraction (-). Name | $\lambda_{1}$ | $\lambda_{2}$ | $\lambda_{3}$ | $\lambda_{4}$
---|---|---|---|---
$E$ | 0 | 0 | $\frac{\sqrt{18x^{3}(Z-1)-9x^{2}(Z-1)(3\mathcal{R}+1)+x(Z-1)\left(9\mathcal{R}^{2}+18\mathcal{R}-1\right)+Z\left(-9\mathcal{R}^{2}+9\mathcal{R}+1\right)+9(\mathcal{R}-1)\mathcal{R}}}{\sqrt{6}\sqrt{Z-1}}$ | $-\lambda_{3}$
$dS_{1\pm}$ | $\mp 4\sqrt{\frac{\mathcal{R}}{3}}$ | $\mp\sqrt{3\mathcal{R}}$ | $\mp 2\sqrt{\frac{\mathcal{R}}{3}}$ | $\mp\sqrt{3\mathcal{R}}(1-\mathcal{R})$
Table 3:  The eigenvalues for the fixed points of the full system given in Table 2. Name | Stability Character
---|---
$E$ | Centre, Saddle or Cusp
$dS_{1+}$ | Attractor
$dS_{1-}$ | Repellor
$dS_{2\pm}$ | Saddle
$S_{+}$ | Repellor
$S_{-}$ | Attractor
Table 4:  The linear stability character for the fixed points of the full
system given in Table 2.

At an Einstein point, the remaining expression left in the Hubble expansion
variable (30), which we call $f(x)$, is

$f(x)=z(x)+2r(x)+x(1+3\mathcal{R})-3\mathcal{R}-3x^{2}\,,$ (64)

which we need to solve numerically. An Einstein point exists whenever this
expression is equal to zero. Taking the limit of $f(x)$ when
$x\rightarrow\mathcal{R}$ we find $f(x)\rightarrow-2\mathcal{R}$, and taking
the limit of $f(x)$ when $x\rightarrow 1$ we find $f(x)\rightarrow\infty$.
Therefore, between $x=\mathcal{R}$ and $x=1$, $f(x)$ is always equal to zero
at least once, and so at least one Einstein point always exists. The possible
number of Einstein points admitted by the system can be seen in Fig. 7. For
each curve in the plot, which have a specific value of $c_{z}$, an Einstein
point occurs when the curve meets the $x$-axis at $f(x)=0$. Between the
smallest and largest values of $c_{z}$, we find the maximum number of Einstein
points the system may admit is three.

Figure 7: A plot of $f(x)$ (64) with varying values of $c_{z}$ to show the
possible number of Einstein points admitted by the system. Starting from the
top we have $c_{z}>0.38$ (orange), $c_{z}\simeq 0.38$ (green),
$0.20<c_{z}<0.38$ (magenta), $c_{z}\simeq 0.20$ (blue) and $c_{z}<0.20$
(purple). The system always admits at least one Einstein point, and a maximum
of three can exist.

The cases where two Einstein points are admitted are limiting cases, where two
roots coincide, and where the maximum or minimum of the $f(x)$ curve touches
but does not cross $f(x)=0$. To find the values of $c_{z}$ in these limiting
cases, we first need to find the derivative of $f(x)$ (64) with respect to
$x$, which is

$\frac{df(x)}{dx}=\frac{1}{3(1-x)(x-\mathcal{R})}[18x^{3}-3x^{2}(7+9\mathcal{R})\\\
+3x(1+10\mathcal{R}+3\mathcal{R}^{2})-3\mathcal{R}(1+3\mathcal{R})+3z+8r]\,.$
(65)

Setting this equal to zero and rearranging for $c_{z}$, we find

$c_{z}=\frac{1}{3}\left(\frac{1-x}{x-\mathcal{R}}\right)^{\frac{1}{1-\mathcal{R}}}[-18x^{3}+3x^{2}(7+9\mathcal{R})\\\
-3x(1+10\mathcal{R}+3\mathcal{R}^{2})+3\mathcal{R}(1+3\mathcal{R})-8r]\,,$
(66)

which we substitute into $f(x)=0$ (64). Solving this equation for $x$, we find
$x\simeq 0.25$ and $x\simeq 0.60$, which from (66) correspond to $c_{z}\simeq
0.20$ and $c_{z}\simeq 0.38$. These values are shown in Fig. 7 by the blue and
green curves, respectively.

In the general case where there are three Einstein points, there is a limiting
case. In general, there are two separate CFS curves when three Einstein points
are present; one through the lower and one through the upper Einstein point.
In the limiting case, the two CFS curves coincide and the maxima of the
potential are equal. To find where this occurs, we choose an initial value for
$c_{z}$, then equate the potential U (20) at the upper and lower Einstein
points. We iterate until we converge on a value of $c_{z}$, which occurs at
$c_{z}\simeq 0.32$.

Therefore, we find seven different dynamical cases for this system. The open
and flat models in each case all evolve in the same way. Expanding
(contracting) models evolve from (to) an initial (a future) singularity to
(from) a flat de-Sitter spacetime. Closed models within the FFS but outside
the CFS also evolve in this way in each case, as the curvature is never large
enough for a bounce or turn-around to occur. The behaviour of models within
the CFS changes depending on the value of $c_{z}$ as this affects the number
and character of the Einstein points present. We present these sub-cases in
the following subsections, starting with the smallest range of $c_{z}$ and
increasing the value with each case.

###  IV.1 1 Einstein Point

The phase space for the range $c_{z}<0.20$ is given in Fig 8. One Einstein
point exists at $x_{E}\simeq 0.97$ which is a saddle point, and corresponds to
the maximum of the potential. Most of the phase space within the CFS to the
left of the Einstein point consists of bouncing models driven by the dark
energy, while matter and radiation are subdominant, which evolve between two
de Sitter fixed points. For the models which are dominated by matter and
radiation, to the right of the Einstein point, the evolution is between two
singularities with a turn-around. For all trajectories, the expansion
(contraction) is always accelerating ($a^{\prime\prime}>0$) to the left of the
Einstein point. Therefore, for this case, the bouncing models are always
accelerating, and the turn-around models are never accelerating.

Figure 8: $c_{z}<0.20$ case, top panel: the projection of the full 4-D
dynamics on the 2-D $x$-$Y$ plane. Bottom panel: the corresponding potential
as in Eq. (20), where $a$ is given by Eq. (53), $z$ by Eq. (57) and $r$ by Eq.
(61), and trajectories of the same colour in the two panels correspond to each
other. One Einstein fixed point $x_{E}\simeq 0.97$ (a saddle) exists, which
corresponds to the maximum of the potential, therefore we obtain bouncing
trajectories for $x<0.97$ within the CFS, however they are always
accelerating.

###  IV.2 2 Einstein Points

The phase space for the first limiting case with two Einstein points where
$c_{z}\simeq 0.20$ is shown in Fig. 9. The Einstein point at $x_{E}\simeq
0.87$ is a saddle point, which corresponds to the maximum of the potential and
is part of the outer CFS. Within this CFS, dark matter and radiation dominate
for $x>0.87$ where turn-around models evolve between two singularities. For
$x<0.87$, dark energy is dominant and we obtain bouncing models which evolve
between two de Sitter points. The Einstein point at $x_{E}\simeq 0.25$ is a
cusp (corresponding to the two de Sitter points coinciding), which corresponds
to the horizontal point of inflection in the potential, and is part of the
innermost CFS. Within this CFS, the bouncing models also evolve between the
two de Sitter points. To the left of the Einstein point at $x_{E}\simeq 0.87$,
trajectories are always accelerating, therefore the bouncing models are always
accelerating, and the re-collapsing models are always decelerating.

Figure 9: The first limiting case where $c_{z}\simeq 0.20$, top panel: the
projection of the full 4-D dynamics on the 2-D $x$-$Y$ plane. Bottom panel:
the corresponding potential as in Eq. (20), where $a$ is given by Eq. (53),
$z$ by Eq. (57) and $r$ by Eq. (61), and trajectories of the same colour in
the two panels correspond to each other. Two Einstein fixed points occur. The
fixed point $x_{E}\simeq 0.25$ (a cusp) corresponds to the horizontal point of
inflection in the potential, and the fixed point at $x_{E}\simeq 0.87$ (a
saddle) corresponds to the maximum of the potential. We obtain bouncing
trajectories for $x<0.87$ within the CFS, however these are always
accelerating.

###  IV.3 3 Einstein Points

####  IV.3.1 Extra bouncing region

The first case where three Einstein points exist is shown in Fig. 10. The
Einstein point at $x_{E}\simeq 0.19$ is a saddle point, which corresponds to a
local maximum of the potential, and is part of the innermost CFS. Within this
CFS, bouncing models evolve between two de Sitter points. Cyclic models are
present around the Einstein fixed point at $x_{E}\simeq 0.34$, which is a
centre, corresponding to a local minimum of the potential. These models
contract until they bounce, and expand until they reach a turn-around. The
Einstein point at $x_{E}\simeq 0.83$ is another saddle point, corresponding to
the maximum of the potential, and is part of the outermost CFS. Within this
CFS, turn-around models which evolve between two singularities are present for
$x>0.83$, where dark matter and radiation are dominant over the dark energy.
For $x<0.83$, dark energy is the dominant component, and bouncing models are
present between the two CFS curves which evolve between two de Sitter points.
These bounces evolve with an early- and late-time acceleration, connected by a
period of deceleration. These are the physically important models, which we
discuss in section IV.6.

Figure 10: $0.20<c_{z}<0.32$ case, top panel: the projection of the full 4-D
dynamics on the 2-D $x$-$Y$ plane. Bottom panel: the corresponding potential
as in Eq. (20), where $a$ is given by Eq. (53), $z$ by Eq. (57) and $r$ by Eq.
(61), and trajectories of the same colour in the two panels correspond to each
other. The Einstein fixed points $x_{E}\simeq 0.19$ and $x_{E}\simeq 0.83$
(saddle) correspond to the maxima of the potential, and $x_{E}\simeq 0.34$
(centre) to a local minimum. Bouncing behaviour is obtained for $x<0.83$
within the outer CFS, and cyclic behaviour is obtained for $x>0.19$ within the
inner CFS. The cyclic models evolve with early acceleration and late-time
deceleration, and the bouncing models evolve with early- and late-time
acceleration connected by a decelerating period.

####  IV.3.2 Limiting case

A limiting case exists for the general case with three Einstein points, when
the two maxima of the potential have the same value. This case is shown in
Fig. 11. Here, the two saddle Einstein fixed points at $x_{E}\simeq 0.17$ and
$x_{E}\simeq 0.76$ correspond to the two maxima of the potential, and the two
CFS curves corresponding to each of the fixed points merge and coincide.
Within the CFS, turn-around models which evolve between two singularities are
present in the $x>0.76$ region of the phase space, where dark matter and
radiation are dominant. Bouncing models which evolve between two de Sitter
points are present when $x<0.17$. Cyclic models are present around the
Einstein point at $x_{E}\simeq 0.43$ (a centre) which corresponds to the local
minimum value of the potential. These models contract until they bounce, and
expand until they reach a turn-around. In this case, the cyclic models evolve
with an early acceleration and late-time deceleration, the bouncing models
always accelerate, and the re-collapsing models never accelerate.

Figure 11: The limiting case $c_{z}\simeq 0.32$, where the two CFS curves
coincide. Top panel: the projection of the full 4-D dynamics on the 2-D
$x$-$Y$ plane. Bottom panel: the corresponding potential as in Eq. (20), where
$a$ is given by Eq. (53), $z$ by Eq. (57) and $r$ by Eq. (61), and
trajectories of the same colour in the two panels correspond to each other.
The fixed points $x_{E}\simeq 0.17$ and $x_{E}\simeq 0.76$ (saddle) correspond
to the maxima of the potential which have the same value of the potential, and
$x_{E}\simeq 0.43$ (centre) to a local minimum. Cyclic behaviour is obtained
for $0.17<x<0.76$ and bouncing behaviour for $x<0.17$, both within the CFS.

####  IV.3.3 Extra turn-around region

The final case with three Einstein points is given in Fig. 12. The Einstein
point at $x_{E}\simeq 0.72$ is a saddle point, which corresponds to the local
maximum value of the potential, and is part of the innermost CFS. Within this
CFS, turn-around models which evolve between two singularities are present for
the $x>0.72$ region of the phase space. These models never accelerate, and
dark matter and radiation are dominant. For $x<0.72$, dark energy is the
dominant component. Cyclic models are present around the centre Einstein fixed
point at $x_{E}\simeq 0.48$, which has a local minimum value of the potential.
These models contract until they bounce, then expand until they turn-around
and evolve with early acceleration and late-time deceleration. The Einstein
fixed point at $x_{E}\simeq 0.16$ is another saddle point which is part of the
outermost CFS, and occurs at the maximum of the potential. Bouncing models
evolve between two de Sitter points for $x<0.16$, which always accelerate. For
the $x>0.16$ region of the phase space between the two CFS curves, turn-around
models are present which evolve between two singularities. In these models,
radiation and dark matter are initially dominant, but decrease more rapidly
than the dark energy, which becomes dominant as the trajectory passes the
fixed point at $x_{E}\simeq 0.72$. After the turn-around, the radiation and
matter increase more rapidly than the dark energy, and become dominant again
when $x>0.72$, where the models approach a singularity. These re-collapsing
models evolve with an early and late-time deceleration connected by an
accelerating period.

Figure 12: $0.32<c_{z}<0.38$ case, top panel: the projection of the full 4-D
dynamics on the 2-D $x$-$Y$ plane. Bottom panel: the corresponding potential
as in Eq. (20), where $a$ is given by Eq. (53), $z$ by Eq. (57) and $r$ by Eq.
(61), and trajectories of the same colour in the two panels correspond to each
other. The fixed points $x_{E}\simeq 0.16$ and $x_{E}\simeq 0.72$ (saddle)
correspond to the maxima of the potential, and $x_{E}\simeq 0.48$ (centre) to
a local minimum. Cyclic behaviour is obtained for $x<0.72$ within the inner
CFS and bouncing behaviour for $x<0.16$ within the outer CFS. Cyclic models
evolve with early acceleration and late-time deceleration, and the bouncing
models always accelerate.

###  IV.4 2 Einstein Points

The phase space for the system for $c_{z}\simeq 0.38$ is shown in Fig. 13,
which is the second limiting case with two Einstein points. The Einstein point
at $x_{E}\simeq 0.60$ is a cusp (corresponding to the two de Sitter points
coinciding), which corresponds to the horizontal point of inflection in the
potential, and is part of the innermost CFS. Within this CFS, turn-around
models evolve between two singularities, where the dark matter and radiation
are the dominant components, and so never accelerate. The Einstein point at
$x_{E}\simeq 0.16$ is a saddle which is part of the outermost CFS, and
corresponds to the maximum value of the potential. For the $x<0.16$ region of
the phase space, bouncing models evolve between two de Sitter points. In these
models, the dark energy always dominates, and they always accelerate. Between
the two CFS curves, turn-around models are present which evolve between two
singularities, and never accelerate. Dark matter and radiation are dominant
for $x>0.60$, and dark energy is dominant for the $x<0.60$ region of the phase
space. After the turn-around, the radiation and dark matter components
increase more quickly than the dark energy throughout the collapse, and become
dominant on approach to the singularity.

Figure 13: $c_{z}\simeq 0.38$ case, top panel: the projection of the full 4-D
dynamics on the 2-D $x$-$Y$ plane. Bottom panel: the corresponding potential
as in Eq. (20), where $a$ is given by Eq. (53), $z$ by Eq. (57) and $r$ by Eq.
(61), and trajectories of the same colour in the two panels correspond to each
other. The fixed point $x_{E}\simeq 0.16$ (saddle) corresponds to the maximum
of the potential, and $x_{E}\simeq 0.60$ (cusp) to the horizontal point of
inflection in the potential. Bouncing behaviour is obtained within the CFS for
$x<0.16$, however these models always accelerate.

###  IV.5 1 Einstein Point

The final case for this system is given in Fig. 14. Qualitatively, the
behaviour here is the same as in Fig. 8, except that the Einstein point in
this case is much closer to $x=\mathcal{R}$ at $x_{E}\simeq 0.14$. This is a
saddle point that corresponds to the maximum value of the potential for the
system, and is part of the CFS. Within the CFS, bouncing models are present
for the $x<0.14$ region of the phase space, in which dark energy is always
dominant. For the $x>0.14$ region, turn-around models which evolve between two
singularities exist. These models are initially dominated by matter and
radiation, and then by dark energy. After the turn-around, the dark matter and
radiation increase more quickly than the dark energy, and become dominant
again as the singularity is approached. For all trajectories, the expansion
(contraction) is always accelerating to the left of the Einstein point,
therefore the bouncing models always accelerate, and the re-collapse models
always decelerate.

Figure 14: $c_{z}>0.38$ case, top panel: the projection of the full 4-D
dynamics on the 2-D $x$-$Y$ plane. Bottom panel: the corresponding potential
as in Eq. (20), where $a$ is given by Eq. (53), $z$ by Eq. (57) and $r$ by Eq.
(61), and trajectories of the same colour in the two panels correspond to each
other. The fixed point $x_{E}\simeq 0.14$ (saddle) corresponds to the maximum
of the potential. Bouncing models are obtained within the CFS for $x<0.14$,
which always accelerate.

###  IV.6 Acceleration regions

In order to determine which of the cases are viable, we can consider the
acceleration regions in each phase space. For our models to be feasible, they
must have an early and late time acceleration, connected by a decelerating
phase where large-scale structure can form. The acceleration equation for this
system is given by

$\frac{a^{\prime\prime}}{a}=z+2r-3\mathcal{R}+(1+3\mathcal{R})x-3x^{2}\,.$
(67)

We can find the boundaries between accelerating and decelerating regions by
setting the acceleration equation (67) to zero. We find that the boundaries
where acceleration is zero go through each Einstein fixed point parallel the
the $Y$-axis. Therefore, one case remains where a bouncing model exhibits
early and late time acceleration, which is the $0.20<c_{z}<0.32$ case. This
case is shown again in Fig. 15, where we have included the red
$a^{\prime\prime}=0$ boundaries. The $0<x<0.19$ and $0.34<x<0.83$ regions are
accelerating, and the $0.19<x<0.34$ and $0.83<x<1$ regions are decelerating.
Therefore, the bouncing models present between the two CFS curves evolve with
an early- and late-time accelerated phase, and are therefore the models of
interest.

Figure 15: The phase space (top panel) and corresponding potential (bottom
panel) for $0.20<c_{z}<0.32$ as in Fig. 10, with the boundaries between
accelerating and decelerating regions in red. The $0<x<0.19$ and $0.34<x<0.83$
regions are accelerating, and the $0.19<x<0.34$ and $0.83<x<1$ regions are
decelerating.

However, we find that for all of the cases we have presented, the bouncing
behaviour is spoiled when dark matter and radiation become dominant, and turn-
around models occur for trajectories in the phase space near $x=1$. Therefore,
these bouncing models are only relevant when matter and radiation appear after
the bounce, and would need a mechanism such as inflation with a subsequent
period of reheating in order to be feasible. These models also only occur when
they have sufficient curvature. Although we cannot make a quantitative
statement from our qualitative analysis, in principle, these models would be
in tension with observations Aghanim et al. (2018b).

##  V Upper bounds on $z$ and $r$

In light of the results thus far, we introduce an upper bound on the dark
matter and radiation, such that their energy densities cannot reach infinity.
The idea is that matter and radiation are not always present in standard
cosmology, but that they are created in a period of reheating. It is beyond
the scope of our paper to introduce a reheating phase, perhaps through an
interaction between radiation, dark matter and the dark energy component, so
we simply assume that matter and radiation have an upper bound. In the
following, we focus on investigating whether bouncing models can be the
general closed model when matter and radiation appear after the bounce has
occurred, even for models that are close to being spatially flat.

In order to implement an upper bound on matter and radiation, we assume that
their EoS’s are modified by a non-linear term at high energies in a similar
way to the high energy bound for the dark energy component. Physically, these
EoS’s are not particularly meaningful, however this is a simple way of
imposing an energy scale at which matter and radiation appear such that the
bounce is not spoiled. As far as our qualitative analysis is concerned, the
specific physical mechanism for having an upper bound for matter and radiation
is not that relevant. We therefore re-write Eq. (11) and Eq. (12) as

$z^{\prime}=-\frac{3Yz(1-\frac{z}{z_{*}})}{\sqrt{1-Y^{2}}}$ (68)

$r^{\prime}=-\frac{4Yr(1-\frac{r}{r_{*}})}{\sqrt{1-Y^{2}}}\,,$ (69)

where $z_{*}$ and $r_{*}$ are the characteristic energy scales of dark matter
and radiation, respectively. Should we model a proper inflationary era, we
could set $z_{*}$ and $r_{*}$ at the energy scale of reheating. Here, we will
simply assume that the upper bounds imposed on matter and radiation are at the
highest possible energy scale, coinciding with the high energy cosmological
constant of the dark energy, somewhere between the Planck scale and that of
inflation. This assumption does not qualitatively change the dynamics, and is
the most robust way of investigating the survival of the bounce, therefore we
set $z_{*}=r_{*}=1$.

The effective EoS’s for dark matter and radiation can be found by equating
$z^{\prime}$ (68) and $r^{\prime}$ (69) to the continuity equation. We find

$w_{z}=-z$ (70)

$w_{r}=\frac{1}{3}-\frac{4}{3}r\,.$ (71)

For $z>1/3$ and $r>1/2$, we find their respective effective EoS’s are less
than -1/3, and so the dark matter and radiation can contribute to acceleration
in these regions of the phase space. We find that the Raychaudhuri equation
for $Y$, the compact variable representing the Hubble expansion scalar, now
becomes

$Y^{\prime}=-Y^{2}(1-Y^{2})^{\frac{1}{2}}-\frac{(1-Y^{2})^{\frac{3}{2}}}{6}[z(1-3z)\\\
+2r(1-2r)-3\mathcal{R}+x(1+3\mathcal{R})-3x^{2}]\,,$ (72)

and $x^{\prime}$ is still as in Eq. (29). As previously, we project the 4-D
dynamics onto the $x$-$Y$ plane. Integrating $z$ (68) with respect to $x$
(29), we find

$z(x)=\frac{c_{z}\left(\frac{x-\mathcal{R}}{1-x}\right)^{\frac{1}{1-\mathcal{R}}}}{1+c_{z}\left(\frac{x-\mathcal{R}}{1-x}\right)^{\frac{1}{1-\mathcal{R}}}}\,,$
(73)

and integrating $r$ (69) with respect to $x$ (29), we obtain

$r(x)=\frac{c_{r}\left(\frac{x-\mathcal{R}}{1-x}\right)^{\frac{4}{3(1-\mathcal{R})}}}{1+c_{r}\left(\frac{x-\mathcal{R}}{1-x}\right)^{\frac{4}{3(1-\mathcal{R})}}}\,.$
(74)

We can then calculate the first integral $c_{r}$ in same way as before. Using
the Planck 2018 density parameters Aghanim et al. (2018b), along with our
definition of $\mathcal{R}$ (10) and the low energy cosmological constant
$\rho_{\Lambda}$ (42), we find

$c_{r}=\frac{\mathcal{R}\Omega_{r}}{\alpha\Omega_{\Lambda}-\mathcal{R}\Omega_{r}}\left(\frac{\alpha-\mathcal{R}}{\mathcal{R}(1-\alpha)}\right)^{\frac{4}{3(1-\mathcal{R})}}\simeq
0.00077$ (75)

for a value of $\mathcal{R}=0.05$.

The fixed points for this system are given in Table 5. In this system, there
are no physical singularities, but the fixed points $dS_{2\pm}$ and
$dS_{4\pm}$ are coordinate singularities of the de Sitter spacetime when
represented as an FLRW, see Hawking and Ellis (1973).

The eigenvalues of the fixed points for this system are shown in Table 6, and
the linear stability classification for each point is given in Table 7. As
previously, we cannot include the eigenvalues of the fixed points at $Y=\pm
1$, however we find from the Raychaudhuri equation (72) that the de Sitter
point at $Y=+1$ is a repeller and the de Sitter point at $Y=-1$ is an
attractor along $x=1$. The de Sitter points along $x=\mathcal{R}$ at $Y=\pm 1$
are saddle points, with trajectories moving away from the $Y=+1$ point and
towards the $Y=-1$ point. Similarly to before, the stability character of the
Einstein fixed points depends on the value of the first integral $c_{z}$ for a
fixed $\mathcal{R}$ and $c_{r}$. Similarly to before, the stability character
of the Einstein fixed points depends on the value of the first integral
$c_{z}$ for a fixed $\mathcal{R}$ and $c_{r}$.

In this case, in order to find the Einstein points from the Raychaudhuri
equation (72), we need to find the zeros of the following function

$f(x)=z(x)[1-3z(x)]+2r(x)[1-2r(x)]\\\
-3\mathcal{R}+x(1+3\mathcal{R})-3x^{2}\,.$ (76)

Name | $x$ | $Y$ | $z$ | $r$
---|---|---|---|---
$E$ | $x$ | $0$ | $z$ | $\frac{1}{4}\left(1-\sqrt{1+4(x+z)+12(x\mathcal{R}-\mathcal{R}-x^{2}-z^{2})}\right)$
dS1± | $\mathcal{R}$ | $\pm\sqrt{\frac{\mathcal{R}}{3+\mathcal{R}}}$ | $0$ | $0$
dS2± | $\mathcal{R}$ | $\pm 1$ | $0$ | $0$
dS3± | $1$ | $\pm\frac{1}{\sqrt{2}}$ | $1$ | $1$
dS4± | $1$ | $\pm 1$ | $1$ | $1$
Table 5:  The fixed points of the system when $z$ and $r$ have an upper bound Name | $\lambda_{1}$ | $\lambda_{2}$ | $\lambda_{3}$ | $\lambda_{4}$
---|---|---|---|---
$E$ | 0 | 0 |  $\frac{1}{\sqrt{6}}\\{-1+18x^{3}-z+9z^{2}(2z-1)+27\mathcal{R}-9x^{2}(1+3\mathcal{R})-x$ $+[1+2z(3z-1)+2x(3x-1)]\sqrt{1+4z(1-3z)-12\mathcal{R}+4x(1-3x+3\mathcal{R})}\\}^{\frac{1}{2}}$ | $-\lambda_{3}$
$dS_{1\pm}$ | $\mp 4\sqrt{\frac{\mathcal{R}}{3}}$ | $\mp\sqrt{3\mathcal{R}}$ | $\mp 2\sqrt{\frac{\mathcal{R}}{3}}$ | $\mp\sqrt{3\mathcal{R}}(1-\mathcal{R})$
$dS_{3\pm}$ | $\pm 4$ | $\pm 3$ | $\mp 2$ | $\pm 3(1-\mathcal{R})$
Table 6:  The eigenvalues for the fixed points of the system when $z$ and $r$ have an upper bound Name | Stability Character
---|---
$E$ | Centre, Saddle or Cusp
$dS_{1+}$ | Attractor
$dS_{1-}$ | Repellor
$dS_{2\pm}$ | Saddle
$dS_{3\pm}$ | Saddle
$dS_{4+}$ | Repellor
$dS_{4-}$ | Attractor
Table 7:  The linear stability character for the fixed points of the system
when $z$ and $r$ have an upper bound

We need to solve numerically for the values $x_{E}$, such that $f(x_{E})=0$.
Taking the limit of $f(x)$ when $x\rightarrow\mathcal{R}$, we find
$f(x)\rightarrow-2\mathcal{R}$, and when $x\rightarrow 1$, we find
$f(x)\rightarrow-6$. Therefore, this system does not necessarily admit any
Einstein points between $x=\mathcal{R}$ and $x=1$ as both limits are negative.
Fig. 16 shows the range of Einstein points which can be admitted by the system
for different values of $c_{z}$.

Figure 16: $f(x)$ with varying values of $c_{z}$ to show the different number
of Einstein points which can be admitted by the system. Starting from the left
we have $c_{z}>11.0$ (orange), $c_{z}\simeq 11.0$ (green), $0.25<c_{z}<11.0$
(magenta), $c_{z}\simeq 0.25$ (blue) and $c_{z}<0.25$ (purple). When
$c_{z}<0.25$ or $c_{z}>11.0$ no Einstein points are admitted. One Einstein
point is admitted in the limiting cases when $c_{z}\simeq 0.25$ and
$c_{z}\simeq 11.0$, and two Einstein points are admitted when
$0.25<c_{z}<11.0$.

In total, we find three different cases for the dynamics, except here more
than one range of $c_{z}$ can give the same qualitative dynamics. Increasing
$c_{z}$ increases the contribution of the $z(1-3z)$ term to $f(x)$ (76) until
$z=1/3$. Once this point is reached, increasing $c_{z}$ then increases a
negative contribution from the $z(1-3z)$ term. There are two limiting cases
where one Einstein point is admitted which have the same qualitative dynamics
that are shown by the green and blue curves in Fig. 16. As before, to find the
values of $c_{z}$ in the limiting cases, we take the first derivative of
$f(x)$ and solve for $c_{z}$ when $df(x)/dx=f(x)=0$, and find $c_{z}\simeq
0.25$ and $c_{z}\simeq 11.0$. Both $c_{z}<0.25$ and $c_{z}>11.0$ admit the
same qualitative dynamics, where no Einstein points are admitted. In the range
$0.25<c_{z}<11.0$ two Einstein points are admitted.

The open and flat models evolve in the same way for each case. Expanding
(contracting) open models evolve between two de Sitter fixed points, from (to)
an open geometry to (from) a flat geometry. Expanding (contracting) flat
models evolve between two de Sitter fixed points along the FFS. The dynamical
behaviour of the closed models changes depending on the value of $c_{z}$. In
each case, we find that the closed models avoid a singularity as $U\rightarrow
0$ as $x\rightarrow 1$, forming a potential barrier. Therefore, there are no
turn-around models that evolve between singularities, and all closed models
admit a bounce. We present these cases in the following subsections, and note
that qualitatively the dynamics is the same as the sub-manifolds in Section
III when only dark energy is present.

###  V.1 No Einstein points

The phase space for the system when $c_{z}<0.25$ is shown in Fig. 17, which is
qualitatively the same as when $c_{z}>11.0$. As before, the green outermost
thick curve is the flat Friedmann Separatrix (FFS), which here has the general
form

$\frac{Y^{2}}{1-Y^{2}}=\frac{x}{3}+\frac{z}{3}+\frac{r}{3}\,.$ (77)

The inner most curve is the closed Friedmann Separatrix (CFS), which has the
general form

$\frac{Y^{2}}{1-Y^{2}}=\frac{x}{3}+\frac{z}{3}+\frac{r}{3}-\frac{1}{a^{2}}\left(\frac{x_{E}}{3}+\frac{z_{E}}{3}+\frac{r_{E}}{3}\right)\,.$
(78)

$z$ and $r$ are given by Eq. (73) and Eq. (74), respectively, and $a$ is given
by Eq. (53). As previously, the horizontal lines along $x=\mathcal{R}$ and
$x=1$ are the de Sitter lines where $x^{\prime}=z^{\prime}=r^{\prime}=0$.

In this case, no Einstein points are admitted, and all closed models bounce,
evolving between two de Sitter points. Therefore, when matter and radiation
have an upper bound, the bounce is no longer spoiled and a singularity is
avoided. When $c_{z}<0.25$, dark energy is the dominant component for most of
the phase space, until $x=1$ is approached where dark matter and radiation
become comparable to the dark energy. When $c_{z}>11.0$, dark matter is the
dominant component for most of the phase space, except sufficiently close to
$x=1$ where all three components are comparable, and close to $x=\mathcal{R}$
where dark energy becomes dominant. All trajectories are always accelerating
in this case.

Figure 17: $c_{z}<0.25$ case, which is qualitatively the same as the
$c_{z}>11.0$ case. Top panel: the projection of the full 4-D dynamics on the
2-D $x$-$Y$ plane. Bottom panel: the corresponding potential as in Eq. (20),
where $a$ is given by Eq. (53), $z$ by Eq. (73) and $r$ by Eq. (74).
Trajectories of the same colour in the two panels correspond to each other. No
Einstein points are admitted, and all closed models bounce, which always
accelerate. Figure 18: $c_{z}\simeq 0.25$ case, which is qualitatively the
same as the $c_{z}\simeq 11.0$ case. Top panel: the projection of the full 4-D
dynamics on the 2-D $x$-$Y$ plane. Bottom panel: the corresponding potential
as in Eq. (20), where $a$ is given by Eq. (53), $z$ by Eq. (73) and $r$ by Eq.
(74). Trajectories of the same colour in the two panels correspond to each
other. One Einstein point exists at $x_{E}\simeq 0.23$ (cusp) corresponding to
the horizontal point of inflection in the potential. All closed models bounce,
and these are always accelerating.

###  V.2 One Einstein point

The limiting case where $c_{z}\simeq 0.25$ and one Einstein point is admitted
is shown in Fig. 18. The qualitative dynamics in this case is the same when
$c_{z}\simeq 11.0$. One Einstein point is admitted at $x_{E}\simeq 0.23$
(cusp) which is part of the CFS, and corresponds to the horizontal point of
inflection in the potential. All closed models bounce, evolving between two de
Sitter fixed points. When $c_{z}\simeq 0.25$, dark energy is dominant for most
of the phase space, and when $c_{z}\simeq 11$, dark matter is mostly dominant,
except close to $x=\mathcal{R}$ where dark energy is dominant. All three
components become comparable near $x=1$. Here, all trajectories are always
accelerating.

Figure 19: $0.25<c_{z}<11.0$ case, top panel: the projection of the full 4-D
dynamics on the 2-D $x$-$Y$ plane. Bottom panel: the corresponding potential
as in Eq. (20), where $a$ is given by Eq. (53), $z$ by Eq. (73) and $r$ by Eq.
(74). Trajectories of the same colour in the two panels correspond to each
other. Two Einstein points exist at $x_{E}\simeq 0.15$ (saddle), which
corresponds to a local maximum of the potential, and at $x_{E}\simeq 0.31$
(centre) corresponding to a local minimum of the potential. Trajectories
around the $x_{E}\simeq 0.31$ fixed point within the CFS are cyclic, and all
other closed models bounce. Bouncing models outside the CFS evolve with an
early- and late-time acceleration, connected by a period of deceleration.

###  V.3 Two Einstein points

The case where two Einstein points are admitted is shown in Fig. 19. The
Einstein point at $x_{E}\simeq 0.15$ (saddle) corresponds to a local maximum
of the potential, and is part of the CFS. There are bouncing models which
evolve between two de Sitter points within the CFS, and cyclic models around
the Einstein point at $x_{E}\simeq 0.31$ (centre), which corresponds to a
local minimum of the potential. These bouncing models always accelerate, and
the cyclic models evolve with an early-time acceleration and a late-time
deceleration. Outside of the CFS, all closed models bounce, evolving between
two de Sitter points. These models evolve with an early- and late-time
accelerated expansion, connected by a decelerating period. If the value of
$c_{z}$ is closer to 0.25, most of the phase space will be dark energy
dominated, however if its value is closer to 11.0, then dark matter will be
mostly dominant.

Figure 20: The phase space (top panel) and corresponding potential (bottom
panel) for $0.25<c_{z}<11.0$ as in Fig. 19, with the boundaries between
accelerating and decelerating regions in red. The $0.15<x<0.31$ region of the
phase space is decelerating, and the $x<0.15$ and $x>0.31$ regions are
accelerating.

###  V.4 Acceleration

As before, we calculate where the acceleration is zero to find the boundaries
between accelerating and decelerating regions of the phase space. In this
system, the acceleration equation is given by

$\frac{a^{\prime\prime}}{a}=z(1-3z)+2r(1-2r)-3\mathcal{R}+(1+3\mathcal{R})x-3x^{2}\,.$
(79)

The case of interest in this new system is $0.25<c_{z}<11.0$, in which two
Einstein points are admitted. This case with the boundaries between
accelerating and decelerating regions of the phase space can be seen in Fig.
20. The two curves through the Einstein fixed points parallel to the $Y$-axis
denote where acceleration is zero. The $0.15<x<0.31$ region of the phase space
is decelerating, and the $x<0.15$ and $x>0.31$ regions are accelerating.
Therefore, bouncing models outside the CFS evolve with early- and late-time
acceleration, with a decelerating phase in between, and are therefore the
models of interest.

Finally, using $a^{\prime}=ay$ and $a^{\prime\prime}/a=y^{\prime}+y^{2}$ in
Eq. (79) (where each variable now is a function of $a$), in Fig. 21 we
illustrate this model with the phase space for the scale factor $a$
(normalised to $a_{o}=1$ today) and non-compactified Hubble expansion scalar
$y$. Again, we see in this plot that at high energies, i.e. small $a$, all
models have an accelerated phase, followed by a decelerated one and a final
accelerated one at recent times, i.e. when $a\rightarrow 1$.

Figure 21: The phase space (top panel) and corresponding potential (bottom
panel, where now $U=U(a)$) for $0.25<c_{z}<11.0$ in terms of the scale factor
$a$ and non-compactified Hubble expansion scalar $y$. This plot is equivalent
to Fig. 20, with the boundaries separating the accelerating and decelerating
regions shown in red. The region in between the two Einstein fixed points is
decelerating, and the regions outside are accelerating.

###  V.5 A note on $\mathcal{R}$

As we mentioned previously, in order to alleviate the old cosmological
constant problem we would set $10^{-120}<\mathcal{R}<10^{-60}$, however purely
for the purpose of illustrating the full range of dynamics in readable phase
plots, we have set $\mathcal{R}=0.05$. If we had set
$10^{-120}<\mathcal{R}<10^{-60}$, only one case for the dynamics would remain
in each system. For both systems we have considered, we find from Eq. (63) and
Eq. (75) that for $10^{-120}<\mathcal{R}<10^{-60}$, we obtain
$10^{16}<c_{r}<10^{36}$. The effect of increasing radiation has qualitatively
the same effect as increasing matter.

For the first system of equations (27) - (30) where $z$ and $r$ are unbounded,
we find the only dynamical case left, regardless of the value of $c_{z}$, is
as in Fig. 14. The Einstein point would be pushed further towards
$x=\mathcal{R}$ and the closed models are mostly turn-around models which
always decelerate. The bouncing models in this case are always accelerating,
which is not the evolution we require for a realistic model.

For the second system of equations (29), (68), (69) and (72) where $z$ and $r$
have an upper bound, we find that the only dynamical case is as in Fig. 17.
All closed models bounce, however they are always accelerating, which again is
not the evolution we require.

As a final remark, we note that in order for a bouncing model to be
theoretically robust and realistic, the value of $\mathcal{R}$ would need to
satisfy the bounds $10^{-120}<\mathcal{R}<10^{-60}$, producing a bounce with
an early acceleration, followed by a decelerated era, and finally a late-time
acceleration as in Fig. 20. The decelerated era should be a standard matter
and radiation dominated phase, in order to satisfy observational constraints,
but a contribution from the homogeneous dark energy remains a possibility.

##  VI Conclusions

In this paper, we have studied the dynamics of FLRW models containing dark
matter, radiation and dark energy with a quadratic EoS. This is an extension
of Ananda and Bruni (2006a), who studied a more general quadratic EoS in the
high, low and full energy regimes without the inclusion of matter and
radiation, and found bouncing and cyclic models were possible with certain
combinations of parameters. A quadratic EoS is the simplest nonlinear EoS and
the qualitative analysis of the dynamics that it generates serves as guidance
for more complicated nonlinear models; its study finds motivations in brane
models, k-essence and loop quantum cosmology Shiromizu et al. (2000); Bridgman
et al. (2002); Gergely and Maartens (2002); Langlois (2003); Maartens and
Koyama (2010); Vandersloot (2005); Parisi et al. (2007); Scherrer (2004);
Giannakis and Hu (2005) (see also Ananda and Bruni (2006a, b) and references
therein). We have restricted the EoS here so that the dark energy evolves
between a high energy effective cosmological constant, which could be of order
of the Planck energy, and a low energy effective cosmological constant close
to the observed dark energy density today. Our aim was to investigate the
effect of matter and radiation on the bouncing and cyclic models, to find
whether these were still possible in this more realistic scenario. In
particular, our focus was on all closed models avoiding singularities, instead
having a bounce or cycles. Because the EoS we consider is barotropic,
naturally the evolution of the system is adiabatic, so that in phase space the
expansion is a mirror image of the contraction, as in Ananda and Bruni (2006a,
b); Bruni et al. (2022).

In Section II, we have presented the system of equations for compactified
variables, with the energy densities of dark matter and radiation satisfying
the standard energy conservation equations, and as such are in principle
unbounded and able to become infinite. In Section III, we present the sub-
manifolds of the system. In Section IV, we found one case of interest, shown
in Fig. 15, where closed bouncing models can evolve with early and late time
acceleration. However, we found that this bouncing behaviour is spoiled when
matter and radiation become dominant for certain values of the initial
conditions, and instead the system evolves with a turn-around between two
singularities. In comparison to Fig. 5 (with the same parameters but no matter
and radiation), in which all closed models bounce, we concluded that these
closed models could only be viable, i.e. always have a bounce independently of
the initial conditions, if matter and radiation are only present at energies
below the energy scale of the bounce. This would require a process such as
reheating during the expansion phase after the bounce.

In light of this, we modelled the process by simply introducing upper bounds
on dark matter and radiation in Section V, in order to avoid any models
becoming singular. To keep things simple, with no need of extra parameters, we
set this upper energy scale in line with the high energy effective
cosmological constant of the dark energy component. In this case, all closed
models have a bounce. The case of interest, shown in Fig. 20, admits bouncing
models that evolve with early- and late-time acceleration which are not
spoiled by the inclusion of matter and radiation, as they only appear at the
same energy scale as the dark energy.

However, we noted that for the specific case of a quadratic EoS, this case of
interest is not possible with a model that can alleviate the old cosmological
constant problem Weinberg (1989); Straumann (1999); Ellis et al. (2012). For a
model where the two cosmological constants differ by 60 to 120 orders of
magnitude, we found that the only case for the dynamics of the system is as in
Fig. 17. All closed models bounce, however they always accelerate and so the
evolution is undesirable.

Overall, our qualitative analysis shows that a model with nonlinear dark
energy and radiation and matter can be a realistic scenario like the one shown
in Fig. 20 and 21, with an early- and a late-time acceleration, where all
closed models have a bounce, assuming that matter and radiation appear only
after the bounce through a process such as reheating. However, while the
simple specific model for dark energy considered here helps our understanding
of precisely what the qualitative realistic scenario should be, it does not
give the right quantitative behaviour. Learning from the present analysis, in
future work we will consider models for dark energy, or for an interacting
vacuum Bruni et al. (2022), giving the qualitative behaviour of Fig. 20 and
21, but with a realistic quantitative behaviour that will be worth
constraining with data. Another interesting possibility that we leave for a
future analysis will be to include a reheating phase through an interaction
term between the dark energy and radiation components, with a possible role of
dark matter.

## Acknowledgements

We thank Johannes Noller and David Wands for useful discussions and comments.
This work has been supported by UK STFC grants ST/S000550/1 and ST/T506345/1.

## References

  * Aghanim et al. (2018a) N. Aghanim, Y. Akrami, F. Arroja, M. Ashdown, J. Aumont, C. Baccigalupi, M. Ballardini, A. J. Banday, R. B. Barreiro, N. Bartolo, et al. (2018a), URL http://arxiv.org/abs/1807.06205.
  * Aghanim et al. (2018b) N. Aghanim, Y. Akrami, M. Ashdown, J. Aumont, C. Baccigalupi, M. Ballardini, A. J. Banday, R. B. Barreiro, N. Bartolo, S. Basak, et al. (2018b), URL http://arxiv.org/abs/1807.06209.
  * Penrose (1965) R. Penrose, Phys. Rev. Lett. 14, 57 (1965).
  * Hawking (1965) S. W. Hawking, Phys. Rev. Lett. 15 (1965).
  * Hawking (1966a) S. W. Hawking, Proc. Roy. Soc. Lond. A 294, 511 (1966a).
  * Hawking (1966b) S. W. Hawking, Proc. Roy. Soc. Lond. A 295, 490 (1966b).
  * Hawking (1967) S. W. Hawking, Proc. Roy. Soc. Lond. A 300, 187 (1967).
  * Hawking and Penrose (1970) S. Hawking and R. Penrose, Proc. Roy. Soc. Lond. A 314, 529 (1970).
  * Hawking and Ellis (1973) S. Hawking and G. Ellis, _The large scale structure of space-time_ (Cambridge University Press, 1973).
  * Wald (1984) R. M. Wald, _General Relativity_ (The University of Chicago Press, 1984).
  * Joshi (1993) P. Joshi, _Global Aspects in Gravitation and Cosmology_ (Oxford University Press Inc., 1993).
  * Hawking and Israel (1979) S. Hawking and W. Israel, _General Relativity: An Einstein Centenary Survey_ (Cambridge University Press, 1979).
  * Gibbons et al. (2003) G. W. Gibbons, E. P. S. Shellard, and S. J. Ranking, eds., _The Future of Theoretical Physics and Cosmology_ (Cambridge University Press, 2003).
  * Ashtekar et al. (2015) A. Ashtekar, B. K. Berger, J. Isenberg, and M. MacCallum, eds., _General Relativity and Gravitation: A Centennial Perspective_ (Cambridge University Press, 2015).
  * Gasperini and Veneziano (1993) M. Gasperini and G. Veneziano, Astropart. Phys. 1, 317 (1993).
  * Khoury et al. (2001) J. Khoury, B. A. Ovrut, P. J. Steinhardt, and N. Turok, _The ekpyrotic universe: Colliding branes and the origin of the hot big bang_ (2001).
  * Kallosh et al. (2001) R. Kallosh, L. Kofman, and A. Linde, Phys. Rev. D. 64, 123523 (2001).
  * Steinhardt and Turok (2002) P. J. Steinhardt and N. Turok, Physical Review D 65, 126003 (2002).
  * Novello and Bergliaffa (2008) M. Novello and S. E. P. Bergliaffa, Phys. Rept. 463, 127 (2008).
  * Wands (2008) D. Wands, Adv. Sci. Lett. 2, 194 (2008).
  * Maier et al. (2012) R. Maier, S. Pereira, N. Pinto-Neto, and B. B. Siffert, Phys. Rev. D. 85, 023508 (2012).
  * Maier et al. (2013) R. Maier, N. Pinto-Neto, and I. D. Soares, Phys. Rev. D. 87, 043528 (2013).
  * Battefeld and Peter (2014) D. Battefeld and P. Peter, Phys. Rept. 571, 1 (2014).
  * Cai et al. (2017) Y. Cai, Y. Wan, H.-G. Li, T. Qiu, and Y.-S. Piao, JHEP 1, 90 (2017).
  * Brandenberger and Peter (2017) R. Brandenberger and P. Peter, Found. Phys. 47, 797 (2017).
  * Cai and Piao (2017) Y. Cai and Y.-S. Piao, JHEP 9, 27 (2017).
  * Ijjas and Steinhardt (2018) A. Ijjas and P. J. Steinhardt, Class. Quant. Grav. 35 (2018).
  * Ijjas and Steinhardt (2019) A. Ijjas and P. J. Steinhardt, Phys. Lett. B 795, 666 (2019).
  * Ganguly and Bruni (2019) C. Ganguly and M. Bruni, Phys. Rev. Lett. 123, 201301 (2019).
  * Ilyas et al. (2020) A. Ilyas, M. Zhu, Y. Zheng, Y.-F. Cai, and E. N. Saridakis, JCAP 9, 2 (2020).
  * Battista (2021) E. Battista, Class. Quant. Grav. 38, 195007 (2021).
  * Zhu et al. (2021) M. Zhu, A. Ilyas, Y. Zheng, Y.-F. Cai, and E. N. Saridakis, JCAP 11, 45 (2021).
  * Riess et al. (1998) A. G. Riess, A. V. Filippenko, P. Challis, A. Clocchiatti, A. Diercks, P. M. Garnavich, R. L. Gilliland, C. J. Hogan, S. Jha, R. P. Kirshner, et al., Astron. J. 116, 1009 (1998).
  * Perlmutter et al. (1998) S. Perlmutter, G. Aldering, G. Goldhaber, R. A. Knop, P. Nugent, P. G. Castro, S. Deustua, S. Fabbro, A. Goobar, D. E. Groom, et al., Astrophys. J. pp. 565–586 (1998).
  * Weinberg (1989) S. Weinberg, Rev. Mod. Phys. 61 (1989).
  * Straumann (1999) N. Straumann, Eur. J. Phys. 20, 419 (1999).
  * Ellis et al. (2012) G. Ellis, R. Maartens, and M. MacCallum, _Relativistic Cosmology_ (Cambridge University Press, 2012).
  * Lemaître (1934) G. Lemaître, Proceedings of the National Academy of Sciences of the United States of America 20, 12 (1934).
  * Wands et al. (2012) D. Wands, J. De-Santiago, and Y. Wang, Class. Quant. Grav. 29 (2012).
  * Bruni et al. (2022) M. Bruni, R. Maier, and D. Wands, Phys. Rev. D 105, 063532 (2022).
  * Ananda and Bruni (2006a) K. N. Ananda and M. Bruni, Phys. Rev. D 74, 023523 (2006a).
  * Shiromizu et al. (2000) T. Shiromizu, K. I. Maeda, and M. Sasaki, Phys. Rev. D. 62, 024012 (2000).
  * Bridgman et al. (2002) H. A. Bridgman, K. A. Malik, and D. Wands, Phys. Rev. D 65, 043502 (2002).
  * Gergely and Maartens (2002) L. Gergely and R. Maartens, Class. Quant. Grav. 19, 213 (2002).
  * Langlois (2003) D. Langlois, Astrophys. Space Sci. 283, 469 (2003).
  * Maartens and Koyama (2010) R. Maartens and K. Koyama, Living Rev. Relativity 13 (2010).
  * Vandersloot (2005) K. Vandersloot, Phys. Rev. D. 71, 103506 (2005).
  * Parisi et al. (2007) L. Parisi, M. Bruni, R. Maartens, and K. Vandersloot, Class. Quant. Grav. 24, 6243 (2007).
  * Scherrer (2004) R. J. Scherrer, Phys. Rev. Lett. 93, 011301 (2004).
  * Giannakis and Hu (2005) D. Giannakis and W. Hu, Phys. Rev. D. 72, 063502 (2005).
  * Ananda and Bruni (2006b) K. N. Ananda and M. Bruni, Phys. Rev. D 74, 023524 (2006b).
  * Burkmar (2023) M. Burkmar (2023), URL https://github.com/MollyBurkmar1/Dynamics_Quadratic_EoS.git.
  * Hogg and Bruni (2021) N. B. Hogg and M. Bruni, Mon. Not. Roy. Astron. Soc. 511, 4430 (2021).
  * Arrowsmith and Place (1982) D. Arrowsmith and C. Place, _Ordinary Differential Equations_ (Chapman and Hall, 1982).
  * Arrowsmith and Place (1992) D. K. Arrowsmith and C. M. Place, _Dynamical systems: differential equations, maps and chaotic behaviour_ (Chapman and Hall, 1992).
  * Strogatz (2015) S. Strogatz, _Nonlinear dynamics and chaos: with applications to physics, biology, chemistry, and engineering_ (Westview Press, 2015), 2nd ed.
  * Boehmer and Chan (2014) C. G. Boehmer and N. Chan (2014), URL https://arxiv.org/abs/1409.5585.