Bouncing cosmology from nonlinear dark energy with two cosmological constants

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 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. We find that bounces are always the general closed model, 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 Λ-cold-dark-matter (ΛCDM) model has provided a successful framework to describe the history of the Universe, and is consistent with a wealth of observational data [1,2]. 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 [3][4][5][6][7][8][9][10], however their current interpretation is that they represent points where General Relativity breaks down [11][12][13][14]. 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 [15][16][17][18][19][20][21][22][23][24][25][26][27][28][29][30][31][32].
Another problem facing ΛCDM is that of Λ itself. Observations have provided strong evidence for the accelerated expansion of our Universe [33,34], which requires a component that violates the strong energy condition [9], known as dark energy. The cosmological constant Λ is the simplest explanation we have for dark energy, however the observed value of Λ 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) [35][36][37].
Taking into account a dynamic dark energy component beyond Λ, it is worth considering whether a bounce can be produced. A non-interacting vacuum is equivalent to Λ in General Relativity [38], so a simple extension is to consider a vacuum that interacts with other components, thereby becoming dynamical [39]. 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 [40].
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 [41]. In addition, a quadratic EoS finds motivations from e.g. brane-world models [42][43][44][45][46], loop quantum cosmology [47,48] and kessence [49,50] (see also [41,51] 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 arXiv:2302.03710v2 [gr-qc] 18 May 2023 cosmological constants and not on the specific model used to produce them.
In this paper, we analyse a subset of the quadratic EoS from [41,51] 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 [41,51] 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 Λ. 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 [52].
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 Λ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 postbounce process such as reheating, when matter and radiation are created. In this spirit, in Section V we reconsider 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 → 0 and w r → 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 latetime 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πG = c = 1.

II. COSMOLOGICAL DYNAMICAL SYSTEM
A. 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 tȯ where ρ x and P x are the energy density and isotropic pressure of dark energy, ρ m is the energy density of dark matter and ρ 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 iṡ We assume that the dark energy EoS is barotropic, (1) has two fixed points (ρ x = 0), corresponding to two effective cosmological constants: It follows from the assumption P x (ρ x ) ≥ −ρ x that the dark energy density decreases/increases during a period of expansion/contraction. Then, during expansion (H > 0), we take ρ Λ to represent the attractor at low energy close to the dark energy density we observe today, and ρ * the repellor at high energy; their roles are inverted during contraction. We assume ρ * 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 [41,51], which is the simplest non-linear case to study and, as we said in the introduction (see also [41,51]), has other motivations from various scenarios [42][43][44][45][46][47][48][49][50]. With respect to [41,51], we restrict parameters such that the linear term is always positive, α > 0, the quadratic term is always negative, β = −1/ρ * , and the constant pressure term is negative, P 0 = −ρ Λ . With this choice of parameters our EoS is This secures the two effective cosmological constants in Eq. (5), as now the energy conservation equation of the dark energy (1) can be written aṡ 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 (ρ x ) + ρ x > 0 for ρ Λ < ρ x < ρ * , i.e. between the two fixed points (5) satisfying P x (ρ x ) = −ρ x , as these conditions are enough to describe any monotonically decreasing ρ x between ρ * and ρ Λ during expansion, see e.g. the EoS in [53].
The effective EoS parameter w x = P x /ρ x for the dark energy is as 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 where k is the curvature and a is the cosmic scale factor, connected to the Hubble expansion function through the expression H =ȧ/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.

B. 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 [41,51], we define these variables as: (10) The variable R is the ratio of the low energy effective cosmological constant ρ Λ to the high energy effective cosmological constant ρ * 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 [R, 1]. The variable y is the normalised Hubble parameter, z the normalised dark matter energy density, r the normalised radiation energy density and η 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 where the primes indicate differentiation with respect to η. The Raychaudhuri equation (4) can now be expressed as and the Friedmann equation (9) becomes The effective EoS parameter for the dark energy (8) in dimensionless variables becomes therefore, the necessary condition for acceleration is We can also express the evolution of the cosmic scale factor a using our dimensionless variables a = ay .
The Friedmann equation (15) can then be rearranged for a kinetic term a 2 /2, a potential U and total energy E where U is given by and E by 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) ∼ a −3 and r(a) ∼ a −4 , respectively. For x (13), we obtain where the constant of integration is In the following, however, we will express z, r and a as functions of x.

C. Compactified Variables
In our system, y takes values in the range (−∞, ∞), and z and r in the range [0, ∞). x is already limited to the range [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: and the Friedmann equation (15) can be expressed as Now, setting a 0 = 1, we find Finally, the potential U (20) in terms of compact variables is given by

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 [54][55][56]. 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 Sub-manifold Fixed Points Name Eigenvalues Stability Character  The sub-manifolds of the system with their fixed points, eigenvalues and stability character. E denotes an Einstein universe, dS± an expanding (+) or contracting (-) de-Sitter universe, and S± a singularity with infinite expansion (+) or contraction (-). The fixed points at Y = ±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.
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 nonzero 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 eigenval-ues 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 = Z = R = 0. Y can vary for positively and negatively curved de Sitter models. For flat models, Y is constant, Y = 0, giving rise to de Sitter fixed points. These occur at the effective cosmological constants x = 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 = 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 submanifolds of the system, which simplifies the dynamics by reducing the dimensions of the system [57]. 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 I. The fixed points labelled E represent a static Einstein universe and the fixed points labelled dS ± represent spatially flat expanding (+) and contracting (-) de-Sitter models. S ± denotes singularities with infinite expansion (+) or contraction (-) and infinite energy density. Note we cannot compute the eigenvalues of any fixed point at Y = ±1, as the Jacobian becomes singular and we therefore cannot Taylor expand around them. However, from the Raychaudhuri equation (30), we find that Y is negative along the de Sitter lines x = R and x = 1 around the fixed points at Y = ±1.

A. ΛCDM Dynamics
We begin with the sub-manifolds with an effective cosmological constant and so are analogous to ΛCDM, except here x = 1 and x = R are asymptotic values rather than true constants. In each of these cases, x = 0 (29).

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 In this case, the Friedmann equation (31) becomes (36) and the potential (33) is For the sub-manifolds where x is a constant, we express a in terms of Z. Integrating a (18) with respect to Z (27), we find We can then solve for the fixed points of the system. There are de Sitter fixed points where Z = R = 0, and Y = ±1/2 and Y = ±1. The fixed points at Y = ±1 and at Z = R = 0 are coordinate singularities of the de Sitter spacetime when represented as an FLRW, see [9]. At Z = R = 1, (27) and (28) have fixed points, and at Y = ±1 (35) has fixed points. Together, the asymptotic points at Z = R = 1 and Y = ±1 represent singularities, with infinite expansion (+) or contraction (-) and infinite energy density. At the Einstein fixed point, Y = 0 and the condition Y = 0 reduces Eq. (35) to a constraint between Z and R: 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 where H 0 , Ω r and Ω m can be fixed using Planck values [2]. This results in (27) and (28) admitting a first integral c rz , which can be used to write a relation between R and Z where c rz depends on Ω m , Ω r and ρ * . 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 ρ * = ρ Λ /R, where the future asymptotic value ρ Λ is by definition some fraction of our current observed dark energy density, and 0 < α < 1. For the purpose of our qualitative analysis, we fix α = 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 < 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 R. We therefore choose R = 0.05 in order to illustrate the evolution of our variables, given that the value of R will not qualitatively change the dynamics. We therefore find where we have used Ω Λ = 0.6889 and Ω m = 0.3111, and calculated Ω r = 9.1824 × 10 −5 using Ω m and the redshift of matter-radiation equality z = 3387 [2]. 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 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 where again here x = x E = 1 and the relations in Eq. (38) and Eq. (41) apply. 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 = Z = R = 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 > 0) to the left of the Einstein point. Therefore in a ΛCDM cosmology, such as the one we are illustrating here, bouncing models are always accelerating, and re-collapsing models are never accelerating.
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.

x = R
Next we consider the case where x = 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 The Friedmann equation (31) (47) and the potential (33) is The dynamics of this system is qualitatively the same as in Fig. 1. The singularities are present at Z = R = 1 and Y = ±1, and there are de Sitter points at Y = R/(3 + R) and Y = ±1 along Z = R = 0. At the Einstein fixed point, the Hubble expansion variable (46) reduces to 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 = R/(1 + R), and in the R = 0 case at Z E = 2R/(1 + 2R).

B. Dynamic Dark Energy
We now consider the sub-manifolds with dynamic dark energy so that x is not fixed, and therefore x = 0 in Eq. (29).
We consider the sub-manifold where matter and radiation are not present, and there is only dark energy. In this case Z = 0 in Eq. (27) and R = 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 (50) When Z = R = 0, the Friedmann equation (31) becomes and the potential (33) is 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 (18) with respect to x (13), we find i.e. the inverse of Eq. (22), where c a is as in Eq. (23).
There are de Sitter points at Y = R/(3 + R) and Y = ±1 along x = R, and at Y = ±1/2 and Y = ±1 along x = 1. The fixed points at Y = ±1 and at x = R and x = 1 are coordinate singularities of the de Sitter spacetime when represented as an FLRW, see [9]. At the Einstein fixed point, Y = 0, Eq. (50) reduces to Rearranging for x E , we find Given the condition 0 < 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 The R > R M case is shown in Fig. 3, where the system admits no Einstein fixed points. The two vertical black lines along x = 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 R = 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 R < R M is shown in Fig. 5. Two Einstein points are admitted: the fixed point at x E 0.07 is a saddle point, which is part of the CFS (the black loop in Fig. 5), and the other at x E 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 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).

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.
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 where the constant of integration is Here, to express c z , we set z 0 using the Friedmann equation (31) when k = 0: We can then express z 0 as Then to express r in terms of x, from Eq. (12) and Eq. (13) we find where the constant of integration is We set c r in a similar way to c rz , using the density parameters [2], and our definitions of R (10) and ρ Λ (42). Keeping our choice of R = 0.05, we then find 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  (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.
de Sitter lines along x = 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.

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 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 II. As before, E represents a static Einstein universe, dS ± represents spatially flat expanding (+) and contracting (-) de-Sitter models, and S ± represents singu-larities with infinite expansion (+) or contraction (-) and infinite energy density. Note that the dS 2± fixed points are coordinate singularities of the de Sitter spacetime when represented as an FLRW, see [9]. Table III shows the eigenvalues of the fixed points for this system, and the linear stability classification for each point is given in Table IV. Again, we cannot include the eigenvalues of the fixed points at Y = ±1 in Table III, 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 = ±1 are saddle points, with trajectories moving away from the Y = +1 point and towards the Y = −1 point along x = 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. The fixed points of the full system. E denotes an Einstein universe, dS± an expanding (+) or contracting (-) de-Sitter universe, and S± a singularity with infinite expansion (+) or contraction (-). The eigenvalues for the fixed points of the full system given in Table II.
At an Einstein point, the remaining expression left in the Hubble expansion variable (30), which we call f (x), is 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 → R we find f (x) → −2R, and taking the limit of f (x) when x → 1 we find f (x) → ∞. Therefore, between x = 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.
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 Setting this equal to zero and rearranging for c z , we find which we substitute into f (x) = 0 (64). Solving this equation for x, we find x 0.25 and x 0.60, which from (66) correspond to c z 0.20 and c z 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 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.

A. 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 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 turnaround. For all trajectories, the expansion (contraction) is always accelerating (a > 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.  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 xE 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.

B. 2 Einstein Points
The phase space for the first limiting case with two Einstein points where c z 0.20 is shown in Fig. 9. The Einstein point at x E 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 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 0.87, trajectories are always accelerating, therefore the bouncing models are always accelerating, and the re-collapsing models are always decelerating.

Extra bouncing region
The first case where three Einstein points exist is shown in Fig. 10. The Einstein point at x E 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 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 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 F.

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 0.17 and x E 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 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.

Extra turn-around region
The final case with three Einstein points is given in Fig. 12. The Einstein point at x E 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 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 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 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.

D. 2 Einstein Points
The phase space for the system for c z 0.38 is shown in Fig. 13, which is the second limiting case with two Einstein points. The Einstein point at x E 0.60 is a cusp (corresponding to the two de Sitter points coincid-ing), 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 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.

E. 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 = R at x E 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.

F. 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 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 = 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 in-  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 xE 0.14 (saddle) corresponds to the maximum of the potential. Bouncing models are obtained within the CFS for x < 0.14, which always accelerate.

terest.
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 turnaround 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 [2].  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.

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 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 (68) and r (69) to the continuity equation. We find 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 and x 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 and integrating r (69) with respect to x (29), we obtain . (74) We can then calculate the first integral c r in same way as before. Using the Planck 2018 density parameters [2], along with our definition of R (10) and the low energy cosmological constant ρ Λ (42), we find (75) for a value of R = 0.05.
The fixed points for this system are given in Table V. In this system, there are no physical singularities, but the fixed points dS 2± and dS 4± are coordinate singularities of the de Sitter spacetime when represented as an FLRW, see [9].
The eigenvalues of the fixed points for this system are shown in Table VI, and the linear stability classification for each point is given in Table VII. As previously, we cannot include the eigenvalues of the fixed points at Y = ±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 = R at Y = ±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 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 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  We need to solve numerically for the values x E , such that f (x E ) = 0. Taking the limit of f (x) when x → R, we find f (x) → −2R, and when x → 1, we find f (x) → −6. Therefore, this system does not necessarily admit any Einstein points between x = 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 .
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 0.25 and c z 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 → 0 as x → 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.  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 xE 0.23 (cusp) corresponding to the horizontal point of inflection in the potential. All closed models bounce, and these are always accelerating.

B. One Einstein point
The limiting case where c z 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 11.0. One Einstein point is admitted at x E 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 0.25, dark energy is dominant for most of the phase space, and when c z 11, dark matter is mostly dominant, except close to x = R where dark energy is dominant. All three components become comparable near x = 1. Here, all trajectories are always accelerating.

C. Two Einstein points
The case where two Einstein points are admitted is shown in Fig. 19. The Einstein point at x E 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 0.31 (centre), which corresponds to a local minimum of the potential. These bouncing models always accelerate, and the cyclic models evolve with an earlytime 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 earlyand 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.

D. 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 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 = ay and a /a = y + 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 noncompactified 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 → 1.

E. A note on R
As we mentioned previously, in order to alleviate the old cosmological constant problem we would set 10 −120 < R < 10 −60 , however purely for the purpose of illustrating the full range of dynamics in readable phase plots, we have set R = 0.05. If we had set 10 −120 < 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 < 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 = 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. for 0.25 < cz < 11.0 in terms of the scale factor a and noncompactified 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.
As a final remark, we note that in order for a bouncing model to be theoretically robust and realistic, the value of R would need to satisfy the bounds 10 −120 < 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 [41], 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 [42][43][44][45][46][47][48][49][50] (see also [41,51] 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 [40,41,51].
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 turnaround 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 [35][36][37]. 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 qual-itative 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 [40], 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.