{'fraction_non_alphanumeric': 0.03388153555252627, 'fraction_numerical': 0.012774284332263313, 'mean_word_length': 4.216538593875445, 'pattern_counts': {'":': 0, '<': 14, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A new approach to the analysis of the phase space of f(R)-gravity

We propose a new dynamical system formalism for the analysis of f(R) cosmologies. The new approach eliminates the need for cumbersome inversions to close the dynamical system and allows the analysis of the phase space of f(R)-gravity models which cannot be investigated using the standard technique. Differently form previously proposed similar techniques, the new method is constructed in such a way to associate to the fixed points scale factors, which contain four integration constants (i.e. solutions of fourth order differential equations). In this way a new light is shed on the physical meaning of the fixed points. We apply this technique to some f(R) Lagrangians relevant for inflationary and dark energy models.


Introduction
Since the first formulation of General Relativity (GR), many extensions of the original Einstein equations have been investigated. The reasons of the interest in such theories are quite disparate: from the first attempts to unify geometrically the electromagnetic and gravitational interaction started by Weyl [1], to the understanding of the corrections to the gravitational action typical of quantum field theory in curved spacetime and fundamental unification schemes [2], to the attempts to give a complete geometric explanation of the dark phenomenology.
Among these extensions, the class of theories called f (R)-gravity [3] is the simplest realisation of higher order gravity (order four) and has an important role as a natural model for inflation [4]. More recently, it was shown to have also an interesting (and very debated) role as geometrical dark energy model [8]. In dealing with these theories, the necessity to solve fourth order differential equations is the source of the difficulties in the true understanding of the features of their cosmology. Such problems were the origin of the development of a series of methods to indirectly analyse the physical properties of f (R) cosmologies.
Dynamical System Approach (DSA) has proven to be one of the most effective of such methods. DSA has a number of different realisations [5,6]. In the following we shall consider the one proposed in the 1970's by Collins [7] and successively developed by Wainwright, Ellis and Uggla. This version of DSA is defined in terms of dimensionless, expansion normalised variables and allows a very clear physical interpretation of the results. DSA has had a key role in analysing Bianchi models [9] and minimally coupled scalar tensor [10] in the context of GR and has allowed the exploration of the cosmology of a number of modifications of Einstein's theory [11,12]. In [11], the model f (R) = χR n was analysed with this method for the fist time. Later on the method was extended in [13,14] to the case of a generic f (R). DSA allowed for the first time to analyse in detail the phase space for f (R) cosmologies making some generals statements on these cosmologies and revealing a number of interesting attractor solutions.
In spite of this success the method above has some unsatisfactory aspects. First of all the possibility to analyse a given f (R) theory depends on the exact resolution of an algebraic equation of generic order, or, in the most complicated cases, a transcendental equation. Such operation not only limits the set of possible Lagrangians which can be analysed with DSA, but also introduces a number of singularities, so that in many cases the dynamical system is not of class C 1 .
Another important problem is that the choice of the dynamical system variables was made in such a way to obtain in the fixed points solutions with only two integration constants. This implies that these solutions correspond to the general solutions of the fourth order cosmological model where two integration constants have been set to zero. This is problematic because without knowledge of the full solution in a fixed point it is impossible to characterise the correct behaviour of these cosmologies when these fixed points are nodes.
In addition, since the dynamical system variables are not always independent from each other, the dynamical system might present fixed points which can correspond to inconsistent conditions in the cosmological equations i.e. to have f (R) = 0 and R = 0 for a function f for which f (0) = 0. The presence of these points must be a spurious effect due to the way in which he DSA is constructed.
In this paper we propose a new DSA to deal with the cosmology of f (R) gravity. This method is built in such a way to avoid the necessity of solving exactly algebraic/transcendental equations and to give, in the fixed points, solutions of the cosmological equations which contain four integration constants. We will show that the new formulation contains the results of the old method but reveals unsuspected additional features of the evolution of f (R) cosmologies.
The paper is organised as follow. In section 2 we will give the basic equations. In section 3 we will introduce the original DSA for f (R)-gravity. Section 4 is dedicated to the construction of the new DSA, and in Section 5 the new method is applied to some interesting model of f (R) gravity. Section 6 is dedicated to the conclusions.
Unless otherwise specified, natural units ( = c = k B = 8πG = 1) will be used throughout this paper, Latin indices run from 0 to 3. The symbol ∇ represents the usual covariant derivative and ∂ corresponds to partial differentiation. We use the (−, +, +, +) signature and the Riemann tensor is defined by where the W a bd are the Christoffel symbols (i.e. symmetric in the lower indices), defined by W a bd = 1 2 g ae (g be,d + g ed,b − g bd,e ) .
The Ricci tensor is obtained by contracting the first and the third indices Finally the Hilbert-Einstein action in the presence of matter is given by

Basic Equations
In this paper we are going to deal only with metric f (R) theories, which are characterized by the Action where L m represents the matter contribution. In general the function f is considered an analytic function of the Ricci scalar R and contains a set of additional dimensional parameters indicated by barred Greek letters.
Varying the action with respect to the metric gives the generalisation of the Einstein equations: δg ab represents the stress energy tensor of standard matter. These equations reduce to the standard Einstein field equations when f (R,ᾱ,β...) =ᾱR withᾱ = 1/2.
Our treatment will consider only homogeneous and isotropic spacetimes i.e. Friedmann Lemaître Robertson Walker (FLRW) metrics: where a is the scale factor and k the spatial curvature. We also assume that the cosmic fluid is a prefect fluid with equation of state p = wµ with 0 ≤ w ≤ 1. It is common to write the field equations (2.2) in the metric (2.3) as two equation resembling the Raychaudhuri and the Friedmann equations in GR

5)
H ≡ȧ/a, the prime represents the derivative with respect to R and the "dot" is the derivative with respect to t. The two equations (2.4) are not independent: the second can be obtained deriving the first with respect to the cosmic time t once the Bianchi identities for T m µν is considered. In FLRW these identities take the forṁ which is the same as the GR energy conservation equation for the cosmic fluid.
3 The Original Dynamical Systems Approach for f (R) gravity (in brief ).
Using the equations above we can formulate the cosmic evolution in terms of dynamical systems [13,14] (beware of the differences in signature!). Introducing the general dimensionless variables : where µ m represents the energy density of a perfect fluid that is present in the model. As customary, we also define the logarithmic (dimensionless) "time variable" N = ln a. Note that in choosing this time variable we are assuming that we represent the phase space for H > 0 i.e. we are considering only expanding cosmologies 1 . In this variables the cosmological equations (2.4) are equivalent to the autonomous system: The quantity Υ is defined as Since Υ is a function of R only, the problem of obtaining Υ = Υ(x, y, z, Ω) is reduced to the problem of writing R = R(x, y, z, Ω). This can be achieved from the definitions (3.1): Solving the above equation for R allows one to write R in terms of y and z and close the system (3.2). It is clear that the properties of (3.4) determine the possibility of closing (and therefore analysing) the system (3.2) as well as some of the properties of this system i.e. the differential structure. One example is f (R) = R p exp(qR) for which In this case it is evident that the system is not C(1) as the curve y 2 − pz 2 = 0 is singular. This fact can have serious repercussions on the properties of the flow [14]. 4 The new Dynamical Systems Approach.
We will now start to construct the new approach. Before we define the dynamical system variables we will need, however, two preliminary steps. The first one will concern the form of the action and the second one will be the introduction of new (cosmic) parameters. These (re-)definitions will be the cornerstones of the new method.

The form of the action.
In dealing with dynamical system it is crucial to gain an understanding (and control) over the dimensional structure of the theory we are considering. The reason is that the number of dynamical system variables needed for DSA will also depend on the number of dimensional constants present in the theory. To construct the new method, therefore, we will rewrite the action (2.1) in a special form. In particular we will introduce a constant R 0 such that the product RR 0 is dimensionless. In addition, we will also introduce some dimensionless parameter in the form of Greek letters which will represent the ratio between the coupling constant of the additional invariants in the theory and (a power of) R 0 . In this way the (2.1) can be written as where f has the same properties of the one in (2.1) and R 0 will be assumed non-negative.
The main reason behind the formulation above is that in this way any f (R) action contains only one dimensional constant. Therefore, instead of defining a dynamical variable for each dimensional constant (ᾱ,β, ...), we only need one dynamical variable related to R 0 to analyse the phase space of actions of any complexity. In addition, this setting prevents the appearance of fixed points not consistent with the cosmological equations which are typical of the original DSA. We will use this formulation (4.1) of the action as a starting point in the construction of our new DSA.

New cosmic parameters.
Looking for a different way to constrain the properties of the cosmic fluids, Visser proposed a set of cosmic (or cosmographic) parameters [15][16][17][18] with which one is able to characterize completely a cosmological model. These quantities are directly related with the Taylor development of the scale factor and this property determines also our capability to measure them. With few exceptions [16], in the case of GR only the lowest order cosmic parameter have been so far fully exploited, but in the case of f (R)-gravity the situation is different: higher order parameters become crucial and can be used to characterize the evolution of many important cosmological phenomena e.g. structure formation [19]. However, after a quick look to the f (R) cosmological equations written in terms of these quantities (e.g. [19,20]), one soon realizes that this type of cosmographic parameters are not always the ideal objects to work with. The same happens when one tries to use them to formulate a dynamical system approach: although in principle the (4.2) are the ideal objects to construct the DSA they do not constitute always an advantageous set of variables.
For this reason, it is necessary to look for new sets of parameters which share the structure of the original cosmographic parameters (4.2) but, at the same time, are more suitable to deal with our specific problem. One possibility is to use the Hubble parameter to define the variables:q The variables above have been used in [21] to propose new ways to perform the reconstruction of exact cosmological solutions. Unfortunately also the (4.3) do not prove useful to implement a more powerful dynamical system approach. Another possibility, which will be adopted in the following, is to define This choice, apparently cumbersome, appears much simpler in terms of the logarithmic time N : The (4.5) will be the second cornerstone of the mode we intend to propose. In terms of q, j, s the Ricci scalar and its derivatives read R = 6H 2 s + 4j(q + 1) + (q + 8)q 2 H 2 + 2(2 − q) k a 2 , (4.8) and the cosmological equations (2.4) can be written as respectively. Note that now the cosmological equations are equations for j and s instead of H andḢ. In fact, in (4.9), this last quantity has been substituted according to (4.4).
In the system (4.9) the Ricci scalar only remains present in the function f . This means that the equations above should be supplied with an additional constraint given by the (4.6). This relation will allow to simplify considerably the final dynamical system.

The General Method
The first step to obtain and autonomous system of first order differential equations corresponding to the (4.9) is the definition of the dynamical variables. In the case of single fluid cosmology 2 , we choose the set of variables: (4.10) Note that the variable associated to matter does not coincide exactly with the matter density parameter. This is a manifestation of the non-minimal coupling of matter and gravitation typical of f (R)-gravity. Also, because of our definition of R 0 , the variable A will be always non-negative. In addition to the above variables, we will introduce, like in the original DSA, the logarithmic time variable N . With this choice the cosmological equations (4.9) are completely equivalent to the autonomous system: together with the two constraints The first corresponds to the Hamiltonian (Friedmann) constraint, which guarantees the conservation of matter energy, and the second is simply the definition of the Ricci scalar. The (4.14) These quantities represent the part of the system which depends on the form of the Lagrangian (4.1). Note that, differently from the approaches presented in [13,14], in this version of the Dynamical Systems Approach no resolution of algebraic equation is required to close the system. This means that with the technique proposed here we can in principle analyze all f (R) models. The two constraints allow us to eliminate two variables (J and Q) and reduce the total system to: It is clear that this system posses a minimum of three invariant submanifolds: (i) Ω = 0, (ii) Depending on the form of Y one can also have a fourth invariant submanifold in R = 0. These invariant submanifolds can be imagined as "parts" of the phase space which have the property that any orbit that starts in them is trapped. They can have a very specific physical meaning. For example, the existence of the invariant submanifold K = 0 tells us that if we start in an orbit with zero spatial curvature we cannot evolve towards a positive or negative curvature parts of the phase space and that the existence of Ω = 0 implies that a vacuum cosmology remains vacuum (i.e. standard matter cannot be created or destroyed). The submanifold R = 0 represents, instead, the case in which the Ricci scalar is identically zero. Since the equation for R contains terms with Y at denominator, this submanifold can be singular. Finally, the submanifold A = 0 corresponds to the case in which the constant R 0 is zero and it is of more difficult physical interpretation. It can be thought to correspond to the case in which the gravitational part of the action is identically zero i.e. there is no gravitational interaction 3 . Since A appears also in Y, the submanifold A = 0 can also be singular. The singular character of the invariant submanifolds R = 0 and A = 0, however, does not imply the absence of fixed points. The existence of fixed points that belong to a singular manifold is one of the exotic properties of dynamical systems which are not of order C(1). The coordinates of (and therefore the solution associated with) these points can be explained in terms of the limiting form of f for R, R 0 → 0. In fact, close to the fixed points in the R = 0 and A = 0 manifolds the phase space of a given theory f and the one of its limit for R, R 0 → 0 can be considered isomorphic and will admit the same fixed points. For these fixed points the stability analysis will not necessarily be the standard one. For example, some of the eigenvalues of these points might diverge. We will impose the condition that actual fixed points will exist if the equations at least of class C(2) i.e. that the dynamical system equations, the cosmological equation and the eigenvalues will be finite for a given fixed point.
With the above setting, we are ready to analyse the phase space for a given form of f (R) 4 (and therefore a form of the functions X, Y, Z). The analysis can be done using the classical dynamical system theory i.e. obtaining the fixed points imposing that the Nderivatives of the dynamical variables are zero and using the Hartman-Grobman theorem to determine their stability 5 . Naturally, since the phase space is not compact one will need to perform an asymptotic analysis in order to obtain a full description of the dynamics. Such task, however, will be left for a future work.
The solutions associated to the fixed points can be derived writing the modified Raychaudhuri equation in a fixed point, where the asterisk indicates the value of a variable in a fixed point 6 . In spite of being a third order differential equation the (4.16) is always solvable exactly for H. In addition, since (4.16) is a fourth order equation in a, the solution associated to the fixed points solution contains four integration constants (and up to three different "regimes"). As we will see, in the case in which a theory can be treated with both the new and the original method, there is correspondence between the fixed points. This suggests that the original method implicitly sets to zero some integration constants. We will discover that this can hide information on the actual meaning of the fixed point, particularly when it is a sink or a source. The other key cosmological quantities can be deduced by the cosmological equations once the (4.16) is solved. Specifically one has, from (2.6), In literature one often defines the barotropic factor of the high order corrections considered as an effective fluid. Unfortunately, it will not always be possible to integrate the (4.16) exactly up to an expression for the scale factor. In this case we will rely on considerations on the structure of the equation for a and numerical integrations to obtain the behaviour of the scale factor and the other key quantities. For this reason we will limit ourselves to give the solutions for the scale factor in terms of plots.

Examples
In this section we are going to apply the method described above to some specific models. We will start with two simple ones (f (R) = R n , f (R) = R + αR n ) that are treatable also with the original DSA highlighting common features and differences between these two methods. After that, we will explore other two models (the Starobinsky and the Hu-Sawicki models) that are more physically interesting, but cannot be analysed with the original DSA.

The case of R n -gravity
As a first check for our new formalism let us consider the first model that has been analysed with the original DSA [11]. This model, which is also called sometimes "R n -gravity", is characterized by an action in which the Ricci scalar appears as a generic power rather than linearly and it constitutes the simplest fourth order modification of GR.
Putting the action in the form of Section 4.1 we have Using the original DSA it was realised that for specific values of the parameter n there could be orbits which naturally present a transition between decelerated and accelerated expansion [11]. This model was subsequently subjected to more detailed studies which involved cosmological perturbations [19,[22][23][24] as well as astrophysical and cosmological tests (e.g. [25][26][27]). The result of these investigations and other physical considerations points to the fact that R ngravity is inconsistent with multi-scale observations and should be considered only a toy model. It is known that in terms of the original DSA, the case f ∝ R n is degenerate: two dynamical system variables (y and z) coincide. The special character of R n -gravity is present also in the new approach: since the term R n 0 appear as a factor of R n the equation for A is decoupled and can be excluded. Substituting in the remaining equations the expression for X, Y and Z we obtain The system presents in general three invariant submanifolds: K = 0, Ω = 0 and R = 0. Table  1 contains the standard fixed points of (5.3), together with their associated solution. The stability of the fixed points is shown in the case w = 0 in Table 2.
The solutions associated to the fixed points deserve further discussion. Since now we are solving the full (4.16) these solutions will be specified by a linear differential equation in a which is more complex than the one of the original DSA. For example for the points A, B and C we haveȧ This equation can be solved exactly (but almost always implicitly) only in the case in which two of the constants H i are zero. For H 2 and H 3 zero one has, for example, is clear that the solution has sigmoid behaviour: after a power law growth of the type (5.5), the expansion rate starts, at first, to increase and then to decrease, to approach eventually a constant. This results is confirmed by the numerical check of the first derivative of the solution.
The solutions associated to the points D and E are given by the equatioṅ  which can be integrated exactly to give a pure power law behaviour. Table 1.
It is instructive to compare the results above with the ones of [11]. It is evident that the fixed points we have obtained, correspond to the ones found in [11]. In particular the solutions associated to the fixed points are coincident when H 2 and H 3 are set to zero. This result shows how the original DSA would give at best an incomplete result. For example in (5.4), since the a −1 term is only relevant at small a, if A, B, C are attractors the cosmology Table 2. Stability of the fixed points of f (R) = χR n in the case w = 0. Here A stays for attractor, R for repeller, S for saddle.
will tend to become static after a phase of accelerated expansion. In terms of the stability, instead, the two methods show a complete consistency: the nature of the fixed points is the same as the one in [11]. Specifically, comparing the ranges of stability the possibility of a transition between almost Friedmann (point F) and power law inflation/dark energy era (point G) is present also with the new method. In fact, since both the fixed points are on the K = 0 invariant submanifold, we can check this result explicitly plotting this part of the phase space (see Figure 5). It is evident that there is, in complete agreement with [11], a set of initial condition in which the transition appears. The stability analysis results in Table 2 guarantees that this is the case also for orbits in the full phase space. Let us now consider another important model for fourth order gravity: the one in which a generic power of the Ricci scalar is added to the Hilbert-Einstein term. Historically this is the most studied form of fourth order gravity because of its appearance in several quantum gravity calculations [28]. In cosmology the model with n = 2 has gained particular attention as a model of geometric inflation [4]. The Lagrangian for this class of theories can be written in the form showed in Section 4.1 as f (R 0 R, α) = R 0 R + αR n 0 R n . (5.10) Note that for n < 0 this theory is not defined in R = 0 and in R 0 = 0, but in fact we will see that, due to the presence of derivatives of f with respect to R in the cosmological equations (and therefore in the dynamical system) we will have divergences also in other intervals of n. Substituting the f above in the expression for X, Y and Z we obtain Substituting in the general system (4.15) we obtain (5.12) This system admits four invariant submainfolds (K = 0, Ω = 0, A = 0, R = 0). The last two invariant submanifolds can be singular, but they can still contain some fixed points. As said, the presence of fixed points is based on the requirement of convergence of the cosmological equations and the eigenvalues of the fixed points. This implies that some fixed points will only exist for certain values of n.
The list of the fixed points, their associated solutions and the interval of existence of the fixed points are given in Table 3. Note that the points A − G are the same of the ones found in the case of R n -gravity. The presence of this type of fixed points can be understood thinking that nearby A = 0 the function f can be approximated with its limit for small R 0 . For the values of n for which these fixed points exist the approximate function is R n 0 R n and therefore the fixed points of the case f ∝ R n 0 R n appear also in the phase space of this theory. On top of the points A − G of the previous case we have some additional ones, which we will name H i . These points will exist if α(n − 2) > 0 (remember that the variable A is defined to be non negative) and their number depends on the value of n.
In the H i the scale factor is described by the equatioṅ a a = H 0 + 3H 1 log a + 9H 2 log 2 a, (5.13) which admits the exact solution (5.14) For 4H 0 H 2 − H 2 1 > 0 this solution is monotonically growing with two inflection points at and presents a discontinuity in For t →t + this solution approach to zero, whereas when t →t − the solution presents a vertical asymptote. Table 3. Fixed points of f (R) = R + αR n with their interval of existence and their associated solutions.

Point Coordinates {R, K, Ω, A} Scale Factor Existence
In the case 4H 0 H 2 −H 2 1 < 0 the (5.13) is instead not periodic and approaches a constant. The expanding or contracting character depends on the sign of the quantity P > 0 implies a growing scale factor and P < 0 a decaying one. Plots of the solution (5.13) can be found in Figure 6. In the case 4H 0 H 2 − H 2 1 > 0 the solution (5.14) presents features which are physically very interesting. For times close tot + the solution grows exponentially, then around t * 1 it changes into a decelerated expansion with non-constant deceleration factor. After a coasting phase around t * 2 , the decelerated expansion is followed by a new accelerated expansion phase. In a finite time (at t =t − ), however, the solution becomes singular in the sense that a and its derivatives as well as the Ricci scalar diverges at this specific time. This is a well known property of f (R)-gravity 7 [31,32], but in this context one is able to appreciate both this drawback of the theory and its potential as a model that unify inflation and dark energy.
Note that for H 1,2 = 0 the solution (5.14) reduce to the standard de Sitter solution. This fact on one hand connects the points H i to the de Sitter fixed point of the original DSA. On the other hand gives a hint of the true meaning of de Sitter solutions in the framework of f (R)-gravity.     The stability of the fixed points can be calculated as in the previous case using the Hartman Grobman theorem, and it is illustrated in Table 4. It is tempting to attempt a general derivation for the stability of points H i , however numerical inspection shows that such analysis might be unreliable. The calculations show that when points H i exist there are only two options for their stability: they can be either a saddle or an attractor (or their focus counterpart). In Table 4 we report some samples of the stability of points H i for different values of the parameters α and n. Examples of the phase space for these theories in the form of the invariant submanifold R, A is given in Figures 8. Since this model can be treated also with the original DSA it is useful to make a comparison between the results we obtained above and the ones in [14]. Differently from the case of R n gravity the phase space obtained by the two methods is not the same. In particular, the original DSA returns a phase space with many more fixed points. The origin of this difference is probably to be attributed to the choice of variables of the original DSA. For the common points one can compare the results on the stability and it is easy to verify complete consistency. For example, the de Sitter solution that in [14] is associated with E * has a stability that coincide with the point H of the present analysis.
In the case n = 2 the new DSA returns a phase space with no finite fixed points. It is known that the case n = 2 in the theory (5.10) presents significant physical differences with respect to the other model of this class [28], and it is to be expected that the phase space will reflect these differences. The fact that the phase space does not present a point of type H does not necessarily imply that the model does not have de Sitter solutions (we know in fact that they are present [30]). A careful analysis of the equations shows that the fixed point in this case is asymptotic (A → ∞) and it is therefore excluded by the present analysis.
It is interesting to note that one of the conclusions in [31] is that to avoid the singularity one either has to recur to special initial conditions or to add additional curvature invariants. This result seems consistent with our findings. In order to avoid the singularity one either has to control the initial conditions (by setting to zero some of the constant H i ), the values of the parameters like in Figure 5.2, or hope that adding additional curvature invariants the fixed points H will become irrelevant (in the same way of what happens with the case n = 2 above). Table 4. Stability of the fixed points of f (R) = R + αR n in the case w = 0. Here n 1 is the smallest real solution of the equation 256n 3 − 608n 2 + 417n − 81 = 0 , A stays for attractor, R for repeller, S for saddle, F S for saddle focus. The fixed points appear in the table only if they exist in at least one of the intervals of the parameter n indicated. Table 5. Stability of the fixed points H i of f (R) = R + αR n in the case w = 0. Here A stays for attractor, S for saddle, F A for saddle focus and NA represents the absence of fixed points. In the cases considered there is only one fixed point H.  The issue is that the feature of the inflation dark energy connection seems to be inextricably tied to the approach to the singularity. Thus, in order to generate cosmic acceleration without incurring in the singularity requires that other dynamical mechanisms/forms of the function f have to be found.

The Starobinsky model.
The Starobinsky model is one of the most important class of models of Dark Energy based on fourth order gravitation [35]. The basic idea is to construct, via a function of the Ricci scalar, an effective cosmological constant which is relevant in curved spacetimes, but approaches zero in the case of flat spacetime. The function f of this model is given by with n and λ positive andR 0 of the order of the inverse of the present value of the cosmological constant. The parameter λ and the value of n are related by an algebraic equation. In [35] considerations on the Solar System constraints and cosmological linear perturbation theory were used to find that one should have n ≥ 2 and λ 0.94. Following Section 4.1 we can write the Lagrangian (5.18) for this class of models as where the parameters α, β are related to the ones in (5.18) by the relations λ = α |β| R 0 and The functions X, Y, Z are 20) and the dynamical system (4.15) becomes (5.21) The system presents three invariant submanifolds: K = 0, Ω = 0, A = 0, although the last submanifold is singular. Note that, differently form the previous case, none of the fixed points A − G is present here. This can be explained looking at the limit of the action for R 0 → 0. In this case in fact the action reduces to R 0 R rather than R n 0 R n and therefore the fixed points of R n -gravity are not present. Even if this system is valid for any value of the parameters, it is not possible to find its fixed points analytically 8 : this task would entail the resolution of an algebraic equation of order n for which no general analytical solution is known. We therefore refer to the bounds in [35] and we set from now on n = 2 and β, α = 1. In this case the system admits three fixed points with real coordinates (see Table 6). One of them is on the singular submanifold and presents singular eigenvalues so that it does not appear in Table 6. The other two (H 1 and H 2 ) are both associated to a solution of the type (5.13).
The stability analysis reveals that the character of the points H i is in general different. In our particular case, only one of these point is an attractor whereas the other is unstable. Since the phase space contains invariant submanifolds, we can conclude that for these values of the parameters there is only a specific set of initial conditions which lead to the attractor H 2 . Therefore, also in the case of the Starobinsky model there is the possibility that the cosmology will evolve towards a singularity, but this occurrence depends strictly on the choice of the initial conditions. In orbits which do not approach H, the fate of the cosmological models depends on the presence of asymptotic attractors.
A plot of the section of the invariant submanifold K = 0, Ω = 0 which contains H 1 and H 2 is given in Figure 9.

The Hu-Sawicki model
As a last example we consider a model for geometric DE proposed by Hu and Sawicki [36]. The model was designed to be able to reproduce cosmic acceleration without the explicit introduction of a cosmological constant and, at the same time, to be compatible with cosmological and Solar System tests.

(5.23)
Substituting in (4.15), we obtain the dynamical system This system can present four invariant submanifolds (K = 0, Ω = 0, A = 0 and R = 0), but the last one exists only for 0 < n < 2. The A = 0 submanifold can be singular. This happens for n > 1. Differently from the Starobinsky model, the phase space in this case contains different fixed points depending on the value of the parameter n. In particular, the conditions β > 0 and A > 0 limits strongly the number of fixed point belonging to the class H. The list of fixed points for n = 1 is given in Table 7. The case n = 1 requires a special treatment and will be covered in a separate subsection. As usual, the conditions of existence for the fixed points in Table 7 have been determined asking that the dynamical system should be at least of class C(2) in the fixed point. In one case (the point C) there are values of n for which the Jacobian is divergent, but the eigenvalues converge. We refer as interval of existence of the point C the one of convergence of the eigenvalues. Within the interval of existence of the fixed points the stability can be determined with the standard methods. The results are given in Table 8. As in the case of f (R) = R + αR n the stability of H i needs to be calculated case by case, but it can only be a saddle or an attractor of their focus counterparts. An example of the stability of points H i for different combinations of the values of α, β and n in the case of dust (w = 0) is given in Table 8 In this case the points C and H are the only possible finite attractors for the cosmology. The first point, however, only exists for 0 < n < 1. Therefore for n > 1, we are in the same situation of the Starobinsky model: only a careful choice of the initial conditions could avoid the singularity of solution (5.13).

The case n = 1.
It is interesting to note that the dynamical system (5.24) is not defined for n = 1, which means that the dynamical system in this case needs to be obtained re-deriving X, Y, Z and Table 7. Fixed points of the Hu-Sawicki model for n = 1 and their associated solutions. Here A 0 represents the positive solutions of the equation given in the last line of the Table. Point Coordinates {R, K, Ω, A} Scale Factor Existence the equations (4.15). This procedure yields , (5.25) Point 0 < n < 1/2 1/2 < n < 1 n > 1 and, therefore, In particular, for α = 1 the phase space contains three fixed points of the type H and an additional one H Ω = {R = 2, K = 0, Ω = −(4 + 3w), A = − 1 12β }. However, all of these points have A < 0 for α and β positive and have to be discarded. Therefore, this version of the Hu-Sawicki model presents a finite phase space structure analogous to the one of f (R) = R+αR 2 . The difference is that the Hu-Sawicki model dose not necessarily have asymptotic fixed points of the type H.
For α = 1, the situation is very different. In this case the function f reduces to f (R, R 0 ) = βR 0 R 1 + βR 0 R , (5.27) i.e. the Hilbert Einstein term is eliminated. The phase space for this case contains the points A (which for n = 1 coincides with C) to E. Three fixed points of the type H exist, but two of them have A < 0 for β positive and have to be discarded. The third one has coordinates {R = 2, K = 0, Ω = 0, A = 0}. In addition, two new fixed points with coordinates I = {R = 4, K = 3, Ω = 0, A = 0} and L = {R = 1 4 (5 − 3w), K = 0, Ω = − 9 8 (1 + w), A = 0} are present.
The modified Raychaudhuri equation (4.16) returns the same solutions for the fixed points A-E. The solution associated to I is (5.4) whereas the one associated to L is (5.28) The stability of the fixed points in this case can be calculated in the standard way and with the exception of H all appear to be unstable (see Table 9). Point H has instead one zero eigenvalue and the analysis of its stability requires the use of the Center Manifold Theorem (CMT) [37]. To apply this theorem the first step is to write the system (5.26) in the form where a is a vector that depends on the parameter β. Since the first non-zero term in the previous expression is quadratic, (5.31) implies that point H has the stability of a saddle-node. In Figure 10 we plot an example of the invariant submanifold K = 0, Ω = 0 corresponding to this case.

Conclusions
In this paper we have presented a new approach to analyse the finite phase space of the cosmology of f (R)-gravity. The new method can be applied to any form of the function f  Figure 10. The section of the invariant submanifold K = 0, Ω = 0 in the case n = 1, α = 1, β = 1 12 for the Hu-Sawicki model. Note that on the A > 0 part of this invariant submanifold H appears an attractor, but has a saddle character in the unphysical A < 0 of the phase space (not represented here).
which is analytical in R without the need of cumbersome inversions. The phase space obtained is more regular than the one obtained with the original dynamical systems approach, although it is not possible to eliminate singularities in a complete way as they are due to the essential form of the cosmological equations. The new method also naturally excludes the fixed points which represent states incompatible with the definition of the dynamical system variables. Therefore the new DSA returns a phase space which matches in a closer way the actual evolution of the cosmological equations.
Using the idea of higher order cosmological parameters (like "jolt" and "snap"), the new DSA is able to associate to the fixed points full solutions of the cosmological equations (in the sense of solutions with four integration constants). Among the fixed points found in our analysis, the points labeled H i are surely the most interesting. They correspond to the dominance of the curvature terms and their number depends on the value of key parameters appearing in the function f . The scale factor in H is a transcendental function with the remarkable property to be able to combine an initial exponential expansion a phase of decelerated expansion and a final accelerated expansion phase. Its existence on one hand describes clearly the role of f (R)-gravity as model for double inflation or a unified model for inflation, standard cosmology and dark energy. On the other, however, it confirms clearly that f (R) cosmologies with these properties can run into finite time singularities. This is indicated by the fact that in many of the cases we have analysed, the solution H is the only attractive fixed point of the phase space (although not a global attractor). It is easy to prove with our method, and in accordance to the results in literature (see e.g. [31,33]) that adding a special set of additional invariants one might avoid this scenario, but in general such avoidance requires fine-tuning.
It is important to be careful in considering the nature of the points H. Their associated solution (5.13) contains integration constant which can take any value and it is in general related to the constants appearing in f . This means that the fact that a theory posses an attractor H does not mean that the full behaviour (5.13) is necessarily realised: one could have, for example, that for specific values of the coupling constants and of the H i vanishes. If H 2 and H 3 are zero then H represents a standard de Sitter solution. This depends on the value of the coupling constant of the model as well as the initial condition for the model. However it is evident that the nature of (5.13) has repercussions on the understanding of the actual meaning and role of the de Sitter solutions in the context of f (R)-gravity.
The DSA proposed has first been verified on the simple case of f (R) = R n and we found a complete agreement in terms of fixed points and stability with the original DSA in [11]. However, the solutions associated to the fixed points present significative differences. In fact, even if the solutions obtained with the new method reduce to the ones of the old DSA setting to zero two of the four integration constants, it appears clear that their (time-)asymptotic behaviour is different. This means that, for example, when the fixed points are attractors, the parts of the solution excluded in the original DSA can be dominant. In this respect, therefore, the original DSA returns incomplete information on the cosmology of f (R)-gravity.
As another test of the new DSA, we have considered another form of f that was analysed with the original dynamical system method: f (R) = R + αR n . In this case, it turns out that the only difference between the two treatment is in the number of fixed points: the new method excludes fixed points that give conditions inconsistent with the definition of the variables. The remaining fixed points, however, present a stability that is completely equivalent. The theory f (R) = R + αR n is also the simplest theory that shows points of the type H in the finite phase space. An interesting exception in this respect is the case n = 2 and it is natural to expect that this difference is related somehow to the special properties of this model.
As last steps we have applied the new DSA to two models which were not analysable with the original method, but at the same time constitute important theoretical models for inflation and/or dark energy: the Starobinsky and the Hu-Sawicki models.
The Starobinsky model, unfortunately, cannot be treated in general due to the complexity of the algebraic equations needed to determine the fixed points. We have therefore limited our analysis to the case in which the value of the parameters are chosen to be compatible with Solar System and cosmological perturbations constraints given in [35]. In this specific case, we obtain a phase space in which only two fixed points of the type H appear. Since only one of the points is an attractor we can conclude that the scenario of solution (5.14) is possible in this model, but it is not achievable for general initial conditions.
The case of the Hu-Sawicki model is more involved. Only for specific intervals of the parameter n the phase space admits fixed points different form the type H. The situation is complicated by the fact that the points belong to singular submanifold. The case n = 1 has to be treated separately and reserves a number of surprises. For example, in the case n = 1, α = 1 the phase space has no finite fixed points, much in the same way of the case f (R) = R + αR 2 . The two models, however, differ in the asymptotic structure of the phase space. In the case n = 1, α = 1 we found, instead, that only the (unique) point H is always effectively an attractor.
The new DSA seems therefore to be very efficient in uncovering the features of interesting f (R) cosmologies. However, as all methodologies, the new DSA presents also a series of drawbacks. For example, one would like to be able to consider in an easier way the GR-like states for the cosmology to be able to find (if they exist) "mimicking behaviours" of these theories, not dissimilar to the isotropization mechanism already found in scalar tensor gravity [38,39]. In addition, the impossibility to define a set of compact variables and the fact that our results point clearly to the presence of asymptotic fixed points makes the present analysis incomplete. The resolution of these issues in the context of the new DSA will be the focus of future studies.