{'fraction_non_alphanumeric': 0.04093567251461988, 'fraction_numerical': 0.010220747062852326, 'mean_word_length': 4.1471086705986036, 'pattern_counts': {'":': 0, '<': 193, '<?xml version=': 0, '>': 75, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Cosmo-dynamics and dark energy with non-linear equation of state: a quadratic model

We investigate the general relativistic dynamics of Robertson-Walker models with a non-linear equation of state (EoS), focusing on the quadratic case P = P_0 + \alpha \rho + \beta \rho^2. This may be taken to represent the Taylor expansion of any arbitrary barotropic EoS, P(\rho). With the right combination of P_0, \alpha and \beta, it serves as a simple phenomenological model for dark energy, or even unified dark matter. Indeed we show that this simple model for the EoS can produce a large variety of qualitatively different dynamical behaviors that we classify using dynamical systems theory. An almost universal feature is that accelerated expansion phases are mostly natural for these non-linear EoS's. These are often asymptotically de Sitter thanks to the appearance of an effective cosmological constant. Other interesting possibilities that arise from the quadratic EoS are closed models that can oscillate with no singularity, models that bounce between infinite contraction/expansion and models which evolve from a phantom phase, asymptotically approaching a de Sitter phase instead of evolving to a"Big Rip". In a second paper we investigate the effects of the quadratic EoS in inhomogeneous and anisotropic models, focusing in particular on singularities.


I. INTRODUCTION
Model building in cosmology requires two main ingredients: a theory of gravity and a description of the matter content of the universe. In general relativity (GR) the gravity sector of the theory is completely fixed, there are no free parameters. The matter sector is represented in the field equations by the energy-momentum tensor, and for a fluid the further specification of an equation of state (EoS) is required. Apart from scalar fields, typical cosmological fluids such as radiation or cold dark matter (CDM) are represented by a linear EoS, P = w ρ.
The combination of cosmic microwave background radiation (CMBR) [1,2], large scale structure (LSS) [3] and supernova type Ia (SNIa) [4] observations provides support for a flat universe presently dominated by a component, dubbed in general "dark energy", causing an accelerated expansion. The simplest form of dark energy is an ad hoc cosmological constant Λ term in the field equations, what Einstein called his "biggest blunder". However, although the standard ΛCDM "concordance" model provides a rather robust framework for the interpretation of present observations (see e.g. [2,5]), it requires a Λ term that is at odds by many order of magnitudes with theoretical predictions [6]. This has prompted theorists to explore possible dark energy sources for the acceleration that go beyond the standard but unsatisfactory Λ. With various motivations, many authors have attempted to describe dark energy as quintessence, kessence or a ghost field, i.e. with scalar fields with various properties. There have also been attempts to describe dark energy by a fluid with a specific non-linear EoS like the Chaplygin gas [7], generalized Chaplygin gas [8], van der Waals fluid [9], wet dark fluid [11] and other spe-cific gas EoS's [10]. Recently, various "phantom models" (w = P/ρ < −1) have also been considered [13,14]. More simply, but also with a higher degree of generality, many authors have focused on phenomenological models where dark energy is parameterized by assuming a w = P/ρ = w(a), where a = a(t) is the expansion scale factor (see e.g. [15,16]).
Another possibility is to advocate a modified theory of gravity. At high energies, modification of gravity beyond general relativity could come from extra dimensions, as required in string theory. In the brane world [17,18,19,20] scenario the extra dimensions produce a term quadratic in the energy density in the effective 4-dimensional energy-momentum tensor. Under the reasonable assumption of neglecting 5-dimensional Weyl tensor contributions on the brane, this quadratic term has the very interesting effect of suppressing anisotropy at early enough times. In the case of a Bianchi I braneworld cosmology containing a scalar field with a large kinetic term the initial expansion is quasi-isotropic [21]. Under the same assumptions, Bianchi I and Bianchi V brane-world cosmological models containing standard cosmological fluids with linear EoS also behave in a similar fashion 1 [23], and the same remains true for more general homogeneous models [24,25] and even some inhomogeneous exact solutions [26]. Finally, within the limitations of a perturbative treatment, the quadraticterm-dominated isotropic brane-world models have been shown to be local past attractors in the larger phase space of inhomogeneous and anisotropic models [27,28]. More precisely, again assuming that the 5-d Weyl tensor contribution to the brane can be neglected, perturbations of the isotropic models decay in the past. Thus in the brane scenario the observed high isotropy of the universe is the natural outcome of generic initial conditions, unlike in GR where in general cosmological models with a standard energy momentum tensor are highly anisotropic in the past (see e.g. [29]).
Recently it has been shown that loop quantum gravity corrections result in a modified Friedmann equation [30], with the modification appearing as a negative term which is quadratic in the energy density. Further motivation for considering a quadratic equation of state comes from recent studies of k-essence fields as unified dark matter (UDM) models 2 [31,32]. The general k-essence field can be described by a fluid with a closed-form barotropic equation of state. The UDM fluid discussed in [31] has a non-linear EoS of the form P ∝ ρ 2 at late times. More recently, it has been shown that any purely kinetic kessence field can be interpreted as an isentropic perfect fluid with an EoS of the form P = P (ρ) [33]. Also, low energy dynamics of the Higgs phase for gravity have been shown to be equivalent to the irrotational flow of a perfect fluid with equation of state P = ρ 2 [34].
Given the isotropizing effect that the quadratic energy density term has at early times in the brane scenario this then prompts the question: can a term quadratic in the energy density have the same effect in general relativity. This question is non-trivial as the form of the equations in the two cases is quite different. On the brane, for a given EoS the effective 4-dimensional Friedmann and Raychaudhuri equations are modified, while the continuity equation is identical to that of GR. With the introduction of a quadratic EoS in GR, the Friedman equation remains the same, while the continuity and Raychaudhuri equations are modified 3 .
Taking into account this question (to be explored in detail in Paper II [48]), the diverse motivations for a quadratic energy density term mentioned above and with the dark energy problem in mind, in this paper we explore the GR dynamics of homogeneous isotropic Robertson-Walker models with a quadratic EoS, P = P 0 + αρ + βρ 2 . This is the simplest model we can consider without making any more specific assumptions on the EoS [35]. It represents the first terms of the Taylor expansion of any EoS function P = P (ρ) about ρ = 0. It can also be taken to represent (after re-grouping of terms) the Taylor expansion about the present energy density ρ 0 , see [35]. In this sense therefore the out-coming dynamics is very general. Indeed it turns out that this simple model can produce a large variety of qualitatively different dynamical behaviors that we classify using dynamical systems theory [36,37]. An outcome of our analysis is that accelerated expansion phases are mostly natural for nonlinear EoS's. These are in general asymptotically de Sitter thanks to the appearance of an effective cosmological constant. This suggests that an EoS with the right combination of P 0 , α and β may provide a good and simple phenomenological model for UDM, or at least for a dark energy component. Other interesting possibilities that arise from the quadratic EoS are closed models that can oscillate with no singularity, models that bounce between infinite contraction/expansion and models which evolve from a phantom phase, asymptotically approaching a de Sitter phase instead of evolving to a "big rip" or other pathological future states [13,38,39].
As mentioned before, the question of the dynamical effects the quadratic energy density term has on the anisotropy in GR is explored in Paper II [48]. There we analyze Bianchi I and V models with the EoS P = αρ + βρ 2 , as well as perturbations of the isotropic past attractor of those models that are singular in the past. We anticipate that Bianchi I and V non-phantom models with β > 0 have an isotropic singularity, i.e. they are asymptotic in the past to a certain isotropic model, and that perturbations of this model decay in the past. Phantom anisotropic models with β > 0 are necessarily asymptotically de Sitter in the future, but the shear anisotropy dominates in the past. For β < 0 all models are anisotropic in the past, while their specific future evolution depends on the value of α.
The paper is organized as follows. In section II we outline the setup and the three main cases we will investigate. In section III, we study the dynamics of isotropic cosmological models in the high energy limit (neglecting the P 0 term). We find the critical points, their stability nature and the occurrence of bifurcations of the dynamical system. In section IV, we consider the low energy limit (neglecting the ρ 2 term). The full system is then analyzed in section V, showing the qualitatively different behavior with respect to the previous cases. We then finish with some concluding remarks and an outline of work in progress in section VI. Units are such that 8πG/c 4 = 1.

A. Dynamics with non-linear EoS
The evolution of Robertson-Walker isotropic models with no cosmological constant Λ term is given in GR by the following non-linear planar autonomous dynamical system:ρ where H is the Hubble expansion function, related to the scale factor a by H =ȧ/a. In order to close this system of equations, an EoS must be specified, relating the isotropic pressure P and the energy density ρ. When an EoS P = P (ρ) is given, the above system admits a first integral, the Friedman equation where K is the curvature, K = 0, ±1 as usual for flat, closed and open models.
Here we are interested in exploring the general dynamical features of a non-linear EoS P = P (ρ). Before considering the specific case of a quadratic EoS, we note some important general points.
First, it is immediately clear from Eq. (1) that an effective cosmological constant is achieved whenever there is an energy density value ρ Λ such that P (ρ Λ ) = −ρ Λ . More specifically: Remark 1. If for a given EoS function P = P (ρ) there exists a ρ Λ such that P (ρ Λ ) = −ρ Λ , then ρ Λ has the dynamical role of an effective cosmological constant.
Remark 2. A given EoS P (ρ) may admit more than one point ρ Λ . If these points exist, they are fixed points of Eq. (1).

Remark 4.
Remark 3 is only valid in GR, and a different condition will be valid in other theories of gravity. Remarks 1 and 2, however, are only based on conservation of energy, Eq. (1). The latter is also valid (locally) in inhomogeneous models, provided that the time derivative is taken to represent the derivative along the fluid flow lines (e.g. see [40]), and is a direct consequence of T ab ;b = 0. Thus Remarks 1 and 2 are valid in any gravity theory where T ab ;b = 0, as well as (locally) in inhomogeneous models.
Remark 6. Any point ρ Λ is a shunt 4 of the autonomous system Eq. (4) if either ρ + P (ρ) < 0 on both sides of ρ Λ , or ρ + P (ρ) > 0 on both sides of ρ > ρ Λ . In this case the fluid is respectively phantom or standard on both sides.
Let's now consider the specific case of a general quadratic EoS of the form: The parameter β sets the characteristic energy scale ρ c of the quadratic term as well as it's sign ǫ Remark 7. Eq. (5) represents the Taylor expansion, up to O(3), of any barotropic EoS function P = P (ρ) about ρ = 0. It also represents, after re-grouping of terms, the Taylor expansion about the present energy density value ρ 0 [35]. In this sense, the dynamical system (1,2) with (5) is general, i.e. it represents the late evolution, in GR, of any cosmological model with non-linear barotropic EoS approximated by Eq. (5).
The usual scenario for a cosmological fluid is a standard linear EoS (P 0 = β = 0), in which case α = w is usually restricted to the range −1 < α < 1. For the sake of generality, we will consider values of α outside this range, considering dynamics only restricted by the request that ρ ≥ 0. The first term in Eq (5) is a constant pressure term which in general becomes important in what we call the low energy regime. The second term is the standard linear term usually considered, with If it is positive, α has an interpretation in terms of the speed of sound of the fluid in the limit ρ → 0, α = c 2 s . The third term is quadratic in the energy density and will be important in what we call the high energy regime.
In the following, we first split the analysis of the dynamical system Eqs. (1,2,5) in two parts, the high energy regime where we neglect P 0 and the low energy regime where we set β = 0, then we consider the full system with EoS (5). Using only the energy conservation Eq. (1) we list the various sub-cases, also briefly anticipating the main dynamical features coming out of the analysis in Sections III, IV and V.

B. Quadratic EoS for the high energy regime
In the high energy regime we consider the restricted equation of state: The energy conservation Eq. (1) can be integrated in general to give: where ρ o , a o represent the energy density and scale factor at an arbitrary time t o . This is valid for all values of ǫ, ρ c and α, except for α = −1. In the case α = −1 the evolution of the energy density is: The EoS with this particular choice of parameters has already been considered as a possible dark energy model [38,44]. We will concentrate on the broader class of models where α = −1.
In Section III we will give a dynamical system analysis of the high energy regime, but it is first useful to gain some insight directly from Eq. (9).
We start by defining noticing that this is an effective positive cosmological constant point only if ǫ(1 + α) < 0. It is then convenient to rewrite Eq. (9) in three different ways, defining a ⋆ = |A| 1/3(α+1) , each representing two different subcases.
A1: ǫ > 0, (1 + α) > 0. In this case a ⋆ < a < ∞, with ∞ > ρ > 0. Further restrictions on the actual range of values that a and ρ can take may come from the geometry. For a subset of appropriate initial conditions closed (positively curved) models may expand to a maximum a (minimum ρ) and re-collapse, and for α < −1/3 not all closed models have a past singularity at a = a ⋆ , having instead a bounce at a minimum a (maximum ρ).
A2: ǫ < 0, (1 + α) < 0. In this case 0 < a < a ⋆ , with 0 < ρ < ∞, and the fluid exhibits phantom behavior. All models have a future singularity at a = a ⋆ , but in general closed models contract from a past singularity, bounce at a minimum a and ρ, then re-expand to the future singularity (we will refer to this as a phantom bounce).
B1: ǫ > 0, (1 + α) < 0, A > 0. In this case a ⋆ < a < ∞, with ∞ > ρ > ρ Λ . As in case A1, further restrictions on the actual range of values that a and ρ can take may come from the geometry. For a subset of initial conditions closed models may expand to a maximum a (minimum ρ) and re-collapse, while for another subset closed models don't have a past singularity at a = a ⋆ , having instead a bounce at a minimum a (maximum ρ).
B2: ǫ < 0, (1 + α) > 0, A < 0. In this case 0 < a < a ⋆ , with ρ Λ < ρ < ∞. As in the case A2, the fluid has a phantom behavior. All models have a future singularity at a = a ⋆ , with closed models contracting from a past singularity to a minimum a and ρ before re-expanding.

C. Low energy regime: affine EoS
In the low energy regime we consider the affine equation of state: This particular EoS has been investigated as a possible dark energy model [11,12], however, only spatially flat Friedmann models where considered. The scale factor dependence of the energy density is: This is valid for all values of P o and α except α = −1. In the case α = −1, the evolution of the energy density is: As in the high energy case, we will concentrate on the broader class of models where α = −1.
In Section IV we present the dynamical system analysis of the low energy regime, but first let us gain some insight from Eq. (17). As with the high energy case, in many cases the fluid violates the null energy condition (ρ+P < 0) and exhibit phantom behavior. Defining we see that a positive effective cosmological constant point exists,ρ Λ > 0, only if P o /(1 + α) < 0. Eq. (17) can be rewritten in three different ways, defining a ⋆ = |B| 1/3(α+1) , each representing two different subcases.
D1: P o > 0, (1 + α) > 0. In this case 0 < a < ∞, with ∞ > ρ > −|ρ Λ |. The geometry places further restrictions on the values that a and ρ can take. The subset of open models (negative curvature) are all non-physical as they evolve to the ρ < 0 region of the phase space. The spatially flat models expand to a maximum a (when ρ = 0) and recollapse. The closed (positively curved) models expand to a maximum a (minimum ρ) and recollapse, and for −1 ≤ α < −1/3 a subset of closed models oscillate between a maximum and minimum a (minimum and maximum ρ).
D2: P o < 0, (1 + α) < 0. In this case 0 < a < ∞, with −|ρ Λ | < ρ < ∞. In this case the fluid exhibits phantom behavior. The subset of open models are all non-physical as they evolve from the ρ < 0 region of the phase space. The spatially flat models contract, bounce at a minimum a when ρ = 0 and re-expand in the future. The closed models contract, bounce at a minimum a and ρ, then re-expand in the future.
E1: P o > 0, (1 + α) < 0, B > 0. In this case 0 < a < ∞, withρ Λ < ρ < ∞. As in the case D2, the fluid behaves in a phantom manner. The flat and open models are asymptotically de Sitter in the past, when their energy density approaches a finite value (ρ →ρ Λ as a → 0), and whenρ Λ becomes negligible in Eq. (22) they evolve as standard linear phantom models, reaching a future singularity in a finite time (ρ → ∞ as a → ∞). The closed models contract to a minimum a (minimum ρ), bounce and re-expand.
E2: P o < 0, (1 + α) > 0, B > 0. In this case 0 < a < ∞, with ∞ > ρ >ρ Λ . All flat and open models expand from a singularity and asymptotically evolve to a de Sitter model, with ρ =ρ Λ . The closed models evolve from a contracting de Sitter model to minimum a (maximum ρ), bounce and then evolve to an expanding de Sitter model.
The subset of open models are all non-physical as they evolve to the ρ < 0 region of the phase space. The flat models evolve from an expanding de Sitter phase to a contracting de Sitter phase. The closed models oscillate between a maximum and minimum a (minimum and maximum ρ).
F2: P o < 0, (1 + α) > 0, B < 0. In this case 0 < a < ∞, with −∞ < ρ <ρ Λ . The fluid exhibits phantom behavior. The open models are all non-physical as they evolve from the ρ < 0 region of the phase space. The flat and closed models evolve from a contracting de Sitter phase, bounce at minimum a and ρ, then re-expand, asymptotically approaching a expanding de Sitter phase.

D. The full quadratic EoS
In Section V we present the dynamical system analysis of the full quadratic EoS models given by Eq. (5), but again we first study the form of ρ(a) implied by conservation of energy, Eq. (1). As with the previous cases the fluid can violate the null energy condition (ρ+P < 0) and therefore may exhibit phantom behavior. The system may admit two (possibly negative) effective cosmological constant points: if is non negative. Clearly, the existence of the effective cosmological points depends on the values of the parameters in the EoS. This in turn affects the functional form of ρ(a). In order to find ρ(a) the following integral must be evaluated: This is done separately for the cases when no effective cosmological points exist (∆ < 0), when one cosmological point exist, ρ Λ,1 = ρ Λ,2 =ρ Λ = 0 (∆ = 0) and when two cosmological points exist, ρ Λ,1 = ρ Λ,2 = 0 (∆ > 0). We now consider these three separate sub-cases.
G1: β > 0, P o > 0. In this case a 1 < a < a 2 (where a 1 < a 2 ), with ∞ > ρ > −∞. The fluid behaves in a standard manner and all models have a past singularity at a = a 1 . All open models are non-physical as they evolve to the ρ < 0 region of the phase space. The flat models expand to a maximum a (ρ = 0) and then re-collapse. The closed models can behave in a similar manner to flat models except they reach a minimum ρ before re-collapsing. Some closed models oscillate between maxima and minima a and ρ.
G2: β < 0, P o < 0. In this case a 1 < a < a 2 (where a 1 < a 2 ), with −∞ < ρ < ∞. The fluid behaves in a phantom manner. All open models are non-physical as they evolve from the ρ < 0 region of the phase space. The flat and closed models represent phantom bounce models, that is they evolve from a singularity at a = a 1 (ρ = ∞), contract to a minimum a (minimum ρ) and then re-expand to the future singularity at a = a 2 .
The subset of open models are all non-physical as they evolve to the ρ < 0 region of the phase space. The flat models evolve from an expanding de Sitter phase to a contracting de Sitter phase. The closed models oscillate between maxima and minima a and ρ.
H2: β > 0, P o > 0, ρ >ρ Λ . In this case a 1 < a < ∞ with ∞ > ρ >ρ Λ and the fluid behaves in a standard manner. Ifρ Λ > 0, the open and flat models evolve from a past singularity (a = a 1 ) and evolve to a expanding de Sitter phase. For a subset of initial conditions closed models may expand to a maximum a (minimum ρ) and re-collapse, while for another subset closed models avoid a past singularity, instead having a bounce at a minimum a (maximum ρ). Ifρ Λ < 0, the open models are non-physical, while flat and closed models represent recollapse models.
The fluid behaves in a phantom manner. The open models are all non-physical as they evolve from the ρ < 0 region of the phase space. The flat and closed models evolve from a contracting de Sitter phase, bounce at minimum a and ρ, then re-expand, asymptotically approaching an expanding de Sitter phase.
H4: β < 0, P o < 0, ρ >ρ Λ . In this case 0 < a < a 1 withρ Λ < ρ < ∞ and the fluid behaves in a phantom manner. All models have a future singularity at a = a 1 . Ifρ Λ > 0, closed models contract from a past singularity to a minimum a and ρ before re-expanding (phantom bounce), while flat and open models are asymptotic to generalized de Sitter models in the past. Ifρ Λ < 0, open models are non-physical, while flat and closed models contract from a past singularity to a minimum a and ρ before re-expanding.
Note that β > 0 (< 0) implies ρ Λ,2 < ρ Λ,1 (ρ Λ,1 < ρ Λ,2 ), and In this case 0 < a < a 1 with ρ Λ,2 > ρ > −∞ and the fluid behaves in a standard manner. The open models are all non-physical as they evolve to the ρ < 0 region of the phase space. The flat models evolve from an expanding de Sitter phase to a contracting de Sitter phase. The closed model region contains a generalized Einstein static fixed point and models which oscillate indefinitely (between minima and maxima a and ρ).
I2: β > 0, P o > 0, ρ Λ,2 < ρ < ρ Λ,1 . In this case 0 < a < ∞ with ρ Λ,2 < ρ < ρ Λ,1 and the fluid behaves in a phantom manner. The open models evolve from one expanding de Sitter phase (ρ = ρ Λ,2 ) to more rapid (greater ρ and H) de Sitter phase (ρ = ρ Λ,1 ), however the spatial curvature is negative in the past and asymptotically approaches zero in the future. The flat models behave in a similar manner except that the curvature remains zero. The closed models undergo a phantom bounce with asymptotic de Sitter behavior, that is they evolve from a contracting de Sitter phase, reach a minimum a, minimum ρ and then evolve to a expanding de Sitter phase.
In this case a 1 < a < ∞ with ∞ > ρ > ρ Λ,1 and the fluid behaves in a standard manner. All flat and open models expand from a singularity at a = a 1 and asymptotically evolve to a expanding de Sitter phase (ρ = ρ Λ,1 ). A subset of closed models evolve from a contracting de Sitter phase to minimum a (maximum ρ), bounce and then evolve to an expanding de Sitter phase. Another subset of closed models expand from a singularity at a = a 1 , reach a maximum a and minimum ρ, only to re-collapse.
In this case a 1 < a < ∞, with −∞ < ρ < ρ Λ,1 and the fluid behaves in a phantom manner. The open models are all non-physical as they evolve from the ρ < 0 region of the phase space. The flat and closed models evolve from a contracting de Sitter phase, bounce at minimum a and ρ, then re-expand, asymptotically approaching a expanding de Sitter phase.
In this case 0 < a < ∞ with ρ Λ,2 > ρ > ρ Λ,1 and the fluid behaves in a standard manner. The open models evolve from a expanding de Sitter phase (ρ = ρ Λ,2 ) to less rapid (lower ρ and H) de Sitter phase (ρ = ρ Λ,1 ) with the spatial curvature being negative in the past and zero asymptotically in the future. The flat models behave in a similar manner, except that the curvature remains zero throughout the evolution. The closed models can undergo a phantom bounce with asymptotic de Sitter behavior in the future and past, a subset of these models enter a loitering phase both before and after the bounce. There are a subset of closed models which oscillate indefinitely.
In this case 0 < a < a 1 with ρ Λ,2 < ρ < ∞ and the fluid behaves in a phantom manner. All models have a future singularity at a = a 1 , with closed models contracting from a past singularity to a minimum a and ρ before re-expanding (phantom bounce).

E. The Singularities
In general, singularities may behave in qualitatively different ways. The singularities present for the non linear EoS are quite different from the standard "Big Bang"/"Big Crunch" singularity. The standard singularities are such that: • "Big Bang"/"Big Crunch" : For a → 0, ρ → ∞.
If the singularity occurs in the past (future) we refer to it as a "Big Bang" ("Big Crunch"). In order to differentiate between various types of singularities, we will use the following classification system for future singularities [38] (cf. also [45]): Here t ⋆ , a ⋆ , ρ ⋆ and |P ⋆ | are constants with a ⋆ = 0. The main difference in our case is that the various types of singularities may occur in the past or the future. The future singularity described in case A2 falls into the category of Type III, however, the past singularity mentioned in case A1 is also a Type III singularity. In the case of the full quadratic EoS, all singularities which occur for a finite scale factor (a = a 1 ) are of Type III.

III. HIGH ENERGY REGIME DYNAMICS
A. The dimensionless dynamical system It is convenient to describe the dynamics in terms of dimensionless variables. In the high energy regime these are: The system of equations (1)-(2) then changes into: and the Friedman equation (3) gives The discrete parameter ǫ denotes the sign of the quadratic term, ǫ ∈ {−1, 1}. The primes denote differentiation with respect to η, the normalized time variable.
The variable x is the normalized energy density and y the normalized Hubble function. We will only consider the region of the phase space for which the energy density remains positive (x ≥ 0). The system of equations above is of the form u ′ i = f i (u j ). Since this system is autonomous, trajectories in phase space connect the fixed/equilibrium points of the system (u j,o ), which satisfy the system of equations f i (u j,0 ) = 0. The fixed points of the high energy system and their existence conditions (the conditions for which x ≥ 0 and x, y ∈ R) are given in Table I.   TABLE I: Location and existence conditions (x ≥ 0 and x, y ∈ R) of the fixed points of the high energy regime system.
The first fixed point (M) represents an empty flat (Minkowski) model. The parabola y 2 = x/3 is the union of trajectories representing flat models, K = 0 in Eq. (34) (see Figs. 1 and 3-7). The trajectories below the parabola represent open models (K = −1), while trajectories above the parabola represent closed models (K = +1). The second fixed point (E) represents a generalized static Einstein universe. This requires some form of inflationary matter and therefore may only exists when The last two points represent expanding and contracting spatially flat de Sitter models (dS ± ). These points exist when the fluid permits an effective cosmological constant point, x Λ := ρ Λ /ρ c = −ǫ(α + 1); in addition x Λ > 0 must be true for the fixed points to be in the physical region of the phase space. There are further fixed points at infinity, these can be found by studying the corresponding compactified phase space. The first additional fixed point is at x = y = ∞ and represents a singularity with infinite expansion and infinite energy density. The second point is at x = ∞, y = −∞ and represents a singularity with infinite contraction and infinite energy density.

B. Generalities of stability analysis
The stability nature of the fixed points can be found by carrying out a linear stability analysis. In brief (see e.g. [37] for details), this involves analyzing the behavior of linear perturbations u j = u j,o + v j around the fixed points, which obey the equations v ′ i = Mv j . The matrix M is the Jacobian matrix of the dynamical system and is of the form: The eigenvalues λ i of the Jacobian matrix evaluated at the fixed points tell us the linear stability character of the fixed points. The fixed point is said to be hyperbolic if the real part of the eigenvalues is non-zero (R(λ i ) = 0). If all the real parts of the eigenvalues are positive (R(λ i ) > 0) the point is said to be a repeller. Any small deviations from this point will cause the system to move away from this state. If all the real parts are negative (R(λ i ) < 0), the point is said to be an attractor. This is because if the system is perturbed away from this state, it will rapidly return to the equilibrium state. If some of the values are positive, while others are negative then the point is said to be a saddle point. If the eigenvalues of the fixed point are purely imaginary then the point is a center. If the center nature of the fixed point is confirmed by some non-linear analysis, then the trajectories will form a set of concentric closed loops around the point. If the eigenvalues do not fall into these categories, we will resort to numerical methods to determine their stability.
The eigenvalues for the fixed points of the system (Eq.'s (33)) are given in Table II and the linear stability character is given in Table III.
We first consider the system when we have a positive quadratic energy density term (ǫ = +1) in the high energy regime EoS. We will concentrate on the region around the origin as this is where the finite energy density fixed points are all located. The plots have been created using the symbolic mathematics application Maple 9.5. The individual plots are made up by three layers, the first is a directional (represented by grey arrows) field plot of the state space. The second layer represents the separatrices and fixed points of the state space . A separatrix (black lines) is a union of trajectories that marks a boundary between subsets of trajectories with different properties and can not be crossed. The fixed points are represented by black dots. The final layer represents some example trajectories (grey lines) which have been calculated by numerically integrating the system of equations for a set of initial conditions. The character of the fixed point M is undefined and so is determined numerically. The fixed point representing the generalized Einstein static solution is a saddle point. The fixed points representing the generalized expanding (contracting) de Sitter points are attractor (repeller) points. The trajectories or fixed points in the y > 0 (y < 0) region represent expanding (contracting) models. We will mainly discuss the right hand side of the state space (expanding models) as in general the corresponding trajectory on the left hand side is identical under time reversal.
The phase space of the system is considered when α < −1 and is shown in Fig. 1. The lowest horizontal line (x = 0) is the separatrix for open models (K = −1) and will be referred to as the open Friedmann separatrix (OFS). The trajectories on the separatrix represent Milne models (x = 0, K = −1 and a(η) ∝ η) which are equivalent to a Minkowski space-time in a hyperbolic co-ordinate system. The second higher horizontal line (x Λ = −(α + 1)) is the separatrix which is the dividing line between regions of phantom (x < x Λ ) and non-phantom/standard behavior(x > x Λ ), we will call this the phantom separatrix (PS). The standard region corresponds to the case B1, while the phantom region corresponds to the case C1. In the phantom region the fluid violates the Null Energy Condition (ρ + P < 0). This means the energy density is increasing (decreasing) in the future for an expanding (contracting) universe. In the standard case of the linear EoS in GR, this occurs when w < −1 and ultimately leads to a Type I singularity [39,46]. The parabola (y 2 = x/3) represents the separatrix for flat Friedmann models (K = 0), we will call this the flat Friedmann separatrix (FFS). The inner most thick curve is the separatrix for closed Friedmann models (K = +1) and will be called the closed Friedmann separatrix (CFS). The separatrix has the form: The constant A is fixed by ensuring that the saddle fixed point coincides with the fixed point representing the generalized Einstein static model (E). The constant is given in terms of the EoS parameters and has the form: The Minkowski fixed point is located at the intersection of the OFS and FFS. The generalized flat de Sitter fixed points are located at the intersection of the PS and FFS. The generalized Einstein static fixed point is located on the CFS. The trajectories between the OFS and the PS (0 < x < x Λ ) represent models which exhibit phantom behavior (the case C1). The open models in the phantom region are asymptotic to a Milne model in the past and to a generalized flat de Sitter model (dS + ) in the future. The closed models in the phantom region evolve from a contracting de Sitter phase, through a phantom phase to an expanding de Sitter phase (phantom bounce). It is interesting to note that unlike the standard GR case the phantom behavior does not result in a Type I singularity but asymptotically evolves to a expanding de Sitter phase. This is similar to the behavior seen in the phantom generalized Chaplygin gas case [39]. such models to give: The region above the PS represents models which evolve in a non-phantom/standard manner (the case B1). The trajectories in the expanding region (y > 0) of the phase space are asymptotic to a Type III singularity in the past 5 . The trajectories outside the FFS represent open models which evolve from a Type III singularity to a flat de Sitter phase, as do the trajectories on the FFS. The trajectories in between the CFS and the FFS evolve from a Type III singularity to a flat expanding de Sitter phase but may enter a phase of loitering. Loitering is characterized by the Hubble parameter dipping to a low value over a narrow red-shift range, followed by a rise again. In order to see this more clearly, we have plotted the normalized Hubble parameter (y) as a function of scale factor for three different trajectories in Fig. 2. The top two curves represent the open and flat models, with 5 The Type III singularity appears to be a generic feature of the high energy regime EoS and can occur both in the future and the past. the Hubble parameter dropping off quicker for the flat Friedmann model. The lower most curve is the Hubble parameter for the closed model. The plot shows that the closed model evolves to a loitering phase. Loitering cosmological models in standard cosmology were first found for closed FLRW models with a cosmological constant. More recently, brane-world models which loiter have been found [47], these models are spatially flat but can behave dynamically like a standard FLRW closed model. The interesting point here is that the models mentioned above loiter without the need of a cosmological constant (due to the appearance of an effective cosmological constant), the topology is asymptotically closed in the past and flat in the future. The trajectories inside the CFS can have two distinct types of behavior corresponding to the central regions above and below the generalized Einstein static fixed point. Trajectories in the lower region represent closed models which evolve from a contracting de Sitter phase, bounce and then evolve to a expanding de Sitter phase. The trajectories in the upper region evolve from a Type III Singularity, expand to a maximum a (minimum x) and then re-collapse to a Type III singularity (we will refer to such re-collapsing models as turn-around models).
2. The −1 < α < −1/3 sub-case The phase space for the system when −1 < α < −1/3 is shown in Fig. 3. The fixed points representing the is no longer present and all trajectories represent models with non-phantom/standard fluids (this corresponds to the case A1). The main difference is that the generic future attractor is now the Minkowski model. The trajectories between the OFS and FFS now evolve from a Type III singularity to a Minkowski model, as do the flat Friedmann models. The models between the FFS and CFS now evolve from a Type III singularity to a Minkowski with the possibility of entering a loitering phase (as before the model is asymptotically flat in the future). The trajectories inside the CFS and above the Einstein static fixed point still represent turn-around models. The trajectories inside the OFS and below the Einstein static model now represent standard bounce models, that is they evolve from a Minkowski model, contract to a finite size, bounce and then expand to a Minkowski model.

The α ≥ −1/3 sub-case
Next we consider the system when α ≥ −1/3, the phase space is shown in Fig. 4. The fixed point representing the Einstein static models is now located in the x < 0 region of the phase space. The fluid in the entire physical region behaves in a non-phantom manner and corresponds to the case A1. The OFS and FFS remain the same and the CFS is no longer present. The trajectories between the OFS and FFS evolve from a Type III singularity to a Minkowski model. All trajectories above the FFS now represent turn-around models which start and terminate at a Type III singularity. The behavior of the models is qualitatively the same as that of the standard FLRW model with a linear EoS where w = α, in the linear EoS case the Type III singularity is replaced by a standard "Big Bang".

D. The ǫ = −1 case
We now consider the system when we have a negative quadratic energy density term (ǫ = −1) in the high energy regime EoS. The character of the fixed point M is still undefined. The fixed point representing the generalized Einstein static model is now a center. The fixed points representing the expanding/contracting flat de Sitter points now have saddle stability. As before, the trajectories or fixed points in the y > 0 (y < 0) region represent expanding models (contracting models), the black lines represent separatrix, grey lines represent example trajectories and fixed points are represented by black dots.
The phase space of the system when α < −1 is shown in Fig. 5. The horizontal line (x = 0) is still the open Friedman separatrix (OFS). The parabola is the flat Friedmann separatrix (FFS). The intersection of the OFS and FFS coincides with the Minkowski fixed point. All the trajectories in the physical region of the phase space exhibit phantom behavior (corresponding to the case A2), the energy density increases in an expanding model. The trajectories in the expanding (contracting) region in general evolve to a Type III singularity in the future (past). The trajectories between the OFS and the FFS are asymptotic to a Milne model in the past and are asymptotic to a Type III singularity in the future. The trajectories on the FFS start from a Minkowski model and enter a phase of super-inflationary expansion and evolve to a Type III singularity. Trajectories that start in a contracting phase during which the energy density decreases, reach a minimum a (minimum x) and then expand where the energy density increases represent phantom bounce models. The trajectories above the FFS represent closed models which evolve through a phantom bounce, but start and terminate in a Type III singularity. The behavior of the models is qualitatively the same as that of the FLRW models with a phantom linear EoS (w < −1) except that the Type III singularity is replaced by a Type I ("Big Rip") singularity.

The −1 < α < −1/3 sub-case
The phase space for the system when −1 < α < −1/3 is shown in Fig. 6. The lowest horizontal line (x = 0) is the OFS. The second higher horizontal line, x Λ = (α + 1) is the phantom separatrix (PS), this divides the state space into regions of phantom (x > x Λ ) and nonphantom/standard behavior (x < x Λ ). The phantom region corresponds to the case B2 and the standard region corresponds to the case C2. The flat de Sitter (dS ± ) points are located at the intersection of the FFS and the PS. The open models in the standard matter region (0 < x < x Λ ) are past asymptotic to open expanding de Sitter models in the past and evolve to Minkowski models in the future. The closed models in the region represent the standard bounce models, that is they evolve from a Minkowski model, contract to a minimum a (maximum x) and then expand to a Minkowski model. The trajectories above the PS (x > (α + 1)) all exhibit phantom behavior. The open models in this region are past asymptotic to open de Sitter models in the past and evolve to a Type III singularity in the future. The closed models in the region all represent models which undergo a phantom bounce but start and terminate in a Type III singularity.  We now consider the system when α ≥ −1/3, the phase space is shown in Fig. 7. The OFS, FFS and the PS are all still present, the phantom regions still corresponds to the case B2 and the standard region to the case C2. The trajectories in the phantom region (x > x Λ ) behave in a similar manner to the previous case, as do the open models in the standard matter region (0 < x < x Λ ). The main difference is in the region representing closed models (K = 1) with non-phantom behavior. There is now a new fixed point which represents a generalized Einstein static model (E). The closed models in the region now represent oscillating models. This is represented by closed concentric loops centered on the Einstein static fixed point. These oscillating models also appear in the low energy system and will be discussed in more detail later.

A. The dimensionless dynamical system
We now consider the system of equations for the low energy regime EoS, which can be simplified and expressed in terms of the following dimensionless variables: The system of equations is then: x ′ = −3y (ǫ p + (α + 1)x) .
The discrete parameter ǫ p denotes the sign of the pressure term, ǫ p ∈ {−1, 1}. The primes denote differentiation with respect to the new η. The variables x and y are the new normalized energy density and Hubble parameter. As before only the positive energy density region of the phase space will be considered. The fixed points of the system and the existence conditions are given in Table IV. As before, by existence we mean the conditions on the parameters to insure x ≥ 0 and x, y ∈ R.
The Minkowski model (x = y = 0) is no longer a fixed point of the system. The first fixed point (E) represents a generalized static Einstein model. This requires the overall effective equation of state to be that of inflationary matter and therefore only exists when ǫ p /(3α + 1) < 0. The last two points represent generalized expanding and contracting flat de Sitter models. These points only exist if the fluid permits an effective cosmological constant , alsox Λ ≥ 0 for the points to be in the physical region of the phase space. The eigenvalues of the equilibrium points are given in Table V, while the linear stability character is given in Table VI.
We start by considering the system when we have a positive constant pressure term (ǫ p = +1) in the low energy regime EoS. The Minkowski (x = y = 0) point is no longer present and the Einstein static solution has the stability character of a center. The fixed points representing the generalized expanding/contracting de Sitter points (dS ± ) now have saddle stability. As before black lines represent separatrix, grey lines represent example trajectories and black dots represent fixed points of the system. The phase space for the system when α < −1 is shown in Fig. 8. The open Friedmann separatrix (x = 0) is no longer present, and the x = y = 0 point is no longer a fixed point of the system. The horizontal line (x Λ = −(α+1) −1 ) is the phantom separatrix (PS), dividing the state space into regions with phantom (x >x Λ ) and standard (x <x Λ ) behavior. The phantom region corresponds to the case E1 and the standard region corresponds to the case F1. The parabola y 2 = x/3 is the separatrix representing the flat Friedmann models (FFS), this divides the remaining trajectories into open and closed models. The intersection of the PS and FFS coincides with the fixed points of the generalized flat de Sitter models.
The trajectories in the upper region that start in a contracting phase (during which the energy density decreases), reach a minimum a (minimum x) and then expand, representing phantom bounce models which terminate in a Type I singularity. The closed models in the phantom region (x >x Λ ) represent phantom bounce models which start and terminate in a Type I singularity 6 . The open models in the phantom region are asymptotic to open de Sitter models in the past and evolve to 6 The Type I singularity is a generic feature of the low energy regime EoS and can appear both in the future and the past.
a Type I singularity in the future. The trajectories below the PS (x <x Λ ) represent models which all behave in a standard manner (the F1 case). The open models in this region are all non-physical as they all evolve to the x < 0 region of the phase space. The region corresponding to closed models (above the FFS) contain a fixed point which represents the generalized Einstein static (E) model. The region is filled by a infinite set of concentric closed loops centered on the Einstein static fixed point, the closed loops represent oscillating models. The Friedmann equation for such models is given by: The constant B is fixed by the location of the Einstein fixed point (E). The constant is given in terms of α and ǫ p : These oscillating models appear for −∞ < α < −1/3 when ǫ p = +1 and are qualitatively similar to the oscillating models seen in the high energy case. The exact behavior of the variables for these models can be calculated by fixing the EoS parameter α. The qualitative behavior remains the same for the models, however the maximum and minimum values of the variables change.
In the case of α = −2/3 the equations are greatly simplified, the scale factor oscillates such that: The maximum and minimum scale factor is then: The normalized hubble parameter (y) is: The maximum and minimum y is given by: The normalized energy density (x) is given by: The maximum and minimum x are:

The −1 < α < −1/3 sub-case
We now consider the case when −1 < α < −1/3, the phase space is shown in Fig. 9. All trajectories in the physical region of the phase space exhibit standard behavior and correspond to the case D1. There is only one fixed point in the x ≥ 0 region of the phase space, this point represents the generalized Einstein static model (E). The FFS represent models which evolve from a standard "Big Bang", evolve to a Minkowski model and then to a standard "Big Crunch" (turn around model). The open models (below the FFS) are non-physical as the evolve into the x < 0 region. The trajectories above the separatrix represent closed models (K > 0) which oscillate indefinitely between a maximum and minimum a (minimum and maximum x), as seen in the previous case.

The −1/3 ≤ α sub-case
The phase space for the system when −1/3 ≤ α is shown in Fig. 10. As in the previous sub-case the fluid behaves in a standard manner and corresponds to the case D1. There are no fixed points in the physical region of the phase space. The parabola is the FFS and represents flat model which evolve from a "Big Bang", approach a Minkowski model and then re-collapse (turnaround models) to a "Big Crunch". The open models (below the FFS) are all non-physical as they evolve to the negative energy density region (x < 0) of the phase space. The closed models evolve from a "Big Bang", reach a maximum a (minimum x) and re-collapse to a "Big Crunch".

C. The ǫp = −1 case
We now consider the system when we have a negative constant pressure term (ǫ p = −1) in the low energy regime EoS. As before,the Minkowski (x = y = 0) point is no longer a fixed point of the system and the OFS is not present. The fixed point representing the generalized Einstein static model (E) has the stability character of a saddle. The fixed points representing the generalized expanding (contracting) flat de Sitter points now have attractor (repeller) stability.

The α < −1 sub-case
The phase space for the low energy system when α < −1 is shown in Fig. 11. All the trajectories in the x > 0 region of the phase space now exhibit phantom behavior and correspond to the case D2. The open models are all non-physical as they all evolve from the negative energy density region of the phase space. The flat and closed models represent phantom bounce models which start and end in a Type I singularity. They evolve from a Type I singularity, contract, bounce at a minimum a (minimum x) and expand to a Type I singularity.

The −1 < α < −1/3 sub-case
The phase space for the system when −1 < α < −1/3 is shown in Fig. 12. The horizontal line,x Λ = (α + 1) −1 is the phantom separatrix (PS), dividing the state space into regions with phantom (x <x Λ ) and standard behavior (x >x Λ ). The standard region corresponds to the case E2 and the phantom region corresponds to the case F2. The intersection of the PS and FFS coincides with the fixed points of the generalized flat de Sitter models (dS ± ). The flat expanding (contracting) de Sitter model is the generic future attractor (repeller). The open models in the standard matter region (x >x Λ ) evolve from a standard "Big Bang" to a flat expanding de Sitter phase. The closed models in this region evolve from a contracting flat de Sitter phase, reach a minimum a (maximum x), bounce and then evolve to expanding flat de Sitter phase. These models represent standard bounce models with asymptotic de Sitter behavior. The open models in the phantom region (x <x Λ ) are all non-physical. The flat and closed models in this region represent models exhibiting phantom bounce behavior which avoid the "Big Rip" and instead evolve to an expanding flat de Sitter phase.  We now consider the system in the parameter range −1/3 ≤ α, the phase space is shown in Fig. 13. The PS (x Λ = (α + 1) −1 ), FFS (y 2 = x/3) and generalized flat de Sitter points (dS ± ) still remain. The flat expanding (contracting) de Sitter model is the generic future attractor (repeller). The inner most black curve is the closed Friedmann separatrix (CFS) and coincides with the generalized Einstein static fixed point (E), which has saddle stability. The CFS is given by: . (51) The constant D is a constant of integration and can be fixed by the location of the fixed point E. The constant is given in terms of α and has the form: The region below the PS (x <x Λ ) remains qualitatively the same. The open models in the standard matter region (x >x Λ ) all evolve from a "Big Bang" to a expanding flat de Sitter phase. The trajectories between the FFS and the CFS also evolve from a "Big Bang" to a generalized expanding flat de Sitter model with the possibility of entering a loitering phase. The models inside the CFS can behave in one of two ways. The trajectories above the generalized Einstein static point represent turn-around models which evolve from a "Big Bang", reach a maximum a (minimum x) and then re-collapse to a "Big Crunch". The trajectories below evolve from a contracting de Sitter phase to an expanding de Sitter phase and represent bounce models.

V. THE FULL SYSTEM
A. The dimensionless dynamical system We now consider the system of equations with the full quadratic EoS, this can be simplified in a similar fashion to the previous case by introducing a new set variables: The system of equations then become: x ′ = −3y ν + (α + 1)x + ǫx 2 .
The parameter ǫ denotes the sign of the quadratic term, ǫ ∈ {−1, 1}. The parameter ν is the normalized constant pressure term. The primes denote differentiation with respect to the new normalized time variable η and only the physical region of the phase space is considered (x ≥ 0). The fixed points and their existence conditions are given in Table VII. The phase space undergo's a topological change for special values of the ν parameter, these values can be expressed in terms of α and are: As in the previous case, by existence we mean x ≥ 0 and x, y ∈ R. The general eigenvalues derived from the linear stability analysis are given in Table IX. The linear stability character of the fixed points is given in Table VIII.
The system has six fixed points and the sign of ǫ no longer effects the linear stability character of the fixed point (but changes it's position in the x − y plane). The first fixed point M represents a Minkowski model and is only present if ν = 0, the linear stability character is undefined. The second fixed point E 1 represents an Einstein static model and has the linear stability character of a saddle. The third fixed point E 2 represents an Einstein static model with the linear stability character of a center. In general, this fixed point is surrounded by a set of closed concentric loops representing oscillating models. The next pair of fixed points dS 1,± represents a set of generalized flat de Sitter models, the expanding (contracting) model has attractor (repeller) stability. The next pair of fixed points dS 2,± also represents a set of generalized flat de Sitter models, but now they have saddle stability. The separatrix for open Friedmann models (OFS) is only present if ν = 0. The parabola y 2 = x/3 (FFS) still separates the open and closed models. The separatrix for the closed Friedmann models (CFS) is present for a narrow range of the parameters and always coincides with the fixed points representing the generalized Einstein static model, E 1 . The fluid permits two possible effective cosmological constants points, they are given by: Where δ = (α + 1) 2 − 4ǫν. There is also a separatrix associated with each of the effective cosmological points, which divide the regions of phantom and non-phantom behavior. The separatrices will be referred to as the phantom separatrix (P S i which corresponds to the line x = x Λ,i ), with the appropriate subscript. For special choices of parameters the separatrices coincide. The discussion of the system will be split into the two categories, ǫ = +1 and ǫ = −1. The locations of the fixed points of the full system. The existence conditions are also given, that is x, y, ∈ R and x ≥ 0. To simplify the expressions special values of ν are used which can be expressed in terms of α, these values are ν1 = (3α+1) 2 36 and ν2 = (α+1) 2

B. The ǫ = +1 case
We first consider the system when we have a positive quadratic energy density term (ǫ = +1). The dynamical system can be further sub-divided into sub-cases with different values of the parameters α and ν. The various subcases have been highlighted in Table X. The majority of sub-cases result in a phase space diagram which is qualitatively similar to cases discussed in previous sections. That is the qualitative behavior of trajectories is the same even though the functional form of ρ(a) is different. The figure numbers not in bold (standard text) indicate choices of variable for which the phase space is qualitatively similar to a previous case, with the following differences: • The Type I singularities are now replaced by Type III singularities.
There is a narrow range of the parameters for which the state space is qualitatively different. The figure numbers given in bold in Table X, indicate the choice of variables for which the state space is qualitatively different to previously discussed cases. We will now discuss the four sub-cases which are different to those discussed in previous sections.
The phase space of the system when α < −1 and ν = ν 1 is shown in Fig. 14. As before the black lines represent separatrix, grey lines represent example trajectories and fixed points are represented by black dots. The fluid behaves in a standard manner and corresponds to the case G1. This choice of parameters results in the two Einstein points (E i ) coinciding. The resulting fixed point is highly non-linear and cannot be classified into the standard linear stability categories as in previous cases. The fixed point coincides with the CFS and the parabola is the FFS. The open models are all non-physical as they IX: The Eigenvalues derived from the linear stability analysis of the fixed points for the full system. In order to simplify the form of the eigenvalues we introduce the following variables, γ1 = (3α + 1), γ2 = (α + 1) and δ = (α + 1) 2 − 4ǫν. These eigenvalues are only valid for the choice of parameters consistent with the existence conditions. Name λ1 λ2  evolve to the x < 0 region of the phase space. The models between the FFS and the CFS represent turn-around models which evolve from a Type III singularity 7 , evolve to a maximum a (minimum x) and then re-collapse. The trajectories above the CFS also represent similar turnaround models.
2. The α < −1, ν2 < ν < ν1 sub-case The phase space of the system when α < −1 and ν 2 < ν < ν 1 is shown in Fig. 15. The fluid behaves in a standard manner and corresponds to the case G1. The Einstein fixed point of the previous case splits into two individual Einstein fixed points (E i ) via bifurcation. The first Einstein fixed point (E 1 ) coincides with the 7 As with the case of the high energy EoS, the Type III singularity is a generic feature of the fully quadratic EoS. CFS, while the second Einstein fixed point (E 2 ) is located inside the lower region enclosed by the CFS. Only the trajectories above the FFS differ from the previous case. The trajectories between the CFS and FFS still represent turn-around models which evolve from a Type III singularity but may now enter a loitering phase. The trajectories inside the CFS and above the Einstein fixed point (E 1 ), evolve from a Type III singularity, reach a maximum a and then re-collapse to a Type III singular- The phase space for the full system with ǫ = +1, α < −1 and ν2 < ν < ν1 (additionally when −1 < α < −1/3 and 0 < ν < ν1). The entire region corresponds to the case G1.
ity. The trajectories below E 1 represent closed oscillating models, the closed loops are centered on the second Einstein fixed point (E 2 ).
3. The α < −1, ν = ν2 sub-case The phase space of the system when α < −1 and ν = ν 2 is shown in Fig. 16. The fluid behaves in a standard manner in both regions and the upper (lower) region corresponds to the case H2 (H1). This choice of parameters results in the two sets of generalized de Sitter points (dS i,± ) coinciding. The fixed points coincide with P S i (x = x Λ,i ) and the FFS. The resulting fixed points are highly non-linear, the points have shunt stability along the FFS direction and the generalized expanding (contracting) de Sitter point has attractor (repeller) stability along the P S i direction. The two Einstein points (E i ) and the CFS are still present. In the x < x Λ,i region, the open models are all non-physical as they evolve to the x < 0 region and the closed models represent oscillating models which are centered on the Einstein point (E 2 ) with center linear stability. In the x > x Λ,i region, the open models are asymptotic to a Type III singularity in the past and a expanding flat de Sitter phase (dS i,+ ) in the future. The trajectories between the FFS and the CFS evolve from a Type III singularity to dS i,+ with the possibility of entering a loitering phase. The models inside the CFS and above the E 1 point represent turn- around models which asymptotically approach a Type III singularity. The closed models inside the CFS and below the E 2 point are asymptotic to a contracting de Sitter model phase (dS i,− ) in the past and a expanding de Sitter phase (dS i,+ ) in the future. The generic attractor in the x > x Λ,i region is the dS i,+ fixed point. 4. The α < −1, 0 < ν < ν2 sub-case The phase space of the system when α < −1 and 0 < ν < ν 2 is shown in Fig. 17. The upper (lower) horizontal line is the P S 1 (P S 2 ). The region above P S 1 corresponds to the case I3 and is qualitatively similar to the H2 region in the previous sub-case. The region below P S 2 corresponds to the case I1 and is qualitatively similar to the H1 region in the previous sub-case. The set of generalized flat de Sitter fixed points (dS i,± ) of the previous case split into two sets of generalized flat de Sitter fixed points via bifurcation. The upper (lower) set of generalized de Sitter points, dS 1,± (dS 2,± ) have attractor/repeller (saddle) stability. The region between P S 1 and P S 2 corresponds to the case I2 and the fluid behaves in a phantom manner. The open models in this region are asymptotic to open de Sitter models in the past and flat de Sitter models in the future. The closed models in the phantom region represent phantom bounce models which asymptotically approach a expanding (contracting) de Sitter phases in the future (past). We now consider the system when we have a negative quadratic energy density term (ǫ = −1). As before the system can be sub-divided into various sub-cases with different values of parameters of α and ν. The various sub-cases have been highlighted in Table XI   TABLE XI: The various sub-cases of the ǫ = −1 full system. The figure numbers given in bold, indicate the choice of variables for which the phase space is qualitatively different to previous cases. As before, the figure numbers not in bold (standard text) indicate choices of variable for which the phase space is qualitatively similar to a previous case, with the following differences: • The Type I singularities are now replaced by Type III singularities.
There are choices of parameters for which the phase space is different (figure numbers in bold in Table XI) and these four sub-cases will be discussed in the following sections.
The phase space of the system when α > −1/3 and −ν 1 < ν < 0 is shown in Fig. 18. The upper (lower) horizontal line at x = x Λ,2 (x = x Λ,1 ) is the P S 2 (P S 1 ) (they have swapped position with respect to the ǫ = +1 case). The region above P S 2 corresponds to the case I6, the region below P S 1 corresponds to the case I4 and the fluid behaves in a phantom manner in both regions. The region between P S 1 and P S 2 corresponds to the case I5 and the fluid behaves in a standard manner. The lower set of generalized de Sitter points (dS 1,± -at the intersection of P S 1 and FFS) have attractor/repeller stability, while the upper set (dS 2,± -at the intersection of P S 2 FIG. 19: The phase space for the full system with ǫ = −1, α > −1/3 and ν = −ν1. The upper, middle and lower regions correspond to the cases I6, I5 and I4 respectively. and FFS) have saddle stability. The CFS is located in between P S 1 and P S 2 and coincides with the Einstein point (E 1 ). The open models in the x < x Λ,1 region (the case I4) are all non-physical as they evolve from the x < 0 region of the phase space. The closed models in this region represent phantom bounce models which evolve from a contracting de Sitter phase (dS 1,− ) to a expanding de Sitter phase (dS 1,+ ). The open models in the standard region (x Λ,1 < x < x Λ,2 corresponding to the case I5 ) are asymptotic to a generalized open de Sitter model in the past and generalized flat de Sitter model in the future (the future attractor has lower x and y). The models between the CFS and the FFS in this region represent bounce models which evolve from a contracting de Sitter phase to a expanding de Sitter phase with the possibility of entering a loitering phase. The models enclosed by the CFS can be split into two groups. The models above the fixed point, E 1 represent oscillating models, the closed loops are centered on the fixed point E 2 . The models below the fixed point, E 1 represent bounce models which evolve from dS 1,− to dS 1,+ . In the x > x Λ,2 region (the case I6) the open models are asymptotic to generalized open de Sitter models in the past and a Type III singularity in the future. The closed models in this region represent phantom bounce models which evolve from a Type III singularity, reach a minimum a (minimum x) and then evolve to a Type III singularity. The generalized expanding flat de Sitter model, dS 1,+ (Type III singularity) is the generic future attractor in the region x < x Λ,2 (x > x Λ,2 ). The trajectories in the regions, FIG. 20: The phase space for the full system with ǫ = −1, α > −1/3 and −ν2 < ν < −ν1 (additionally when −1 < α < −1/3 and −ν2 < ν < 0). The upper, middle and lower regions correspond to the cases I6, I5 and I4 respectively.
x < x Λ,1 and x > x Λ,2 remain qualitatively similar in the following two cases (Fig.19, 20). The phase space of the system when α > −1/3 and ν = −ν 1 is shown in Fig. 19. The phase space is equivalent to the previous sub-case, except for the region x Λ,1 < x < x Λ,2 (the case I5). The open models in this region are still asymptotic to generalized open (flat) de Sitter models in the past (future). The behavior of the closed models has now changed, there are no longer trajectories representing oscillating models. The two generalized Einstein fixed points (E i ) have now coalesced to form one fixed point via bifurcation. The closed models above E i represent bounce models which evolve from dS 1,− to dS 1,+ , with the possibility of entering a loitering phase. The closed models below E i represent bounce models which evolve from dS 1,− to dS 1,+ without entering a loitering phase. any fixed points representing generalized Einstein static models and the CFS is no longer present. The open models in the region behave as in previous sub-cases. The closed models in the region represent bounce model, which evolve to a expanding (collapsing) de Sitter phase in the future (past) without the possibility of entering a loitering phase. 4. The α > −1/3, ν = −ν2 sub-case The next case is the phase space of the system when α > −1/3 and ν = −ν 2 and is shown in Fig. 21. The fluid behaves in a phantom manner in both regions and the upper (lower) region corresponds to the case H4 (H3). The two sets of generalized de Sitter points (dS i,± ) have now coalesced into a single set of generalized de Sitter points (dS ± ) which are located at the intersection of the FFS and the P S i which have also coalesced to form a single separatrix (x Λ,1 = x Λ,2 ). The resulting fixed points are highly non-linear, the points have shunt stability along the FFS direction and the generalized expanding (contracting) de Sitter point has attractor (repeller) stability along the P S i direction. The Type III singularity is the generic attractor in the upper region (x > x Λ,i ) and the dS i,+ is the generic attractor in the lower region (x < x Λ,i ).

VI. DISCUSSION AND CONCLUSIONS
In this paper we have systematically studied the dynamics of homogeneous and isotropic cosmological models containing a fluid with a quadratic EoS. This has it's own specific interest (see Section I for a variety of motivations) and serves as a simple example of more general EoS's. It can also be taken to represent the truncated Taylor expansion of any barotropic EoS, and as such it serves (with the right choice of parameters) as a useful phenomenological model for dark energy, or even UDM. Indeed, we have shown the dynamics to be very different and much richer than the standard linear EoS case, finding that an almost generic feature of the evolution is the existence of an accelerated phase, most often asymptotically de Sitter, thanks to the appearance of an effective cosmological constant. Of course to properly build physical cosmological models would require to consider the quadratic EoS for dark energy or UDM together with standard matter and radiation. Our analysis was aimed instead to derive and classify the large variety of different dynamical effects that the quadratic EoS fluid has when is the dominant component. In this respect, it should be noticed that a positive quadratic term in the EoS allows, in presence of another fluid such as radiation, equidensity between the two fluid to occur twice, i.e. the quadratic EoS fluid can be dominant at early and late times, and subdominant in an intermediate era.
In Section II we have made some general remarks, mostly based on conservation of energy only and as such valid independently of any specific theory of gravity. We have also given the various possible functional forms of the energy density as a function of the scale factor, ρ(a), and listed the many subcases, grouped in three main cases, what we call: i) the high energy models (no constant P o term); ii) the low energy affine EoS with no quadratic term; iii) the complete quadratic EoS.
The quadratic term in the EoS affects the high energy behavior as expected but can additionally affect the dynamics at relatively low energies. First, in Section III, we have concentrated on the high energy models. The specific choice of parameters fixes the behavior of the fluid, it can behave in a phantom or standard manner. In the case of phantom behavior, ρ can tend to zero at early times and either tend to an effective cosmological constant (C1) or a Type III singularity (A2) at late times. Alternatively ρ can also tend to an effective cosmological constant in the past (B2) and a Type III singularity at late times. When the fluid behaves in a standard manner, it can tend to a Type III singularity at early times, with ρ either tending to zero (A1) or to an effective cosmological constant (B1) at late times. The fluid can also behave as an effective cosmological constant at early times with ρ decaying away at late times (C2). The effective cosmological constant allows for the existence of generalized Einstein static (E) and flat de Sitter fixed (dS ± ) points which modify the late time behavior. The main new feature is the existence of models which evolve from a Type III singularity and asymptotically approach a flat de Sitter model (dS + ). Of specific interest are the closed models of this type, which can also evolve through an intermediate loitering phase.
Neglecting the quadratic term, in Section IV we have considered the low energy models with affine EoS. As expected, the constant term in the quadratic EoS affects the relatively low energy behavior. It can result in a variety of qualitatively different dynamics with respect to those of the linear EoS case. Again, the fluid can have a phantom or standard behavior. When the fluid behaves in a phantom manner, ρ can tend to an effective cosmological constant (F2), or can tend to a Type I ("Big Rip") singularity (D2) at late times. Alternatively, ρ can also tend to an effective cosmological constant in the past and a Big Rip in the future(E1). When the fluid behaves in a standard manner, we recover the linear EoS at early times and ρ can either tend to zero (D1) or to an effective cosmological constant (E2) at late times. The fluid can also behave as an effective cosmological constant at early times, with ρ decaying away at late times (F1). The effective cosmological constant allows for the existence of new fixed points(E and dS ± ). Comparing with standard linear EoS cosmology, the most interesting differences are new closed models which oscillate indefinitely and new closed models which exhibit phantom behavior which do not terminate in a "Big Rip", but asymptotically approach an expanding flat de Sitter model (flat and closed models where the fluid behaves as case F2).
When we study the dynamics of the system with the complete quadratic EoS, Section V, we see the appearance of new fixed points representing generalized Einstein and de Sitter models which are not present in the high/low energy systems. The various models of the simplified systems are present in the full system (but with differing ρ(a)), but there are also models with qualitatively new behavior. As with the previous cases, in the case of phantom behavior, ρ can tend to zero at early times and either tend to an effective cosmological constant (H3 and I4) or a Type III singularity (G2) at late times. Alternatively ρ can also tend to an effective cosmological constant in the past (H4 and I6) and a Type III singularity at late times. Finally, in the phantom case ρ can also tend to an effective cosmological constant both in the past and future (I2). In the case of standard behavior the fluid can tend to a Type III singularity at early times, with ρ either tending to zero (G1) or to an effective cosmological constant (H2 and I3) at late times. The fluid can also behave as an effective cosmological constant at early times with ρ decaying away at late times (H1 and I1). Finally, in the standard fluid case ρ can also tend to an effective cosmological constant both in the past and future (I5). There are models which evolve from a Type III singularity, reach a maximum a (minimum x) and then evolve to Type III singularity. These also enter a loitering phase before and after the turn around point. We also see bounce models which enter a loitering phase and asymptotically tend to generalized expanding (contracting) de Sitter models at late (early) times.
Of specific interest are models which evolve from a Type III singularity as opposed to the standard "Big Bang" (A1, B1). The simplest models of this type correspond to the high energy EoS with a positive quadratic term (is possible to recover standard behavior at late times). For these models the positive quadratic energy density term has the potential to force the initial singularity to be isotropic. The effects of such a fluid on anisotropic Bianchi I and V models is investigated in Paper II [48]. This is achieved by carrying out a dynamical systems analysis of these models. Additionally, using a linearized perturbative treatment we study the behavior of inhomogeneous and anisotropic perturbations at the singularity. The singularity is itself represented by an isotropic model and, If the perturbations of the latter decay in the past, this model represents the local past attractor in the larger phase space of inhomogeneous and anisotropic models (within the validity of the perturbative treatment). This would mean that in inhomogeneous anisotropic models with a positive non-linear term (at least quadratic) in the EoS isotropy is a natural outcome of generic initial conditions, unlike in the standard linear EoS case where generic cosmological models are, in GR, highly anisotropic in the past.