Cosmological evolution with quadratic gravity and nonideal fluids

Some cosmological models based on the gravitational theory $f(R) = R+\zeta R^2$, and on fluids obeying to the equations of state of Redlich-Kwong, Berthelot, and Dieterici are proposed for describing smooth transitions between different cosmic epochs. A dynamical system analysis reveals that these models contain fixed points which correspond to an inflationary, a radiation dominated and a late-time accelerating epoch, and a nonsingular bouncing solution, the latter being an asymptotic fixed point of the compactified phase space. The infinity of the compactified phase space is interpreted as a region in which the non-ideal behaviors of the previously mentioned cosmic fluids are suppressed. Physical constraints on the adopted dimensionless variables are derived by demanding the theory to be free from ghost and tachyonic instabilities, and a novel cosmological interpretation of such variables is proposed through a cosmographic analysis. The different effects of the equation of state parameters on the number of equilibrium solutions and on their stability nature are clarified. Some generic properties of these models, which are not sensitive to the particular fluid considered, are identified, while differences are critically examined by showing that the Redlich-Kwong scenario admits a second radiation-dominated epoch and a Big Rip Singularity.


I. INTRODUCTION
Despite being a challenging task, the search for a unified cosmological theory accounting for the entire known evolution history of the Universe, or at least providing a smooth transition between two different cosmic epochs, has been attempted both through single fluid approaches and by proposing modifications of the gravity sector beyond general relativity [1,2]. In the former case a certain single cosmic fluid is adopted to describe two different epochs in the limits of high and low energy, while the latter framework postulates some curvature modifications to the Lagrangian which are dominant at a certain cosmic epoch but dilutes at others. For example, Born-Infeld-like theories can lead to an effective description of the cosmic matter interpolating between dark matter and dark energy dominated epochs as a consequence of the Friedman equations, in terms of the Chaplygin Gas [3,4] or of the Anton-Schmidt fluid [5,6]. Other thermodynamically-motivated fluid models like the Dieterici [7] or the Shan-Chen [8] can as well exhibit a phase transition from a decelerating to an accelerating phase of the universe; the former from a matter-dominated epoch to a dark energy epoch, and the latter from an early radiation-dominated epoch to a dark energy epoch. The Shan-Chen model can also be used for describing the exponential expansion occurring during the inflationary epoch with the advantage of exhibiting a graceful exit mechanism, but for a different choice of the free parameters entering its equations of state than in the former analysis [9].
On the other hand, extended gravity theories in which a certain curvature invariant is added to, or used to replace, the Ricci scalar inside the Einstein-Hilbert Lagrangian can provide as well an evolution between different cosmic epochs as a consequence of the modified field equations themselves [10][11][12][13][14][15][16]. In spite of the correspondence between modified gravity theories and non-ideal fluid pictures (i.e. whose pressure and energy density are connected via P = w(ρ)ρ) [17], the former have the advantage of not violating some of the energy conditions which instead are broken when exotic fluids with negative pressure are adopted, and they preserve causality which would be lost when the adiabatic speed of sound squared becomes negative.
In this paper, we will merge the fluid and the modified gravity approaches and propose some cosmological models in which the gravity sector is accounted for by a Lagrangian of the type f (R) = R + ζR 2 , while the matter content is assumed to obey to some non-ideal equations of state with a well-established thermodynamical foundation known under the names of Redlich-Kwong, (Modified) Berthelot, and Dieterici fluid separately. The former assumption will allow us to account for the early-time dynamics, while the latter for the present-day epoch. Both these two models have been investigated separately in a number of literature works [18][19][20][21][22][23][24][25][26][27]. Here, we will obtain a cosmological dynamics with a rich variety of different behaviors like a non-singular bounce, two de Sitter-like epochs (thanks to the non-linear equation of state of the cosmic fluid in which w(ρ) is not a constant), possibly two radiation-dominated epochs, and possibly a phantom regime (the latter only in the Redlich-Kwong scenario). The comparison between three different realizations of the equation of state parameter function w(ρ) (for example which can be either always regular or admitting singularities, which can blow up or not at small or high energy densities, etc...) will give us the opportunity of enlightening which of our findings hold only when a particular fluid modeling is considered, and which instead seem to be a general characteristic of the cosmological dynamics. We must mention here that previously there have been some attempts to unify early and late time cosmology under certain forms of f (R) gravity [28][29][30]. However, it is worthwhile to remark also that the modifications utilized in those works are completely ad-hoc, lacking any motivation from the field theory point of view. The only modifications to the Einstein-Hilbert Lagrangian with some field theoretical motivations are the quadratic gravity theories. It has been known for some time that gravity Lagrangian containing additional quadratic curvature invariant terms are renormalizable [31,32]. Therefore in this work we do not intend to go beyond quadratic modifications. In particular we consider only the simplest case, namely, an R 2 correction term, along with fluids having a well-defined thermodynamic foundation.
We will tackle the technical difficulties arising in a fourth-order gravity theory like this one by adopting the set of dimensionless variables constructed in [33] which allows to cast the dynamical equations into a system of autonomous first-order equations suited for a dynamical system analysis. Such technique constitutes a powerful mathematical tool for describing the qualitative evolution of the the cosmological model under investigation not only in modified gravity [33][34][35][36][37][38][39][40][41][42][43], but also in multi-interacting fluid models [25,[44][45][46][47][48][49][50][51], and in exact or perturbed anisotropic and inhomogeneous cosmological models [52][53][54][55][56][57][58], just to mention a few examples. However, we will also propose a novel cosmologically transparent interpretation for those variables which was still lacking in the literature by deriving the physical restrictions they should obey to for avoiding tachyonic and ghost instabilities and connecting them to the cosmographic parameters, such as the deceleration, jerk and snap parameters which can be astrophysically constrained. Remarkably, we will show that such physical restrictions still allow the existence of a region in the phase space in which the energy density of the matter field is equal to the energy density of the curvature, which may be relevant for addressing the coincidence problem. Furthermore, our choice of variables will be useful also for showing that certain regions of the phase space are free from any of the five known types of cosmological singularities without the need of using the dominant energy balance formalism [59]. Moreover, after compactifying the phase space we will show that the region at infinity does not have only a geometrical meaning but it is such that the cosmic fluid equation of state reduce to the ideal behavior P ∝ ρ in which the interactions between the fluid constituents are suppressed.
One of the most severe shortcomings of the standard cosmological modeling is the Hubble tension, which is the discrepancy between the large and small scale estimates of the Hubble constant from supernova and cosmic microwave background data. Assuming that these predictions are not affected by any systematics, as to gravitational lensing effects on the cosmic microwave background angular spectrum [60] or to calibration and reddening issues for supernovae [61][62][63], an appropriate theoretical framework should be constructed for taming it. Several different proposals have been formulated, but none of them still seem fully satisfactory. For example the presence of a Proca field would reduce the Hubble tension [64], but there are no laboratory evidences of massive electrodynamic effects, and furthermore gauge invariance is lost in this theory [65]. Also, interactions between dark energy and dark matter may alleviate the tension [66], but thermodynamical considerations based on the Le Chatelier-Braun principle suggests that dark energy should decay into dark matter [67] while the fact that the structure formation era should precede the accelerating phase would require otherwise [68]. Our present work is intended as a rigorous dynamical study of a unified cosmic history model, combining two important frameworks one each from the study of early and late-time universe. Although we do not address the issue of H 0 tension here, an interesting scope for further investigation would be whether a unified cosmic history model, like the one we presented here, can provide an alternative to introducing ad-hoc interactions in the dark sector when it comes to alleviating the H 0 tension. Indeed this is not the first time that modified gravity and other ingredients are merged together. For example, anisotropic models in which the Copernican principle is relaxed have already been considered in Einstein-Aether gravity [69] also with a coupling to a scalar field [53], in braneworld cosmologies [70], or in f (R) gravity [71], just to cite a few examples. On the other hand, for a recent phenomenological proposal which may tame some observational challenges invoking two free parameters and requiring only a modification of the gravity sector in terms of a torsional Lagrangian see [72].
Our paper is organized as follows: in Sect. II we will review the field equations of the class of models we want to analyze and exhibit the equations of state of the cosmic fluids we are adopting mentioning their basic features, and we will as well introduce a formalism in which both curvature and matter effects are combined into an effective picture. Sect. III constitutes the main part of our work: in III A we will recast the equations governing the dynamics of our models as a system of autonomous first order equations in terms of a set of dimensionless variables on which we will also derive appropriate physical restrictions; in III B we will identify the cosmologically meaningful equilibrium solutions, explain for which ranges of the matter equation of state parameters they can arise pointing out possible bifurcations among them for particular types of matter contents, and report their stability showing that radiation-dominated, de Sitter-like and power law cosmologies can arise; in III C we will compactify the phase space and perform the analysis at infinity showing that a nonsinglar bounce occurs; in III D, III E and III F we will investigate the dynamics in the invariant submanifolds both numerically by plotting the trajectories in the phase spaces, by deriving analytically their stability, and by finding analytical results for the phase orbits in some specific cases; in III G we will relate the dimensionless variables we have adopted to the deceleration, jerk and snap cosmographic parameters which can be astrophysically measured. Then, in Sect. IV we will explain why some regions of the phase space are not affected by any cosmological singularity, and in Sect. V we will summarize the patterns that have emerged in our analysis by discussing which cosmological features we have discovered are sensitive to the particular modeling of the fluid, and which instead seem to be a general property. We will conclude in Sect. VI by discussing the cosmological relevance of our analysis and by putting the present work in the perspective of possible future projects. In Appendix A we review the applicability of the fluid models considered in this paper for the description of real gases beyond the cosmological context. The analytical computations of the stability of the isolated fixed points and of the invariant submanifolds are reported in the Appendices B, C, D, E which make use of both the standard notion of linear stability and of a much more advanced technique like the "center manifold analysis".

II. BASIC EQUATIONS OF QUADRATIC GRAVITY
The action of quadratic gravity in the Ricci scalar 1 reads as [78] with 2 f (R) = R + ζR 2 and κ = 8πG, G being the Newton's gravitational constant. ζ is a positive parameter quantifying the deviation of the quadratic gravity from the general relativistic Einstein-Hilbert Lagrangian at high curvature. These contributions are supposed to play an important role in the early universe driving the inflationary dynamics but diluting at later epochs [18][19][20][21][22][23][24]. This model constitutes a specific realization of a scalar-tensor theory of gravity because modifications in the gravity sector can be re-interpreted in the Brans-Dicke language as a new degree of freedom associated to a propagating scalar field [17]. Moreover, S m is the aggregate matter action responsible for all the fluid content of the Universe. In this paper we will assume the cosmic matter to be a perfect fluid (it is fully characterized by its pressure P and energy density ρ) obeying to a nonideal equation of state (pressure and energy density are not directly proportional to each other). To be more specific, we will consider some fluid models which constitute examples of evolving dark energy and/or unification of exotic and regular matter since in this latter case the sign of the pressure can change at different cosmic epochs as a consequence of the evolution of the energy density. Thus, our model is intended to study the evolution from inflationary to dark energy epoch by involving both quadratic corrections in the curvature and some nonideal fluid. Furthermore, in light of the Copernican principle, i.e. that the universe is homogeneous and isotropic, and considering an almost spatially flat universe, our geometrical model will be based on the spacetime Defining F := ∂f /∂R, and introducing the Hubble function H :=ȧ/a, where an overdot denotes a derivative with respect to the cosmic time, we can write the field equations for a flat Friedman universe under the action (1) as [10][11][12]: where we have adopted units such that κ = 1, and the Ricci scalar is related to the Hubble function via R = 6(ȧ 2 + aä) a 2 = 6(2H 2 +Ḣ) .
1 In this paper we will restrict ourselves to a modified gravity model quadratic in the curvature. However, other types of corrections have been proposed in the literature, either quadratic or beyond it, as in f (T ) theories with torsion [73], f (Q) with non-metricity [74], or f (G) with a Gauss-Bonnet term [75]. 2 In principle the most generic quadratic Lagrangian in curvature should also contain the terms R 2 αβ ≡ R αβ R αβ and R 2 αβγδ ≡ R αβγδ R αβγδ , which can be rewritten in terms of the Euler density E ≡ R 2 αβγδ − 4R 2 αβ + R 2 and the Weyl curvature invariant C ≡ R 2 αβγδ − 2R 2 αβ + R 2 3 . E does not contribute to the equation of motion due to the Gauss-Bonnet identity whereas C vanishes for FLRW metric [76,77]. Therefore the action (1) can be taken to be the most generic quadratic Lagrangian in terms of the curvature for a homogeneous and isotropic universe.
The field equations should be complemented by the Bianchi identitẏ ρ = −3H(ρ + P ) (5) which governs the energy conservation of the cosmic fluid. Furthermore, combining (3a) with (3b) we get which will be invoked in what follows for providing a transparent physical interpretation to the various quantities governing the cosmological dynamics. In fact, the joint effects of the matter content and of the modifications to the gravity sector can be combined into an effective total energy density and an effective total pressure which read as [13, Eq.(IV.82)]: Along this line of thinking, one can also define an effective equation of state parameter which encodes information about both the actual cosmic fluid and the curvature effects as For the the description of the matter content of the universe, we find convenient to follow the approach of [27] and consider the following modelings for the equations of state of the cosmic fluid separately: The fluid equation of state parameter defined as w := P/ρ takes respectively the forms: Therefore, our class of models is based on three free parameters (ζ, α, β). Different interplay between these free parameters will affect the existence of certain equilibrium configurations and certain types of finite-time singularities that we will classify in this paper with the purpose of constraining the values that these free parameters can assume by requiring these configurations to be physically meaningful. The two free parameters entering the equation of state of the cosmic fluid should be interpreted as: α > 0 is the temperature at which a thermodynamic phase transition occurs within the fluid, and it sets the strength of the interactions between the fluid particles since in the limit α → 0 all these equations of state describe an ideal fluid for which pressure and energy density are directly proportional to each other P βρ. This latter relation also shows the connection between β and the adiabatic speed of sound inside the fluid. The interested reader can find a more detailed review of the thermodynamic foundation of these fluid approaches in the Appendix of [27], and we will as well mention what the original reasons for their introduction for accounting for some features of real gases were in our Appendix A. More in general, these models try to provide a founded thermodynamical description of an evolving dark energy beyond ad hoc redshift parametrizations for helping its possible direct detection in the the far future. In fact, for accounting for both the Planck and weak lensing datasets, a redshift-dependent modeling of the dark energy equation of state parameter has been assumed in the form of w = w 0 + w 1 (1 + z) with w 0 and w 1 free parameters [82,Sect.6.3]. However, in this simple framework the analysis of the cosmic microwave background constraints on the distance to the last scattering surface is problematic, and therefore the refined Chevallier-Polarski-Linder parametrization w = w 0 + w 1 z/(1 + z) has been introduced [83,84]. The Barboza-Alcaniz w = w 0 + w 1 z(1 + z)/(1 + z 2 ) is another proposal which can be used in the whole redshift range z ∈ [1, ∞) [85]. Although these frameworks have been useful for studying the running of the dark energy potential beyond a cosmological constant, they do not try to establish the microscopic properties of such an exotic fluid which remain mysterious, calling for a physically deeper investigation. Finally, the functional w(ρ) can be interpreted also as an energy-dependent chameleon field [86,87].

III. QUALITATIVE ANALYSIS OF THE DYNAMICS OF QUADRATIC GRAVITY WITH NONIDEAL FLUIDS
In this section, we will derive the dynamical equations governing the evolution of the universe (2) in the gravity model (1) including some nonideal fluids by implementing the set of dimensionless variables considered in [33,36]. Particular attention will be devoted to the rewriting of the equation of state parameters (10) as functions of such dimensionless variables which are suited for a dynamical system analysis. Then, we will set some further constraints on the values of the free parameters of our model by requiring it to be free from instabilities. Lastly, we will list the mathematical equilibria and discuss their cosmological significance (which may provide tighter restrictions on the free parameters), possible bifurcations among them and their stability. Then, we will provide a prescription for compactifying the phase space with the purpose of investigating the dynamics at its infinity, and we will reconstruct the cosmological evolution on some invariant submanifolds also by analytically finding the equations of the phase orbits. In this section we also derive a set of relationships between the dimensionless variables employed in the dynamical system analysis and the observationally-relevant cosmographic parameters.
A. Derivation of the autonomous first-order dynamical system in terms of dimensionless variables The evolution equations to investigate in the R + ζR 2 gravity are: where we have obtained the first two by plugging (4) into (3a)-(3b). We can note that the first equation, which constitutes the Generalized Friedman equation, is not sensitive to the specific cosmic fluid modeling P = P (ρ), unlike the other two dynamical equations. Furthermore, eq. (6) can be rewritten in terms of the Hubble function as: 6ζ (2 ... H + 12HḦ + 9Ḣ 2 ) + 2(54ζH 2 + 1)Ḣ + 3H 2 + P = 0 .
Explicitly, for flat Friedman universes filled with the fluids (9a), (9b), (9c) evolving under the action of quadratic gravity, we get the following set of dynamical equations, respectively: • Redlich-Kwong fluid: • (Modified) Berthelot fluid: • Dietrici fluid: These differential equations are third order in the Hubble function (or equivalently fourth order in the scale factor), and non-linear in both the Hubble function and the energy density. Thus, it is convenient to tackle them by adopting dynamical system techniques and searching possible equilibrium configurations for clarifying their cosmological meaning and analyzing their qualitative dynamics [88][89][90][91]. Following the formalism of [33,36], we can recast these differential equations into a first-order autonomous dynamical system in terms of the following dimensionless variables: We introduce also the following auxiliary quantity: .
It is clear from (16) that these dynamical variables are undefined when H = 0. Therefore, this particular choice of variables pushes any possible fixed point corresponding to Minkowski solutions and bounce (a cosmological bounce is an alternative to the inflationary paradigm 3 ) or turnaround scenarios to the infinity of the phase space. Taking into consideration fixed points at infinity requires a global phase space analysis (see e.g. [37,38] in the context of f (R) gravity), which we will investigate separately in Sect. III C. Also, we do not expect any moment of maximum expansion at whichȧ = 0 = H since we are considering a flat ever-expanding universe filled with the effective fluid (8). However restricting to a domain of the full solution space consisting of only ever expanding (or ever contracting) solutions, this choice of variables is very advantageous when looking for a physical interpretation of the solutions and connecting with the cosmological observables. Therefore, the expansion normalized dynamical variables in (16) are appropriate for the consideration of this paper. The first-order autonomous dynamical system governing the evolution of the cosmological variables (16) is 4 : where ρ = ρ(x, y, z, Ω), N = ln(a(t)) denotes the number of e-folds of the universe [93], and where we have exploited the chain rule for any generic quantity χ = χ(t). From (16) we can write the Hubble function, its time derivative, and the fluid energy density in terms of the dimensionless variables as: which, together with the definitions (9a), (9b), (9c), allow us to rewrite the equation of state parameters as Furthermore, the Generalized Friedman equation (11a) is reduced to the constraint which should be used for removing one cosmological variable from the dynamical system (18). We choose to eliminate x for a twofold reason: the x-equation is apparently the most complicated one, and the w(ρ) can be naturally expressed in terms of (y, z, Ω) as done in (21). Keeping in mind eq.(17), the dynamical system (18) becomes: There are three physical viability conditions which should be accounted for when identifying the cosmologically relevant regions inside the full 3-dimensional y-z-Ω phase space 5 . They are the following: • Firstly, absence of ghost instabilities in f (R) gravity requires F (R) > 0, which implies 1 + 2ζR > 0 for our scenario [94]. From (4)- (20) we can write so that the absence of ghost instabilities requires which can be satisfied for These conditions represent two disconnected regions on the first and third quadrant of the y-z plane bounded by the line y = 2z and the z-axis.
• Secondly, absence of tachyonic instabilities for a generic f (R) gravity theory requires f (R) > 0, which in our case simply implies ζ > 0 [94]. From the definition of the dynamical variables (16), we note that 5 To the best of our knowledge this is the first time that these physical viability conditions are used to constrain the viable region of the phase space spanned by the expansion normalized variables (16).
• Finally, the weak energy condition requires the energy density to be locally non-negative: We can observe also that R is non-negative within the semi-infinite y ≥ 0 region. Therefore, where the last equality follows from the observation that in the General Relativity limit (which corresponds to ζ → 0 in the system (11)), one recovers the usual Friedman equation To summarize, there are two disjoint regions of the phase space which are physically relevant: These are two distinct semi-infinite wedge-shaped sectors above the Ω = 0 plane, in the first and third quadrants of the y-z plane, respectively. The region in the first quadrant is confined between the two lines y = z and y = 2z, while the region in the third quadrant is confined between the line y = z and the z-axis. We stress that till now we have not included the boundaries of these two regions (which are represented by equalities rather than inequalities in (31)), because they require a more careful treatment. The plane defined by the equality y = z accounts for the General Relativity limit R + ζR 2 ≈ R (which can be expressed as ζ → 0 thanks to Eq. (27)) in which the quadratic modification in the Lagrangian is negligible with respect to the Einstein-Hilbert contribution. It is not appropriate to consider the plane y = z in the analysis that follows because the dynamical variables are undefined there and the dynamical system formulation that we are adopting becomes singular on the plane y = z. However, this does not prevent the origin (y, z) = (0, 0) to be describe a physically meaningful configuration, as it can be appreciated from lim y→0,z→0 The dynamical system (23) is therefore singular everywhere on the y-z plane except along the line y = z = 0. The other boundary of the acceptable region in the first quadrant is the plane defined by the equality y = 2z, while for the one in the third quadrant is the z-Ω plane defined by the condition y = 0. The plane y = 2z corresponds to the limit R + ζR 2 ≈ ζR 2 (which is equivalent to ζ → +∞, as it can be seen from Eq.(25)) which occurs when the quadratic modification term in the Lagrangian becomes dominant over the Einstein-Hilbert contribution. At this stage, both the planes y = 2z and y = 0 can be safely included in the physically viable region of the phase space, which is thus given by As a consistency check, one can note from (20) that H 2 > 0 in both these regions. The dynamical system (23) admits the two invariant submanifolds y = 0 and Ω = 0. An invariant submanifold divides the entire phase space into two distinct regions on its both sides. Although they can at most reach the boundary, no phase trajectory can cross the invariant submanifold leaving one region and entering the other. Also, any orbit originating from a point on an invariant submanifold will always remain on that submanifold signifying that the reconstruction of the dynamics requires knowledge on the initial data. We note that R is always negative in the third quadrant of the y-z plane (because y is negative), and thus this region cannot contain any fixed point interpreted as a De-Sitter cosmology (or any other cosmology with negative deceleration parameter). Therefore, observational datasets suggest that the cosmological evolution should not occur in this region of the phase space.
B. Qualitative dynamics: equilibria, stability, and bifurcations In Table I we exhibit all the equilibrium points that arise mathematically for the dynamical system (23), along with the corresponding cosmological solution they represent, if any. In fact, some of the fixed points should be ignored on physical and observational grounds: 1. The point P 3 is unphysical for the scenario of a universe filled with the (Modified) Berthelot fluid because it violates the weak energy condition. Should we consider the Redlich-Kwong fluid, it can carry a physical interpretation for 1 ≤ β ≤ 17 9 . Interestingly, the former point corresponds to the β = 0 case of the latter.
Cosmic fluid Fixed point yeq zeq Ωeq w eff Cosmology  (23) once Eqs. (21) have been implemented for the different types of cosmic fluids. The effective equation of state parameter w eff is computed from (16b) and it takes into account the contributions of both the actual matter content and of the curvature effects (see also (8)). We refer to the main text for a detailed explanation of why some equilibria do not represent any meaningful cosmological model.
2. Any orbit that approaches the point P 4 must reside inside the third quadrant of the y-z plane in which the deceleration parameter is always positive. Therefore, this point should be ignored on observational ground.
3. Similarly for the state P 5 : we can note that 0 ≤ Ω eq ≤ 1 delivers a negative z eq implying that any orbit approaching P 5 must reside within the third quadrant of the y-z plane. Therefore, this point should also be ignored on observational ground. 4. The fixed point P 7 is unphysical for a universe filled with the (Modified) Berthelot fluid because it violates the energy condition Ω eq ≤ 1. In the Redlich-Kwong scenario it is physical for The conditions on the model parameters which should be imposed for endowing the remaining mathematical solutions reported in Table I with a cosmological interpretation are listed in Table II. They follow by imposing 0 ≤ Ω ≤ 1. It should be appreciated that this affects only the range of validity of β, while no further constraints other than the already discussed are arising for α and ζ. Among the physically viable fixed points we can identify three distinct types of cosmological solutions: 1. De-Sitter-like cosmology: There are two different possible realizations of a De-Sitter-like cosmology 6 , which are represented by the isolated fixed points P 1 and P 2 . The equilibrium P 1 always constitutes a physical configuration for all the three fluids, whereas P 2 is relevant in cosmology only imposing certain constraints on the model parameter β as shown in Table II. For all the three types of matter, the ideal fluid regime (α → 0) leads to a saddle-node bifurcation between P 1 and P 2 . Furthermore, in the case of the Redlich-Kwong fluid model, also the equilibrium P 3 reduces to P 1 if we fix β = 1. In this latter case a pitchfork bifurcation is possible if we choose simultaneously α = 0 and β = 1 [95].
2. Power law evolution: There are up to two different possible realizations of the power law evolution (a ∼ t 1/2 ), which are represented by the isolated fixed points P 6 and P 7 . P 6 is always physical for all the three fluids whereas P 7 is relevant for cosmology only in the Redlich-Kwong scenario and restricting Big-Rip singularity: The fixed point P 3 , which is physically well-defined only considering the Redlich-Kwong fluid, represents a big-rip singularity which is asymptotically approached at the finite time 7 In fact, we can note that the scale factor is diverging by looking at its time evolution; the energy density is also diverging because of (20) and taking into account that y = 2z = 0, and this comes also with a divergence in the pressure because of the form of the equation of state (9a). Therefore, all the conditions for the occurrence of a Big-Rip singularity are fulfilled [96][97][98][99]. As the limiting case of β = 1 is approached, which we showed corresponds to a bifurcation with the De-Sitter-like cosmology, the time at which this singularity occurs is shifted at infinity. For β = 1, t s ∼ 1/H 0 and the singularity time is comparable to the age of the Universe. Keeping in mind the parameter range in Table II, we can also see from Table I that both the effective and the matter parameters w eff , w < −1 for this point, where for the latter w = −β as from (21a). Therefore this fixed point also corresponds to a phantom dominated phase at which the Redlich-Kwong fluid itself behaves like a phantom fluid. Furthermore, the adiabatic speed of sound for the Redlich-Kwong fluid, which can be computed from (9a), once specified to P 3 via (20) delivers which is smaller than −1 in the range of interest of β. 6 Here, by De-Sitter-like cosmology we mean a cosmology in which the Hubble function is constant. From the general system of equations (11) we can note that also Minkowski can constitute an equilibrium solution when we consider the Redlich-Kwong, (Modified) Berthelot, and Dieterici fluid models; this would correspond to the particular case of H = const. = 0 (and ρ = P = 0). However, the dynamical variables (16) are ill-defined for a Minkowski solution; we will address this limitation by compactifying the phase space in Sect. III C. We also remark that not all the fluid models currently adopted for a dark matter -dark energy unification are compatible with the Minkowski spacetime being an equilibrium solution, with the (Generalized) Chaplygin Gas and the Anton-Schmidt proposals being some examples; see discussion in [25]. 7 For computing ts, note that and that we fixed a(t = 0) = 1.

P6
Always exists Always exists Always exists Unphysical Does not exist  Table I from mathematical to physical are derived demanding 0 ≤ Ω ≤ 1. The limits α → 0 and ζ → 0 correspond to ideal fluid and General Relativity, respectively. The points P 4 and P 5 are not included in this Table because they belong to a region of the phase in which the deceleration parameter is always positive.
The stability nature of the fixed points is listed in Table III and detailed calculation is presented in Appendix B. It is possible to note that under the assumption that α, ζ > 0, only the parameter β, which is related to the adiabatic speed of sound within the fluid, affects the stability nature of the finite isolated fixed point.

C. Phase space analysis at infinity
Compactification of an unbound phase space is necessary to search for any possible fixed point that lies at its infinity: thanks to this procedure the fixed points at infinity are mapped to the boundary of the corresponding compact phase space. In general all the dynamical variables can tend to infinity, which means the phase space of the theory can exhibit a unlimited extent in all the directions. There are different prescriptions for f (R) cosmologies (see e.g. [38] for a generic f (R) theory and [37] for the particular R + ζR n theory) for compactifying the phase space in all the directions. However, in this Sect. we introduce a new compactification technique which directly exploits the physical viability conditions we previously derived in (33). As we will show below, one can use these constraints to define some invariant submanifolds that border the physically viable region of the phase space and then we are left with only one direction in which the phase space need to be compactified.
From a mathematical point of view, the dynamical system (23) is singular on the plane y = z. Since this plane is one of the boundaries of the region of the phase space we are interested in, this singularity can be regularized by introducing a new time variable τ such that in terms of which the dynamical system can be re-written as Now one can write which show that the planes y = z and y = 2z are invariant submanifolds as well. As discussed in Sect. III A, these two planes are equivalent to the two limits ζ → 0 and ζ → +∞ respectively. To the best of our knowledge this is the first time that the physical viability conditions which follow from the absence of ghost and tachyonic instabilities are recast as invariant submanifolds on the phase space of quadratic gravity. Linear stability analysis reveals that the invariant submanifold y = z is always attracting whereas the invariant submanifold y = 2z is attracting (repelling) for ; detailed mathematical analysis is given in Appendix D.
Before proceeding any further, it is important to comment that the dynamical system in Eq. (38) should not be used to determine the fixed points, because time redefinitions like (37) may introduce artificial solutions which are not appearing in the original dynamical system. For example, one can notice that the system in Eq. (38) has two lines of fixed points given by and both of which do not occur in the original dynamical system (23). These fictitious fixed points are a pure mathematical artefact due to the time redefinition (37). We stress that this and the following steps are purely mathematical treatments aimed towards compactifying the phase space by introducing appropriate invariant submanifolds. All the finite fixed point analysis should be carried out before these steps. Along with Ω = 0, the physically relevant region of the phase space is therefore bounded by three invariant submanifolds. Since in this region the dynamical variable Ω is itself bounded (0 ≤ Ω ≤ 1), as demonstrated in Sect. III A, one needs only to compactify the radial direction in the y-z plane. For achieving this goal we first switch to plane polar coordinates in the y-z plane y := r cos θ , z := r sin θ , subject to the restrictions The dynamical system (38) in terms of the r-θ-Ω variables (41) becomes where the fluid equation of state parameters (21) entering the latter equation are given by w(r, θ, Ω) = 2βζr(2 sin θ − cos θ) 2 2ζr(2 sin θ − cos θ) 2 + αΩ(cos θ − sin θ) w(r, θ, Ω) = 2βζr(2 sin θ − cos θ) 2 4ζr(2 sin θ − cos θ) 2 − αΩ(cos θ − sin θ) exp 2 − αΩ(cos θ − sin θ) ζr(2 sin θ − cos θ) 2 (Dietrici) . (44c) As we have previously remarked, the introduction of the artificial line of fixed points L 1 ≡ (r = 0) is clearly confirmed by inspecting the system in Eq. (43). We should remove this fictitious fixed point by another time redefinition so that the dynamical system becomes dr dτ * = r r cos 4 θ + The radial direction can be compactified by introducing the new compact variable [36,39,40] so that r = 0 coincides with R = 0 and r = ∞ is mapped onto R = 1. In terms of R the dynamical system to investigate is with We can note that all the three fluid equations of state remain well-behaved at infinity, i.e. have finite limits as R → 1. The latter dynamical system has a pole at R = 1, i.e. is apparently singular at the boundary. This can again be eradicated by defining a new time variable η as Therefore, the dynamical system governing the evolution of the compatified variables can be written as Since only the r-direction can be infinite, all the asymptotic fixed points should correspond to r → ∞ (or R → 1). Therefore we need to identify the fixed points in the R-θ-Ω phase space which fulfill R = 1. Setting R = 1 in (51), we obtain Interestingly, the evolution at spatial infinity is not explicitly sensitive to the modeling of the cosmic fluid as it was observed in the case of R n gravity [36] because w does not enter anylonger the dynamical system (however we remind that we have used previously our particular equations of state for checking that they well behave at infinity). A further information that can be obtained from the analysis at infinity is that R → 1 is an invariant submanifold. To determine the cosmology corresponding to this submanifold first we note that using (20) one can write within the range tan −1 1 2 ≤ θ < π 4 .Ḣ is positive at all points on this hypersurface whereas H, ρ vanish. This is exactly the condition for a matter-less nonsingular bounce. We remark that had we not compactified the phase space, we would have not been able to discover this bounce solution in our cosmological models for the reasons discussed below eq. (17). Keeping in mind the range of θ given in (42), asymptotic dynamical analysis reveals the following features: • The asymptotic invariant submanifold accounted for by R = 1 is a repelling submanifold. Detailed calculation regarding the stability of this submanifold is presented in Appendix D. Therefore the nonsingular bouncing solutions that lie on this submanifold may constitute past epochs of the universe.
• The point P i ≡ (R, θ, Ω) = (1, tan −1 1 2 , 0) is an isolated fixed point at infinity. This fixed point, although represents a nonsingular bounce, does not necessarily need to be matter-less, as at this point tan θ = 1 2 . In fact, as was pointed out in [92], matter-less nonsingular bounce in f (R) gravity requires the equation RF (R) − f (R) = 0 to have a positive root R b , which is not satisfied in case of R + ζR 2 gravity. Linear stability analysis reveals P i is a saddle point. Stability calculation is presented in Appendix E.

D. Evolution on the y = 2z invariant submanifold
The submanifold y = 2z corresponds to the limit ζ → ∞ which accounts for the high energy regime in which the gravitational field is so strong that the theory is dominated by the R 2 term. From (39b) it is seen that y = 2z is an invariant submanifold for the dynamics of the system. This submanifold is of attracting nature for (detailed calculations presented in Appendix D) and of repelling nature for In terms of the variables r-θ-Ω (or R-θ-Ω for the compact case) this submanifold corresponds to θ = tan −1 (1/2). From (44) or (49), for the fluid equation of state parameter we get lim  Figure 1 for the cases of equations of state corresponding to dark matter (e.g. pressureless dust), stiff fluid and a cosmological constant. We remark that stiff fluids are canonically equivalent to massless scalar fields [100], and some cosmological models indeed predict an epoch of the universe in which they are the dominating energy content [101,102]. On this invariant submanifold the dynamical equations can be reduced to: For the case of stiff matter, we can find the orbit in the phase space by solving the differential equation which delivers the implicit solution where J 1 is a constant of integration. The quantity J 1 (z, Ω) is conserved along a particular orbit, but has different values for different orbits, and therefore it can be interpreted as the total "energy" of the Universe. The cosmological evolution must respect the principle of energy conservation: we can interpret eq. (59) as a sort of "energy conservation equation" which is providing a law describing how the energy of the cosmic fluid accounted for by Ω is converted into the "geometrical energy" accounted for by the Ricci scalar R; this result is especially relevant for the description of the inflationary epoch in which the quadratic term in the curvature is dominating. Furthermore, in the case of a stringy fluid with w = − 1 3 , which may describe some topological defects or monopoles arising in the early universe [103], by integrating the differential equation we obtain the implicit orbit equation we focus our attention on the viable region bounded between the lines 0 < z < y < 2z. This analysis provides a graphical confirmation that the submanifold R = 1 acts as a source for the dynamics, and that the cosmic evolution is indeed contained within the physical region.
where J 2 is a constant of integration. Also for the radiation case w = 1 3 it is possible to integrate analytically the evolution equation and we obtain the implicit orbit equation where J 3 is another constant of integration.
E. Evolution on the Ω = 0 submanifold It appears either from (23c) or from (48c) that the plane Ω = 0 is an invariant submanifold for the cosmic dynamics. Taking into account that the physically viable region is constituted by the wedge 0 < z < y < 2z, we depict the phase orbits in this invariant submanifold in Figure 2 by using the evolution eqs. written in polar coordinates (51a)-(51b). In this way we can get a graphical confirmation that the dynamics is indeed bounded inside this region and that the boundary at spatial infinity R = 1 acts as a source for the cosmic dynamics containing possible past epochs of the universe. Unlike the case of the invariant submanifold y = 2z discussed in Sect. III D, the dynamics on the invariant submanifold Ω = 0 does not depend on the particular modeling of the cosmic fluid. However, the stability nature of this invariant submanifold is sensitive to the value of the parameter β as demonstrated in Appendix D, and more in detail it is attracting (repelling) according to 2 − 3β − 3y + z < 0 (> 0) for the Redlich-Kwong and (Modified) Berthelot fluids and 2 − 3e 2 β/2 − 3y + z < 0 (> 0) for the Dietrici fluid.
F. Evolution on the R = 1 submanifold R = 1 is an invariant submanifold at the infinity of the phase space. We can find the equation for the orbit J = J(θ, Ω) at the infinity of the phase space by solving the partial derivative equation Implementing (52) we find where F can be any arbitrary function. For reasons of mathematical simplicity, we choose: We note that the quantity J(θ, Ω) is a positive quantity within our range of θ, which is conserved along a particular orbit but can have different values for different orbits. This quantity can again be interpreted as the total "energy" of the Universe and the cosmological evolution must respect the principle of energy conservation. Therefore, the orbits on this submanifold are a family of curves obeying to the equation where J is a constant. We stress as a consistency check that the same result also follows by integrating a differential equation for dΩ dθ derived by dividing side by side (52c) with (52b). In terms of the original dynamical variables one can write the equation of the orbits as Finally, by using (16) this condition can be recast in terms of the energy density, of the Hubble function and of its first derivative as: where we have introduced the new constantJ This result allows us to confirm independently what written below eq. (53a): since H = 0 andḢ = 0 we get that the submanifold R = 1 corresponds to a matterless cosmological epoch. However, this should not be taken naively to imply that Ω = 0 because this latter quantity comes with a factor H in the denominator and indeed this is true only on the hypersurface y = 2z as it can be understood from (68).

G. Cosmographic analysis
We will now discuss some observational properties of the universe in correspondence of the physical equilibrium points listed in Table I by computing the corresponding three cosmographic parameters, namely the deceleration, jerk and snap parameters [104,105]: The cosmographic parameters are connected to the luminosity distance via [104,[108][109][110][111]: and to the cosmic history of the universe as: where a subscript '0' denotes that the quantity has been evaluated at the present time. In this Sect. instead we will estimate the cosmographic parameters characterizing the relevant equilibrium configurations. We exhibit our findings in Table IV. We will achieve this goal by recasting the dimensionless cosmographic parameters q, j, and s in terms of the dimensionless variables introduced in (16). Using the inter-relations between the cosmographic parameters (72), we can write Calculating the right hand side of the above equations using the dynamical evolution (23), we can provide explicit expressions for the cosmographic parameters in terms of the phase space coordinates: These expressions can be directly generalized to include the fixed points at the infinity of the phase space by switching to the compact phase space coordinates Substituting in Eq.(76) we get the following explicit expressions for the cosmographic parameters in terms of the compact phase space coordinates: First of all, we easily get that (76a) implies that y =constant submanifolds correspond to cosmic moments with the same value of the deceleration parameter. Possible Minkowski solutions necessarily lie on y = 1, and therefore our models do not contain them as equilibrium configurations (this resolves the ambiguity whether the De-Sitter-like cosmologies we have identified in Sect. III B can come with H = const. = 0). The expression of the cosmographic parameters in terms of compact variables also allows us to show that the cosmographic quantities are diverging at spatial infinity of the phase space which is consistent with having a bounce there characterized by d L → ∞. We would like to mention that a cross-check procedure for computing the jerk parameter which does not rely on inter-relations is the following. We implement (20) into (16a) and then solve for the second time derivative of the Hubble function: from which x can be eliminated thanks to the constraint (22): H = (9 + z − Ω − y)y 2 − 8(z + 2)y + 16z 72ζy(y − 2z) 6(z − y) ζy(y − 2z) .   Table I. We refer to the main text on details about the mathematical steps involved in these computations. We remark that for a correct interpretation of these results it is necessary to take into account the appropriate range of validity for the parameter β for each equilibrium point separately, as summarized in Table II.
Finally, the jerk parameter is obtained just by algebraic manipulations. We get: Interestingly, the jerk parameter is regular on y = 2z because the divergence inḦ has been cured by the likewise divergence in H. For estimating it on y = z = 0 it is appropriate to choose a different set of variables taking into account that in such case we fall back in the General Relativity framework. The values we get for the deceleration parameter imply that phase transitions between epochs in which the expansion of the universe is accelerating and decelerating are allowed in our class of models. In particular, at least one equilibrium point comes with q > 0 and at least two with q < 0 for each fluid model. Comparison between available astrophysical datasets and the predicted values of the cosmographic parameters can constrain the theory parameters of f (R) theories [107,112]. A cosmographic interpretation of the Gold SNeIa dataset suggests that q 0 −0.90 and j 0 2.7 [113,114]. It should be noted that due to the presence of w(y, z, Ω) in the expression for the cosmographic parameter s (76c), the present-day epoch would correspond to different triples (y, z, Ω) in the phase space. However, the phase space point representing today universe is located in the region y < 1. Information on physically relevant trajectories in the phase space can therefore be obtained by noticing that from the expression of the jerk parameter in terms of the dynamical system variables (81) we get implying that the jerk parameter is an increasing function with respect to Ω and decreasing with respect to z.

IV. SINGULARITIES CLASSIFICATION
In this section we will investigate the possible occurrence of finite-time singularities in the class of Friedmannian f (R) cosmologies we have previously introduced for clarifying whether the different modelings of the cosmic fluid and the modifications beyond general relativity to the gravity sector affect them. In what follows we will denote with t s the time at which a singularity may occur. Applying a literature scheme [96,97], we will be interested in the following five different possible types of singularity: 3. Big freeze singularity or Type III is characterized by lim t→ts a(t) = a s , lim t→ts ρ eff (t) = ∞, lim t→ts |P eff (t)|= ∞ [119]; 4. Generalized sudden singularity or Type IV is characterized by lim t→ts a(t) = a s , lim t→ts ρ eff (t) = ρ s , lim t→ts |P eff (t)|= P s , lim t→ts H (i) (t) = ∞, i = 2, ... [117,118,120,121]; 5. w singularity or Type V is characterized by lim t→ts a(t) = a s , lim t→ts ρ eff (t) = 0, lim t→ts |P eff (t)|= 0, lim t→ts w eff = lim t→ts P eff (t) ρ eff (t) = ∞ [122,123].
In this classification, we have denoted with a s , ρ s and P s some finite constant values of the scale factor, the effective energy density and its corresponding pressure at time t s . We recall that in our analysis we will assume positive α and ζ, while we will not make any assumptions on the sign of β. We also remark that we are working with the effective values of the energy density, pressure and equation of state parameter which encode information both on the actual matter fluid and the curvature effects, as done for example in [124][125][126][127]. Before analyzing the possible occurrence of a finite-time singularity in a generic point of the phase space, we investigate the situation in correspondence of the isolated fixed points reported in Table I. By looking at the evolution of the scale factor, they can exhibit three different types of cosmological evolution: de Sitter-like (P 1 and P 2 for all the three types of fluids), radiation (P 6 for all the three types of fluids, and P 7 for Redlich-Kwong), and power-law (P 3 for Redlich-Kwong).
• The de Sitter-like cosmologies do not correspond to any finite-time singularity because the effective energy density, pressure and equation of state parameter are finite constants.
• In the case of an "effective" radiation domination, the scale factor (a ∼ t 1/2 ) would approach a s = 0 at the time t = 0 in correspondence of which ρ eff , P eff ∼ 1/t → ∞, and therefore a finite-time (recalling that the present-day time is t 0 > 0) Type III singularity occurs in the past.
• The isolated fixed point P 3 in the Redlich-Kwong scenario can correspond to a Type I singularity occurring at a finite time t s (34) in future if 1 < β ≤ 17/9. We note that in P 3 which diverges also for β → 1; however this does not imply a finite-time singularity as can be seen from eq. (34).
We will now investigate whether some type of finite-time singularity can occur in some other regions of the phase space. By using the definition of effective energy density (7a), and the relationships between the Hubble function and the dimensionless variables (20), and (41), we have Furthermore, by using eqs. (16b)-(41) we can get the effective equation of state parameter defined in (8), and pressure in terms of dimensionless variables as: First of all, we note that on the planes y = 0 and y = 2z, both the effective energy density (84) and effective pressure (86) are diverging, so that two of the requirements for having either a Type I or a Type III singularity are fulfilled. We also remark that in these regions of the phase space both the Hubble function and its first derivative are diverging, as we can understand from eq. (20), and therefore we have a true curvature singularity in which the Ricci scalar (4) is blowing up. More in detail, everywhere on the plane y = 0 the effective fluid behaves like radiation, implying a Type III singularity since a ∼ t 1/2 (see also the equilibrium points P 6 in Table I for all the three types of fluids, and P 7 for Redlich-Kwong); this implies also that both energy density and pressure are diverging as ρ eff , P eff ∼ H ∼ t −1 ∼ 1/a 2 ∼ (1 + z) 2 (where this latter z denotes the redshift). Therefore, assuming that the present-day is at the finite-time t 0 > 0, a Type III singularity occurs in the past at the time t = 0.
For understanding the behavior of the singularity on the line y = 2z, we recall that a Type I singularity would require w eff < −1 [128], i.e. y > 2. Therefore the plane y = 2 separates the line y = 2z into two parts on whose sides a Type I or a Type III singularity can occur; this finding is consistent with the evolution of the scale factor exhibited in Table I, and the previous discussion about the equilibrium point P 3 for the Redlich-Kwong fluid.
We can provide a rough estimate of the time t s at which these singularities occur by approximating y ≈ y s in a small neighborhood of the line y = 2z assuming that the present-time t 0 configuration is contained there. This implies that d(1/H) dt ≈ 2 − y s . Thus, H(t) ≈ H0 1+(2−ys)(t−t0)H0 which diverges at t s ≈ t 0 + 1 (ys−2)H0 showing that the Type I singularity would be a future singularity, while the Type III a past singularity.
On the other hand, for having a finite energy density, but a diverging pressure we would need a diverging equation of state parameter. By looking at (85), we see that this is possible at and only at infinity, that is for r → ∞. In fact, in such a regime, by using eq. (86) we get which can diverge if and only if θ = arctan(1/2). Thus, a Type II singularity may occur only at the point P i . Moreover, in Sect. III C we have showed that H = 0 there, i.e. we have a well-behaving de Sitter-like scale factor and a finite (zero) effective energy density fulfilling all the conditions for having a Type II singularity. We remark that should we have considered the pressure of the actual matter fluid only, a Type II singularity may have arisen in the Dieterici framework only [25]. Moreover, by looking at the second time derivative of the Hubble function in terms of the dimensionless variables given in eq. (80), we see that a Type IV singularity may occur either along y = 0 or along y = 2z. This is the mildest possible singularity because it does not imply geodesic incompletness nor diverging curvature scalars. However, in these regions of the phase space also the energy density is diverging as it can be understood from eq. (84) violating (at least) one of the requirements in the definition of a Type IV singularity; as previously discussed also the Ricci scalar is diverging in such circumstances violating the conditions for a Type IV singularity. Interestingly, this analysis shows that the effective energy density and pressure arising from gravity modifications cannot mimic those of linearly interacting dark matter -dark energy where the latter is modeled according to the Redlich-Kwong or the (Modified) Berthelot fluid, as in those cases a type IV singularity is allowed for certain strengths of the coupling term [25].
Finally, by looking at (85) we see that a Type V singularity may occur only at spatial infinity for which r → ∞; this would be consistent with having also a diverging deceleration parameter there as we have found in Sect. III C. Then, by recalling (87) we see that the effective pressure can vanish if and only if θ = π 4 . Under these assumptions also ρ eff = 0, and taking into account the discussion of Sect. III F we further have a finite scale factor fulfilling all the requirements for a Type V singularity. This result follows from the gravity modifications and constitutes an important difference than General Relativity in which a Type V singularity has been excluded for the three types of Redlich-Kwong, (Modified) Berthelot and Dieterici fluids [25]. In fact, we can observe that such type of singularity persists also in the limiting case of ρ, P → 0, i.e. of absence of an actual cosmic fluid.

V. DISCUSSION ON GENERIC BEHAVIOR
Out of the global dynamical analysis of the system that we have presented in this paper, we note that the finite fixed points P 1 , P 6 and the asymptotic fixed point P i always exist for all the three fluids irrespective of whatever values we choose for the model parameters α, β, ζ, whereas all the other fixed points either exist for a certain fluid or for a specific range of values for the model parameters, and coincide with either P 1 or P 6 for certain values of those model parameters. The fixed points P 1 , P 6 and P i therefore characterize some generic features of the cosmological model in quadratic gravity consisting of the three fluids under consideration. We note that all these three fixed points lie at the line of intersection of the planes Ω = 0 and y = 2z. We stress that Ω = 0 does not necessarily imply a vacuum solution if either r = 0 or y = 2z, so that these three fixed points, although lying on the Ω = 0 plane, should not necessarily correspond to vacuum solutions of the R + ζR 2 gravity theory. Another point to note is that, as we had discussed before, the plane y = 2z corresponds to the limit ζ → ∞, so that the points lying on this plane can be interpreted to be the solutions of f (R) = R 2 theory of gravity. As shown in [36], irrespective of the fluid under consideration, the phase space of R n (n ≥ 2) gravity is always 2-dimensional, which is consistent with our interpretation.
Below we explicitly point out the generic dynamical features of the scenario that we have considered.
• P 1 is a De-Sitter solution that lies on the line of intersection of the planes Ω = 0 and y = 2z. This point represents the exact De-Sitter solution of R 2 gravity, which is the basis of Starobinski's inflationary scenario [78]. Since it is always a non-hyperbolic fixed point one needs to do a center manifold analysis to determine the stability, which is done in Appendix C. From Eq.(C5) we see that two of the eigenvectors of the Jacobian at that point lie on the Ω = 0 plane. The eigenvector corresponding to the negative eigenvalue is along the line (y = 2z, Ω = 0), which implies that the De-Sitter solution in R 2 gravity is an attractor. The eigenvector corresponding to the zero eigenvalue is along the line y + z = 3, and the center manifold analysis reveals that the dynamics is always away from the fixed point along this direction. In the complete R + ζR 2 theory, this corresponds to an exit from the De-Sitter phase.
• P i is a nonsingular bouncing solution (H = 0,Ḣ > 0) as discussed in Sect. III C. As demonstrated in Appendix E, this point is a saddle: repelling in the direction normal to the surface R → 1 and attracting in the directions normal to the planes Ω = 0 and y = 2z. The trajectories flowing from P i to P 1 can be interpreted as early universe solutions with an inflationary phase following a nonsingular bounce 8 . The flow at P 1 away from it along the line y + z = 3 in this case corresponds to the "graceful exit". This is consistent with the well known result that Starobinski's inflationary scenario is a transient attractor in R + ζR 2 gravity [24].
• P 6 is an "effective" radiation dominated phase (w eff = 1 3 ). The trajectories flowing from P 6 to P 1 can be interpreted as late time solutions with a transition from a radiation dominated epoch to a late time accelerating epoch corresponding to dark energy domination. The flow at P 1 away from it along the line y + z = 3 in this case implies an end to the accelerated phase of expansion, which, in GR, is possible only if the cosmological constant changes sign.
Apart from these generic features, there are some other interesting points worthwhile for explicitly commenting upon: • An interesting thing to note is that the same fixed point P 1 can be interpreted as either an inflationary epoch or a late time acceleration epoch, depending on which of the phase trajectories we choose to consider.
• It is also worth mentioning here that we do not get any fixed point corresponding to a matter dominated epoch because we have not considered any dust fluid that may correspond to the CDM. A matter dominated epoch requires w eff = 0 or equivalently y = 1 2 . We note that although any trajectory flowing from P 6 to P 1 crosses the plane y = 1 2 , there is no actual fixed point with y = 1 2 , and therefore no matter dominated "phase" in the picture. It is however interesting that an "effective" radiation-like epoch is arising even without explicitly including any ultra-relativistic fluid in the picture.
• As clear from Tables II-III, for specific ranges of the values of the model parameter β, the fixed points P 2 and/or P 3 can exist and can also be stable. In such cases there might be more than one De-Sitter phases in the complete evolution of certain cosmological solutions. The trajectories that encounter two De-Sitter fixed points (P 1 and P 2 or P 3 ), with P 1 being saddle and P 2 (or P 3 ) being stable, are particularly interesting. It should also be noted that P 2 (or P 3 ), when exists, can only be reached after P 1 . For such solutions P 1 can represent Starobinski's curvature driven inflation, whereas P 2 (or P 3 ) can represent a future attractor corresponding to the late time acceleration.
• It is worthwhile to note that the other two model parameters, namely α and ζ, do not affect neither the existence nor the stability nature of the fixed points, as long as they are assumed to be positive. These two parameters quantify the deviations from ideal fluid and from GR respectively. Existence and stability of fixed points depend only on the model parameter β, which characterizes the equation of state parameter of the fluid in its ideal limit. The parameter α however is crucial in relation to the bifurcation of the De-Sitter fixed points. It is precisely the non-ideal nature of the fluid (α = 0) that makes it possible to obtain two separate De-Sitter fixed points P 1 and P 2 , hence providing a scope for describing the early and the late time De-Sitter epochs at one go.
• The only case in which a big-rip singularity can arise in finite future is for the Redlich-Kwong fluid with β > 1. In this case the De-Sitter fixed point P 2 is a saddle, implying that the late time De-Sitter phase is an intermediate cosmological phase and not an attractor. In this particular case the true future attractor is P 3 , which is a big-rip singularity.
The generic features and other interesting points listed in this section are the take home messages from our present study.

VI. CONCLUSION
In this paper, we have investigated some cosmological models governed by a modified Friedman and a modified Raychaudhuri equation (11) equivalent to the following algebraic relations between the cosmographic parameters 9 : which can be summarized into the a single expression in which the parameter ζ does not enter directly: Whether this evolution of the rate of expansion can tame some of the problems related to the Hubble tension [131][132][133] is beyond the purpose of the present paper, but already at this stage we have demonstrated that these models come with many desirable features: they exhibit an inflationary epoch admitting a graceful exit, a radiation dominated epoch in which light elements may form [134], and a late-time De Sitter epoch consistent with supernovae observations [135,136]. Furthermore, more than one De Sitter epoch in the cosmological history can also be predicted from thermodynamical arguments [137].
We have obtained these results by applying dynamical system techniques making use of both the linear stability analysis and of center manifold analysis to a Friedman universe filled with three different non-ideal fluids separately in f (R) = R + ζR 2 gravity. We have adopted a set of dimensionless variables proposed in [33] on which we have derived the physical restrictions (33) for preserving the theory from ghost and tachyonic instabilities, obtaining nevertheless a model with a rich variety of cosmological behaviors as previously mentioned. It is also interesting to note that the difference between the curvature energy density and the actual matter content energy density, which can be computed from (20), reads as: Therefore, the two energy densities are equal on the line 2z − y − 2yΩ = 0, which describes a configuration that can actually arise within the physical range (33). Whether this can tame some aspects of the coincidence problem [138] will be explored in future publications, but we should remark that this result has not required to introduce any ad hoc interaction terms between the two fluids by modifying by hands the Bianchi identities unlike in [139][140][141], and therefore we can appreciate already at this stage that this potential solution would not be affected by inconsistent directions of such energy flow. For example as mentioned in [138], the solar system has formed at the epoch in which the abundance of dark energy is of the same order of magnitude of the abundance of regular matter so that a local gravitational collapse can occur in a globally accelerated expanding universe, and in our picture the roles of those two fluids would be played by a gravitational effect and by an actual matter fluid separately. We have as well derived a connection between the dynamical system variables we have adopted and the cosmographic deceleration, jerk and snap parameters. Two equilibria points P 1 and P 2 come with the same values of these cosmographic parameters, and while one of them (P 2 ) admits a well-defined energy density of the cosmic fluid, in the case of the other (P 1 ) it exhibits the indefinite form 0/0. Thus, in future we will investigate whether the same dimensionless variables used here can be connected as well to the positions of the CMBR peaks for removing this ambiguity. We have extended the dynamical system analysis up to infinity by introducing an appropriate compactification of the phase space. As far as the Redlich-Kwong, (Modified) Berthelot, and Dieterici fluids are considered, the region at infinity of the phase space does not carry only an abstract geometrical interpretation, but it corresponds to a regime in which the equation of state for the cosmic fluid reduces to P βρ, as it can be seen from (44). Thermodynamically, this means that the interactions between the fluid constituents are suppressed as it would happen in the limit α → 0. This transition to the ideal behavior of P = w(ρ)ρ fluids has already been met in cosmology [8,9,142], and it has been interpreted as a form of asymptotic freedom analogue to the one which characterizes the quark-gluon plasma [143,144], although in this case is occurring at low rather then high energy densities.
Finally, the dynamical system approach has given us the opportunity of identifying the regions of the phase space which are free from any of the known five finite-time cosmological singularity. In our cosmological models Type II and Type V singularities can occur in the past only in correspondence of the nonsingular bounce at the infinity of the phase space, the latter being a direct consequence of the modifications to the gravity sector. A Type I singularity can occur in the future along the line y = 2z, while a Type III in the past in correspondence of the radiation dominated epochs. Our cosmological models are not affected by a Type IV singularity. Our analysis was completely classical and whether quantum gravity correctionsá la Wheeler-DeWitt affect this picture will be clarified in a future project, as for example done in [145]. Other interesting future projects may consist in analyzing the astrophysical data about recombination epoch, 21-cm line excess at cosmic dawn, and Lyman α forest by exploiting the existence of a radiation-dominated epoch in our models; this can tame the previously mentioned disagreement between the thermodynamical Le Chatelier-Braun principle and the fact that a dark matter epoch should have come before the dark energy one [67,68] since those phenomena are usually addressed via interacting scenarios [146][147][148]. The first attempt of accounting for physical properties of real gases beyond their ideal behavior has been performed by the van der Waals equation of state which implements information about the finite size of the molecules and their mutual interactions assumed to be attractive at large distances and repulsive at short ones via a Lennard-Jones type of potential. Although this proposal came with many desirable features because it can reproduce ideal gas isotherms at high temperature and it exhibits a liquid-gas coexistence phase, the experimental collections of more and more precise data about chemical substances has called for some improved models, as for example the Redlich-Kwong, Berthelot and Dieterici formulations. These models are still based on just two free parameters which are the critical temperature and critical pressure at the coexistence of two phases. Van der Waals' idea of combining the two contributions for the pressure due to the volume occupied by the molecules (which sets a limit on the fluid compressibility), and their internal energy (in the ideal picture molecules only have kinetic energy) simply as an algebraic sum P = P att. + P rep. has been assumed also in the Berthelot and Redlich-Kwong equations of state. They have been proposed as more realistic models for accounting for datasets about the fugacity of hydrocarbons at low (close to the ambient pressure) and high pressure respectively. Intuitively the fugacity quantifies the fleeting properties of a material, while rigorously it is the effective pressure of an ideal gas at the same temperature and with the same molar Gibbs free energy as the real gas; its value for a certain substance is determined from measurements of volume as a function of pressure at constant temperature. The success of the Berthelot and Redlich-Kwong formalisms is grounded in being consistent with experimental data of different substances (methane, ethane, propane, isobutane, etc...) belonging to the family of hydrocarbons just by changing the values of the two free parameters α and β for each of them; before it was necessary to consider a temperature-dependent coefficient in the second-order virial expansion to be empirically reconstructed in each case separately. Thus, this has constituted a great advantage in epochs at which computer simulations were still not widely available. The Redlich-Kwong equation of state has then been further improved by introducing a third parameter known as the acentric factor taking into account non-spherical shapes of the molecules as the Soave-Redlich-Kwong equation for a better description of nonpolar compounds [149]. For a modern treatment of such equations of state we refer to some textbooks as [150,151]. On the other hand, the Dieterici proposal still maintains the idea that two contributions should be included in the pressure (repulsive because molecules are assumed to be hard spheres which cannot penetrate each other, and attractive for having a bound system), but it combines them as P = P rep. e −Patt. improving the agreement with experimental data of the compressibility factor at high pressure than the van der Waals equation [150,151]. In cosmology a similar way of thinking than in chemical thermodynamics has been followed by combining into a single formalism the attractive effects of regular matter and the repulsive one of dark energy: at first the van der Waals equation of state has been chosen for the cosmic fluid [154][155][156][157], and then the Redlich-Kwong, Berthelot and Dieterici ones have been used for enlightening whether those different characteristics which have been observed in a laboratory setting come with specific signatures in cosmology [27].

Appendix B: Stability analysis of finite isolated fixed points
In this Appendix we present in some details the calculations regarding the linear stability analysis for the cosmologically relevant isolated fixed points exhibited in Table I. The stability nature of an isolated fixed point in the linear regime is completely determined by the eigenvalues of the Jacobian matrix evaluated at the fixed point, provided the fixed point is hyperbolic, i.e. none of the eigenvalues is zero. There are four distinct possibilities that may arise for a dynamical system (for the stability classification criteria see for example [88][89][90][91]; for the physical significance of a certain type of stability see instead [158,159]): • If all the eigenvalues have positive real parts, then the fixed point is said to be unstable. An unstable fixed point represents a past attractor in cosmology i.e. an epoch which represents a possible initial state for a cosmological evolution. • If all the eigenvalues have negative real parts, then the fixed point is said to be stable. A stable fixed point represents a future attractor in cosmology, i.e. an epoch which represents a possible final state for a cosmological evolution.
• If two of the eigenvalues are complex conjugate to each other with vanishing real parts, then the fixed point is unstable (stable) whether the third eigenvalue is positive (negative). This represents an oscillatory approach towards the past (future) attractor. The past (future) attractor itself represents an epoch around which the cosmological solution oscillates indefinitely.
If one or more of the eigenvalues of the Jacobian matrix are zero then the fixed point is said to be non-hyperbolic. For non-hyperbolic fixed points Jacobian eigenvalues cannot completely determine the linear stability nature, and center manifold analysis is required to determine the stability of non-hyperbolic fixed points.
In Table V we list the eigenvalues of the Jacobian matrix for the cosmologically relevant isolated fixed points presented in Table I. The eigenvalues are functions of the model parameters, and therefore to determine their signs one must keep in mind that α, ζ > 0, and the existence conditions for the various fixed points from Table  II.

Fixed Points
Redlich-Kwong (Modified) Berthelot Dietrici  (23) and presented in Table I. We remark that the correct physical interpretation of these results require α, ζ > 0, whilst the restrictions on the parameters β can be found in Table II.
For case of Redlich-Kwong fluid with β = − 2 3 , we note however that the fixed point P 7 coincides with P 6 . The fixed point P 6 exists for all values of β, and for β = − 2 3 it is unstable. Therefore one can conclude that the fixed point P 7 with β = − 2 3 for Redlich-Kwong fluid is unstable. Stability analysis in the other cases requires the application of a center manifold analysis. Also we note that for the Redlich-Kwong fluid with β = 1, the fixed points P 2 and P 3 coincide with P 1 , implying that a center manifold analysis for P 1 also allows us to complete the stability analysis of P 2 and P 3 .

Appendix C: Center manifold analysis for P1
Center manifold analysis is significantly mathematically rigorous [95,160], and it has been applied in cosmology in [34,39,[161][162][163][164], just to mention a few examples. We carry out this analysis only for the fixed point P 1 ≡ (2, 1, 0), because the Jacobian at this point has a vanishing eigenvalue irrespective of the model parameters for all the three fluids. Firstly we note that in the cases of Redlich-Kwong and (Modified) Berthelot fluids with β < −1 and in the case of Dietrici fluid with β < − 2 e 2 , P 1 is clearly a saddle and center manifold analysis is not required. In the cases of Redlich-Kwong and (Modified) Berthelot fluids with β = −1 and in the case of Dietrici fluid with β = − 2 e 2 , stability analysis of P 1 requires beyond center manifold analysis than presented here, as two of the eigenvalues vanish, and therefore here we investigate only the cases of Redlich-Kwong and (Modified) Berthelot fluids with β > −1 and Dietrici fluid with β > − 2 e 2 . To perform a center manifold analysis we begin by shifting the fixed point to the origin by applying the coordinate translation In terms of Y, Z, Ω the system (23) becomes The matrix that diagonalizes the Jacobian J(0, 0, 0) is the matrix whose three columns are the three eigenvectors above: We note that there is no linear term in V in any of the equations in the system (C10). This is because by construction the V -axis is along the eigenvector corresponding to the zero eigenvalue. Let us consider the phase trajectories in the neighbourhood of the fixed point P 1 = (0, 0, 0). Considering only the leading contributions at the vicinity of this point, from the system (C10) we can write the following • Redlich-Kwong fluid: • (Modified) Berthelot and Dietrici fluid: Considering the leading order contribution in the vicinity of P 1 ≡ (0, 0, 0), the V −equation from (C10) can be written as withβ = β for Redlich-Kwong, (Modified) Berthelot fluid andβ = e 2 β 2 for Dietrici fluid in this case. Taking one more derivative we get To the leading order approximation, dU d ln|V | , dW d ln|V | are constants depending on β, whose value for different fluids can be calculated by substituting the values of W U from Eq.(C14) into eqs.(C12a),(C12b),(C13a),(C13b). If we define then γ is a β−dependent constant, and the first integral of Eq.(C16) gives It is clear from the above result that irrespective of the sign of γ, evolution of V (τ ) is always away from the origin. The fixed point P 1 is therefore always a saddle.
Appendix D: Stability analysis of invariant submanifolds X i = C is called an invariant submanifold of the dynamical systemẊ = f (X) iḟ Stability of an invariant submanifold is determined by the phase flow in its vicinity. If one considers a point in proximity of the submanifold with a coordinate C + δX i , then the component of the flow normal to the submanifold at that point is determined byδ If ∂fi ∂Xi Xi=C is negative (positive), then the phase flow at that point is towards (away from) the submanifold X i = C, and correspondingly the submanifold is attracting (repelling). If ∂fi ∂Xi Xi=C = 0, further analysis is required. Armed with this concept, we can determine the stability of the invariant submanifolds that arise in our dynamical system: • The submanifold y = 2z can be better specified in the polar coordinate as θ = tan −1 1 2 . From (46) one can compute that Therefore the submanifold θ = tan −1 1 2 is attracting (repelling) for r > √ 5(1 − Ω) (r < √ 5(1 − Ω)). In Cartesian coordinates one can state that the submanifold y = 2z is attracting (repelling) for y 2 + z 2 > 5(1 − Ω) 2 (y 2 + z 2 < 5(1 − Ω) 2 ) respectively. The line r = √ 5(1 − Ω) (y 2 + z 2 = 5(1 − Ω) 2 ) separates the two regions of the submanifold with opposite dynamical characteristics.
Appendix E: Stability analysis of fixed points at infinity The isolated fixed point at infinity P i ≡ 1, tan −1 1 2 , 0 lies at the intersection of three invariant submanifolds, namely Ω = 0, θ = tan −1 1 2 and R = 1. This observation completely determines the stability nature of this fixed point. The submanifold R = 1 is everywhere repelling. The submanifold θ = tan −1 1 2 is attracting at P i (since Ω = 0 and r → ∞ at P i ). The submanifold Ω = 0 is also attracting at P i (since −3y + z = −2y − (y − z) → −∞ at P i , assuming β to be finite). Therefore the fixed point P i is a saddle point in the cases under consideration in this section.