Cosmological dynamical systems in modified gravity

The field equations of modified gravity theories, when considering a homogeneous and isotropic cosmological model, always become autonomous differential equations. This relies on the fact that in such models all variables only depend on cosmological time, or another suitably chosen time parameter. Consequently, the field equations can always be cast into the form of a dynamical system, a successful approach to study such models. We propose a perspective that is applicable to many different modified gravity models and relies on the standard cosmological density parameters only, making our choice of variables model independent. The drawback of our approach is a more complicated constraint equation. We demonstrate our procedure studying various modified gravity models and show how much generic information can be extracted before a specific model is considered.


A. Cosmology
Modern cosmology has seen remarkable advances in recent times, and with it, our understanding of the Universe and the gravitational interaction has improved [1,2]. Despite being the weakest of the fundamental forces, it is responsible for the formation and evolution of structures on the largest scales. Hence, cosmology provides a unique testing ground to study the gravitational interaction where it dominates over the other forces.
The best current description of gravity is given by the theory of General Relativity (GR). Einstein's theory has been vastly successful in its predictive and explanatory power [3], though it does face challenges in the dark sector (i.e. relating to dark matter and dark energy). With this in mind, one is encouraged to consider small modifications of General Relativity that look to solve these dark sector problems. In particular, there is no evidence against the presence of additional degrees of freedom that would reproduce the features of dark matter or dark energy, yet cause no impact within the regimes where GR is successful.
Throughout this work we assume the validity of the cosmological principle, namely, that on sufficiently large scales the Universe is homogeneous and isotropic. Consequently, we can model the Universe using the Friedmann-Lemaître-Robertson-Walker (FLRW) line element. Additionally, we restrict our study to spatially flat models, in agreement with most current observational data [4]. This line element reads where a(t) is the scale factor and N (t) is the lapse function. In some theories of modified gravity, it is not necessarily true that the lapse function can be set to one, which is why we keep it arbitrary for now. In fully diffeomorphism invariant theories one can always rescale the time coordinate to absorb the lapse function into the time coordinate. It may be tempting to go beyond the scope of this work and study anisotropic cosmological models like the various Bianchi type models, as these offer interesting theoretical possibilities and contain a much richer structure that the FLRW model. However, anisotropic models are very tightly constrained [2,5] and it is fair to say that there is no observational evidence that challenges the cosmological principle meaningfully.
FLRW cosmology, despite its mathematical simplicity, provides an excellent test bed for modified theories of gravity, see for instance [6][7][8][9][10][11][12][13][14][15][16][17] and references therein. Within the FLRW framework, competing gravitational theories give rise to differing cosmological predictions. Of those predictions, those relating to the current observational tensions are particularly relevant: specifically, the discrepancies around the velocity of receding astronomical bodies (as encoded a Electronic address: c.boehmer@ucl.ac.uk b Electronic address: erik.jensko.19@ucl.ac.uk c Electronic address: ruth.lazkoz@ehu.es in the H 0 value), and the present-day standard deviation of linear matter fluctuations on the scale of 8h −1 (as best enciphered in the S 8 parameter that supersedes the well known σ 8 quantity). However, the current proposals of modified theories that aim to resolve these tensions, either individually or simultaneously, typically are imperfect in the sense they do not reconcile all cosmological parameters. This is precisely the motivation to continue exploring modified gravity models from as many complementary angles as possible.
When considering a homogeneous and isotropic cosmological background, one can take advantage of the ability to cast the field equations of modified gravity theories into those of dynamical systems. It is then possible to explore the evolutionary consequences of the additional degrees of freedom. To this end, we will specify matter-energy content consisting of a matter and a radiation component along with a cosmological constant. This offers the possibility to compare our new scenarios with the consensus ΛCDM model, which is characterised by exhibiting a de Sitter point (cosmological constant dominated) in the asymptotic final state of its evolution, complemented with a radiation dominated repeller (instability along all eigendirections) and a matter dominated saddle point (instability along some eigendirections only). We present an approach of great generality in the definition of variables that will allow us to discuss in detail the differences from the ΛCDM model phase space pattern, and it may serve in the future as a criterion to discard largely discrepant models.

B. Modified gravity
Theories of gravity beyond General Relativity are almost as old as GR itself. Given that GR is formulated on a four dimensional Lorentzian manifold that is torsion free and comes equipped with a metric compatible connection, one is immediately tempted to allow for additional geometrical structures or to increase the number of dimensions. The Einstein-Hilbert action is linear in the curvature scalar, which motivates models with nonlinearities. Matter is coupled minimally in GR, which means it might also be of interest to consider non-minimal couplings. The vast majority of modified gravity models being considered fall somewhere in the above description [18][19][20][21][22][23][24][25][26][27][28][29][30][31][32][33].
For the purpose of the present work, we are particularly interested in second order modified theories of gravity. By this we mean that the gravitational field equations contain at most second derivatives with respect to the independent variables. Those independent variables are the four coordinates, one temporal coordinate, and three spatial coordinates. Given that the overall scale of the Universe should not enter the dynamical equations, other than through a cosmological constant, the cosmological equations of any such modified theory of gravity must take the form Here E 1 and E 2 are two arbitrary functions. This framework is very general and can be used to describe a plethora of second order modified gravity theories. In fact, without introducing non-minimal couplings, all such theories are probably of that form. Next, one divides these two equations by H 2 and introduces the standard density parameters on the right-hand sides to getẼ where we have assumed linear equations of state. This approach can also include scalar fields, in which case the effective equation of state parameter w would not be a constant. If, at this point, the matter satisfies the usual conservation equation (again neglecting non-minimal couplings), a first order equation in the matter and geometrical variables, it is immediately clear that one deals with a first order system of equations. We continue to elaborate on this point in the following subsection.

C. Dynamical systems
A large body of literature already exists which discusses dynamical systems applications in cosmology, see [34] and the many references therein. In the following, we will mainly discuss the choice of variables and the different approaches taken in the past, as this is one of the most relevant topics in the context of the current work. The starting point of the dynamical systems approach has often been the Friedmann constraint equation, which we write in the form where f (1) and f (2) are arbitrary functions depending on the scale factor and its derivatives and ρ other stands for 'matter' terms other than the standard dust or radiation. The cosmological field equations in f (R) gravity, f (T ) gravity, f (Q) gravity and f (G) gravity can all be put into this form. In fact, it is difficult to envisage a model which could not be brought into this form. We note that this is slightly less generic than the above mentioned equations (4) and (5). Equation (6) reduces to General Relativity when f (1) = 1 and f (2) = 0. When working with dynamical systems it is most convenient to introduce dimensionless quantities. This has motivated the following approach when dealing with the Friedmann constraint (6). The dimensions of H 2 and κρ matter are the same, hence f (1) is dimensionless while f (2) also has dimensions of H 2 . Dimensionless variables now appear naturally when equation (6) is divided by a quantity that has the dimensions of H 2 . In standard cosmology with scalar fields, see for instance [35], it is natural to simply divide by H 2 which leads to the cosmological density parameters and the constraint equation takes a particularly simple form. In f (R) gravity, on the other hand, it appears to be natural to divide by the factor 3f (1) H 2 and to introduce dynamical variables that also depend on the chosen function [36]. This leads to a neat set of equations, but can also introduce problems when f (1) vanishes dynamically during the cosmological evolution, the constraint equation again takes a rather simple form. The mentioned problem has motivated the introduction of other variables which remain regular during the cosmological evolution, see in particular [37] and [38] where specific f (R) gravity models were studied, see also [39][40][41]. Often, these improved variables lead to more complicated constraint equations.
Our approach follows a slightly different route. Instead of looking for particularly chosen variables applicable to certain modified gravity models, we simply divide our Friedmann constraint equation by H 2 and introduce the standard cosmological density parameters. The price to pay for keeping this simplicity is an involved constraint equation. The substantial gain from this is that we can apply this method independently of the functions that determine the model. In fact, we can study any model where the function depends on the Hubble function. Moreover, we can present a partial analysis of the model for all functions and extract the generic dynamics for a large class of models. The existence and location of all critical points and critical/singular lines only depends on the values of the function and its first and second derivatives at a finite set of points. Often, we can determine the stability properties of the critical points without having to specify the model altogether. This happens when the eigenvalues give definite values which are independent of the model.
Irrespective of the specific function chosen, we find that the de Sitter point (when it exists) is a stable late time attractor. This is a surprisingly general result which holds in many different modified theories of gravity. It shows that a late time accelerating cosmological solution is a truly generic property of the cosmological models we study.
The principal ideas of our approach are similar to those used in [42], applied to f (T ) gravity, and predating some of the more recent work on other second order modified gravity theories where their approach can also be used. The authors introduce what they call the Friedmann function and rewrite the dynamical equations using this function and its derivative. Their approach uses a dimensionless variable related to the ratio of matter and radiation. However, to recover the standard cosmological matter and radiation densities, the Friedmann function is needed, meaning that it is somewhat more difficult to obtain model independent statements about the matter content at certain critical points.
D. Modified gravity models -f (G), f (T ) and f (Q) In the following, we deal with three classes of modified gravity theories which are generally known as f (G), f (T ) and f (Q) gravity, all of which can be found in [31,32,[43][44][45] and the previously mentioned overview papers. All these models have in common that the cosmological field equations are compatible with the choice N = 1 for the lapse function. Moreover, the cosmological equations are also of second order only. All models are based on a single scalar (or pseudo-scalar in the case of f (G) gravity), which is indirectly linked to the standard curvature scalar.
When considering spatially flat FLRW models, it turns out that the relevant scalars defining the respective theories are given by Perhaps unsurprisingly, the cosmological field equations implied by these different models all coincide provided the signs are taken care of and some minor re-definitions are made to compensate for unimportant factors. For all three classes of modified gravity models the cosmological field equations can be written as The reason for this equivalence can be understood by carefully studying the various terms which connect these theories which are all defined in different geometrical settings. The relevant terms were identified in [31] where equivalence for more general spacetimes was also discussed. For all three models, the energy-momentum conservation equation is automatically satisfied and does not yield additional conditions to be satisfied. This is due to either the choice of suitable tetrads or coordinates, depending on the model in question. A cosmological constant Λ can always be included in these models by using f → f + 2Λ. Here f is the function specifying the model and depending on the geometric setting. Likewise, F ′ stands for the second derivative with respect to the argument. Note that the field equations of f (R) gravity, for example, cannot be brought into this form as this is a fourth order theory while all the above are second order theories.
The following discussion will be valid for all of the three mentioned models. However, in order not to repeat this fact many times and for ease of notation, we will sometimes use Q = 6H 2 and will often work with Q. All of the following could be formulated equally using either T or G.

A. Towards a dynamical systems formulation
The identification of suitable variables for a dynamical systems formulation is generally challenging for modified gravity models, see [34]. One can often introduce distinct, well-motivated variables which lead to distinct dynamical systems with different features [36][37][38]. Following one of the standard approaches, we begin with the Friedmann equation (8) and divide by 6F H 2 which gives where we assume F = 0. This only excludes constant functions which are of no interest to us. Since the signs of both f and F are undetermined, in general, we begin by introducing the two variables Should the model under consideration contain multiple matter sources like matter and radiation, one would naturally introduce several different variables X i , one for each matter source.
The two variables X and Y are not sufficient to rewrite the field equations in the form of a dynamical system as one is not able to eliminate the Hubble function from the equations. We address this issue by introducing a third variable, related to the Hubble function. As Q is related to the Hubble function via Q = 6H 2 , we have Q ≥ 0 at all times. The new variable is and satisfies 0 ≤ z ≤ 1 for all times, which corresponds to 0 ≤ Q < ∞. Clearly one can always introduce such a variable irrespective of the modified gravity model that is being considered. Here H 0 is a constant Hubble parameter and Q 0 = 6H 2 0 . One would then expect to have two independent variables, either {X, z} or {Y, z}, for one matter source or N + 1 independent variables for N matter (fluid with given equation of state) sources.
However, as any function of H can be rewritten in terms of z, the variable Y is not strictly independent. We shall see in the following section that the system can be reduced by eliminating Y altogether, despite its apparent necessity in (12). Here we shall allow Y to represent the free parameters that may occur in f (e.g. a cosmological constant term) but are not captured in z. The Friedmann constraint (12) then determines the values of the free parameters at each point of the phase space. If there are no free parameters, then (12) singles out the allowed 1D trajectory. This will be elucidated in the concrete examples below.
The introduction of a scalar field would require the introduction of two further variables, one related to its kinetic energy and one for the potential. The dynamical equations for X and z are given by where N = log a. The functions f, g,f ,g are defined through the above equations. The function m(z) is defined by taking into account the relation between z and H as given in (13). One can now make various statement about this system without specifying the function f or equivalently the function m(z). Let us remark that the above system is not well-defined if f = c 1 √ 6H 2 + c 2 as this gives m = −1. In this case, the original equations should be considered separately for this particular choice. We will not discuss this particular case separately as solutions of this type were already addressed in [43], in the context of f (Q) gravity.
The line X = 0 needs to be treated carefully as both dynamical equations vanish along this line. The location of all possible critical points is established as follows. Excluding X = 0 the first equation can vanish when X = (1 + m(z))/(2 + m(z)), assuming 2 + m(z) = 0. This means, any critical point will have this X-coordinate for any value of z that makes the right-hand side of the z-equation vanish. Clearly, we have the two values z = 0 and z = 1 provided that 2 + m(0) = 0 and 2 + m(1) = 0. Moreover, any value z c for which m(z c ) = −1 corresponds to a singular line as the system would be ill-defined. On the other hand, if z ⋆ is such that m(z ⋆ ) → ∞, we would also find critical points. Let us summarise these general finding in Table I as follows: Let us note that for some particular functions, things could be rather contrived. For instance, a function for which m(1) → ∞, would mean that B and C would not exist independently. Other particular scenarios could be set up. When discussing specific models, we highlight these issues explicitly. Lastly, while the range of the variable z is bounded from below and above, things are less clear for X. If F cannot change sign, we would have 0 ≤ X, however, there would be no obvious upper bound. This means that the generic phase space is a strip of height one and possibly infinite width. Certain models will reduce the size of this strip. We will later use a poor mans version of compactification where we introduceX = X/(1 + X) which is bounded from above by 1 as our X will be positive, see also [42].
Let us also note for future reference that the X-nullclines are given by the curve X = (1 + m(z))/(2 + m(z)) while the z-nullclines are located along the lines z = 0 and z = 1.

B. Stability discussion
The critical points of system (14)-(15), as given in Table I, depend on the explicit form of m(z). Nonetheless it is possible to make quite general statements about the nature of the critical points. The stability matrix or Jacobian is given by Let us denote the two eigenvalues of J by λ 1 and λ 2 . A direct calculation shows that for X = 0, which is the line L1, the eigenvalues are (0, 3) for every point on this line. Therefore, the system is not hyperbolic along L1, this is perhaps expected for a critical line. The trajectories or orbits of the system satisfy dX/dz = f /g =f /g, which means these are unaffected when X = 0. However, the direction of the trajectories with respect to time does change, they reverse for X < 0 but do not change for X > 0. This simple observation is sufficient to determine whether the line X = 0 attracts or repels trajectories. The properties of the line L2 are more difficult to establish as they depend on m ′ (z), the first derivative of m with respect to z.
For point A, one directly finds λ 1 = 3 and λ 2 = −6/(2 + m(0)). Hence, point A is never stable. Depending on the value of m(0) this point is either a saddle point or a repeller.
Point B is challenging to study as z = 1 corresponds to 6H 2 → ∞ which can pose problems in the definition of m, see (16). One can show that λ 1 = 3(1+m(1))/(2+m (1)), however the second eigenvalue can diverge in the limit z → 1 if m ′ (1) does not decrease fast enough. Likewise, for point C one can establish λ 1 = 3 but the second eigenvalues will again depend on the derivative of m. This means that point C cannot be stable as one eigenvalue is always positive.
For the parameter range −2 < m(1) < −1 the stability will depend on the other eigenvalue provided it is well defined.
We will look at some examples which discuss the above features, almost all of which are absent when studying the simple case of General Relativity.

C. General Relativity
It makes sense to briefly revisit the simple case of General Relativity, where f (6H 2 ) = 6H 2 + 2λ 0 so that F = 1 and F ′ = 0 which implies m = 0. Here λ 0 stands for a cosmological constant term. This means the line L2 does not exist and the point C also does not exist either. The line L1 is present and we have two critical points A and B. The variable X now becomes X = ρ/(6H 2 ) = Ω m /2 where Ω m is the standard matter density parameter, hence From the Friedmann constraint X + Y = 1, one has that each X value fixes Y . And as mentioned previously, the Y variable can be rewritten in terms of z plus free parameters. Therefore each point {X, z} uniquely determines the value of the free parameters in f . In other words, each trajectory represents the evolution of ρ m and H for a different value of λ 0 .
We can now interpret the two critical points. At point A when z = 0 we have H = 0, which immediately gives us a static solution. The finite value of X = 1/2 implies that ρ must vanish at the same rate as H 2 so that we have a static vacuum solution. This solution is only consistent with λ 0 = 0 and corresponds to Minkowski space.
This leads to the following conclusions: the matter-dominated universe is the early time repeller from which all trajectories start. All trajectories are then attracted towards the line X = 0 where they terminate at some value z. If we denote one such value by z 0 then we will have a corresponding H 0 . When X = 0 we have Y = 1 which means which gives H 0 = λ 0 /3. Consequently, the line X = 0 corresponds to de Sitter type solutions parameterised by different values of H 0 . In this formulation of General Relativity with a cosmological term we find an early time universe that is matter dominated, which then evolves towards a dark energy dominated universe. This is visualised in Fig. 1. One can state various explicit results due to the simplicity of the equations when m = 0. One can immediately integrate the X-equation (14) which gives where c 1 is a constant of integration. In line with the above discussion we have lim N →−∞ X = 1/2 and lim N →∞ X = 0 confirming that all trajectories begin at X = 1/2 and terminate at X = 0. Dividing both equations (14) and (15) yields which can easily be solved using separation of variables. This leads to Consequently, all trajectories starting at X = 1/2 satisfy z = 1, as expected, which corresponds to point B. When X = 0 we have z = 1/(1 + c 2 ) which determines vertical position along the X = 0 axis. As stated earlier, this fixes the Hubble constant and corresponds to a de Sitter type solution.

D. Dimensional reduction
The previous discussion follows the 'standard' approach of writing the cosmological field equations in the form of a dynamical system. This standard approach is based on the idea of rewriting the Friedmann equation (8) as a linear combination of various terms which subsequently serve as the dynamical variables of the system. Above, we had the two variables X and Y and the simple constraint X + Y = 1. As it was not possible to eliminate the Hubble function from the dynamical equations, it became necessary to introduce a third variable z related to the Hubble function. This led to the two-dimensional system which was studied, with the second dimension representing the freedom in the parameters of f .
As we noted, the variables Y and z are not completely independent and it is possible to reduce the dimensionality of the system further, see also [42]. To see this, let us first express Y entirely in terms of Q which gives On the right-hand side we make it explicit that Y can be seen as a function of Q. Using the chain rule we can writė while equation (13) yieldsż Combining these previous two equations and the constraint −Ẋ =Ẏ we arrive at It is important to note that Q can be re-written in terms of z. Consequently, one can now replaceẊ on the left-hand side of (14) using equation (25). Next replace X by Y , which is again a function of z, so that one arrives at a single equation in z. This slightly convoluted derivation gives the final equation where we introduced the new function which mirrors the form of the previous variable Y . Clearly, the entire dynamics for any given f (Q) is now determined by the single first order ODE (26). For General Relativity with a cosmological constant f (Q) = Q + 2λ 0 one finds m = 0, as stated earlier, and h = 3z/(1 − z) + λ 0 /H 2 0 . The resulting first-order ODE becomes It is straightforward to show that the dynamics encoded in (28) matches the previous discussion in two dimensions. Similar to before, the relation H ⋆ = λ 0 /3 emerges immediately from the stationary point of this equation. Here H ⋆ stands for the value of H where the right-hand side of (28) vanishes. This value may differ from H 0 which will determine the corresponding value of z ⋆ . It is, however, clear from (28) that the ratio λ 0 /H 2 0 is the only free parameter in this ODE.
Let us explain in different words why this dimensional reduction is possible. In the preceding sections there were three dynamical variables, which roughly related to ρ, H and f . Together with the Friedmann constraint, one studies a system in two variables, {H, ρ} say. However, since f is a function of H or z, the system is in fact 'over-determined'. This is clearer when considering GR without a cosmological constant, f (Q) = Q. From the Friedmann constraint one immediately has X = 1/2, and we're really dealing with a one-dimensional system (line), not a two-dimensional system. All other points on the two-dimensional phase space would be unphysical. Fig. 1 can be seen as a collection of phase lines (i.e. for each individual trajectory) for different values of the fixed parameter λ 0 . In the remainder of this section and the succeeding ones, we remove Y as dynamical variable and study the systems for fixed values of the constants in f .

E. f (Q) gravity example -Beltrán et al. model
Let us apply this approach to the following model which was suggested in [45] with partial results regarding a dynamical systems analysis. We also included a cosmological constant term λ and introduced factors of H 0 such that both parameters λ p and λ become dimensionless, which simplifies the subsequent discussion. Using the above formalism we developed, it is straightforward to compute explicit expressions for m(z) and h(z). When these are put into the first order ODE (26) one finds For λ p < 0 the denominator cannot vanish and we are dealing with a regular ODE. When setting λ p = 0 we are back at GR with a cosmological constant as discussed above. When λ p > 0 we find that the denominator of (30) can become zero when We note that 0 < z s ≤ 1 for all λ p > 0, which means that for positive λ p one always finds a singular denominator. Other than z = 0 and z = 1, the numerator can also vanish, giving rise to up to two critical points. These are which can give two distinct critical points in the permissible range if λ p > 0 and λ < − 3λ p . A direct calculation shows that b ′ (z c ) = 3, which means that any critical point z c will always be unstable for any choice of parameter. The typical phase space for models with λ p < 0 is shown in Fig. 2 for varying value of the cosmological constant λ. All solutions evolve either towards z = 0 or z = 1, depending on the initial conditions. The following Fig. 3 shows the typical phase space plot for λ p > 0. The dashed (red) vertical line represents the singular line of the ODE. Since this line corresponds to a fixed value of z = z s , we have that the Hubble parameter also approaches a constant values H(z) → H(z s ) = H s . However, we note that lim z→zs H ′ (z) → ±∞, depending on the initial conditions. In the following we will consider models with matter and radiation. These models will be two-dimensional and will be conveniently studied using dynamical systems techniques. This will be analogous to the approach taken in GR, where models with matter and radiation can also be studied using a two-dimensional phase space.

A. General setup with matter and radiation
Generalising the formulation in Section II D, where the dimensionality of the system is one fewer than in the previous discussions, we now look at the case with two fluids, matter ρ m (w = 0) and radiation ρ r (w = 1/3). The cosmological equations can be written as Let us introduce the following dynamical variables, which we note differ from (12) by having no F in the denominator of the fluid variables The first two variables are the canonical EN-variables commonly used in cosmology, and Z is the same as z before. Note that none of the variables (35) depend on the choice of function f . This is possible because any function of H can be written in terms the variable Z, which will only enter the dynamics via the Friedmann constraint itself. The benefit is that X 1 and X 2 are now non-negative and take a familiar form. The Friedmann constraint can then be rewritten solely in terms of the new variables where in the final line we have treated f = f (6H 2 0 Z/(1 − Z)) as a function of Z, and similarly for F . These functions can always be seen in two different ways: either as functions of the scalar that determines the modified theory of gravity; or as functions of the Hubble parameter, which we express in terms of Z. Writing f and F in this way means that the expression for X 1 + X 2 is completely in terms of Z once f is specified. We can therefore use this equation to eliminate either X 1 or X 2 , making the phase space two-dimensional. Let us choose to eliminate X 1 and consider the evolution of {X 2 , Z}, where we have made use of (36). Note that f has dimensions H 2 0 , F is always dimensionless and F ′ has units of H −2 0 . Consequently, equations (37)- (38) are dimensionless, with the numerators and denominators all having equal powers of H 0 .
Looking at equations (37)- (38), one can make some general statements about the systems. For convenience, we introduce the functions such that the dynamical equations can be compactly rewritten as Again we remind the reader that here f is defined as f = f ( Note that once f is specified, we are dealing with a closed system of equations which can be studied. However, it turns out that some properties of this system are independent of this function and hence are valid for all models, making our approach very broad.

B. Fixed Points
The fixed points of the system can be divided into two families of solutions, one with X 2 = 0 and the other with X 2 = (4m + 3n − 2)F =: X * 2 , provided they exist and are within the physical phase space (see discussion below). For the case where X 2 = 0, the points with Z = Z * are solutions of the algebraic equation For the X 2 = X * 2 case, the Z coordinates are Z = 0 or Z = 1. The X 2 coordinates must then be evaluated at Z = 0 or Z = 1 because X * 2 is function of Z. These results have been collected in Table II with the P's representing points along the X 2 = 0 line with Z = Z * satisfying (43). The special cases of Z * when n(Z) = 2 and m(Z) → ∞ are labelled as P n and P m respectively. Note that m(Z) can diverge when F → 0, and F appears in both (41) and (42) and the denominator of n(Z) (40). However, F → 0 alone does not represent a critical point of the system. For this reason, we also require n(Z) to be finite whilst m(Z) → ∞. In other words, we require that F ′ → ∞ or that both f → 0 and F → 0 for some value(s) of Z.
The points A and B are solutions with X 2 = X * 2 evaluated at Z = 0 and Z = 1. Lastly, just like in the previous analysis, there exists a singular line L1 if there exists values Z = Z c such that m(Z c ) = −1. There also exists a singular line L2 when n(Z) → ∞ with m(Z) finite. This means that either f → ∞ or that both F → 0 and F ′ → 0 such that n(Z) diverges but m(Z) does not, similar to the conditions for P m .
The corresponding values of X 1 determined by the Friedmann constraint have also been included. The final column states the specific conditions for each point, and it is assumed that all variables X 1 , X 2 , Z must be finite for existence. Note that the fixed points are not necessarily unique and in many cases they coincide with one another, see Table II.  (42) for arbitrary function f .

C. Physical Parameters
Using the formulation above (39)-(42) along with the cosmological field equations, the deceleration parameter can be expressed in terms of the new variables The value of q for each of the fixed points is given in Table III. For the points P n , P m , A and B one obtains fixed values of q, which are model independent. For the lines L1 and L2, the deceleration parameter diverges. The points P 1 and P 2 depend on the functions m and n evaluated at their respective Z values.
The effective equation of state parameter w eff can be expressed in a similar form where the total energy density is ρ tot = ρ m + ρ r + ρ f and ρ f represents the extra terms in (33) that do not appear in the standard Friedmann equation in General Relativity. Similarly, the total pressure is p tot = p m + p r + p f with p f representing the additional terms in the acceleration equation (34). The values of w eff evaluated at the critical points are also given in Table III. Lastly, the values of the Hubble parameter H(t) are also included, provided they are well-defined. Note again that we require X 1 and X 2 to be non-negative and that 0 ≤ Z ≤ 1. In order to enforce these conditions on the physical phase space, we must use equation (36) for a given f to constrain X 2 and Z. Then one can simply read off the fixed points from Table II (that satisfy their appropriate existence conditions) and determine whether they are located within the physical phase space. For example, for f (Q) = Q equation (36) simplifies to X 1 + X 2 = 1 and the physical phase space in {X 2 , Z} is the unit square. On the other hand, for General Relativity with a cosmological constant f (Q) = Q + 6ΛH 2 0 and Λ > 0 one instead has X 1 + X 2 = 1 + Λ − Λ/Z. The physical phase space is then the region satisfying Z ≥ Λ/(1 − X 2 + Λ) and 0 ≤ Z ≤ 1, X 2 ≤ 1, see Fig. 4. Only the fixed points in Table II that lie within the region constrained by (36) are physical. In Fig. 4 we see the points B, P 2 and P n , representing the radiation repeller, matter saddle and de Sitter attractor, respectively.
Before studying other concrete models in this formulation, we quickly look at the linear stability of the system in full generality.

D. Stability Analysis
As in the previous section, the Jacobian matrix or stability matrix can be analysed at the fixed points in Table II to make some general statements about the stability of critical points. The necessary existence conditions are also assumed, but care must be taken with the derivative of the function m(Z), which will contain third derivatives of f (the derivatives of n(Z) can be rewritten in terms of m and n). It is often necessary to assume that F ′′ , which appears in the numerator of the elements of the Jacobian matrix, remains finite at the fixed point in order to determine the stability for a general f . For P 1 = (0, 0) one obtains λ 1 = (2 − 4m(0) − 3n(0))/(1 + m(0)) and λ 2 = (3n(0) − 6)/(1 + m(0)), which has at least one negative eigenvalue. Conditions on m(Z) and n(Z) could then be imposed to determine its stability. The second point P 2 = (0, 1) has the same first eigenvalue λ 1 = (2 − 4m(1) − 3n(1))/(1 + m(1)) and its second eigenvalue is λ 2 = −(3n(1) − 6)/(1 + m(1)), where the close similarities should be noted but with m and n evaluated at Z = 1. Unsurprisingly, with m and n often undefined at these points, there is little to be said about the stability of these points in general, see again their definitions (39) and (40).
For the point P n = (0, Z * ) with n(Z * ) = 2 one gets constant eigenvalues of λ 1 = −4 and λ 2 = −3, meaning this point is always stable when it exists. This is consistent with its interpretation as being a late-time de Sitter attractor with q = −1. The stability of the point P m = (0, Z * ) with m(Z * ) → ∞ depends on the third derivatives of f , so one cannot make definitive claims about its stability. However, if one checks that both n(Z * ) and F ′′ (Z * ) stay finite whilst m(Z * ) → ∞, the eigenvalues will also be λ 1 = −4 and λ 2 = −3, indicating stability. Point A = (X * 2 , 0), with X * 2 defined previously, has eigenvalues λ 1 = −4, λ 2 = −(2 − 4m(0) − 3n(0))/(1 + m(0)). The second eigenvalue is just minus the first of point P 1 's, therefore both points cannot be stable. For point B the eigenvalues depend on the higher derivatives of f , so again, the stability cannot be determined without specifying a model f . If all derivatives of f stay finite at Z = 1 it is possible to at least conclude that this point must be unstable, as the eigenvalues cannot both be negative.
For the lines L1 and L2, all components of the Jacobian diverge, making an analysis is not possible without defining a specific f (Q). Note that in this analysis we have again assumed that all quantities are finite unless otherwise stated, such that none of the components of the stability matrix diverge at the fixed points considered above. The model we are considering now is given by and we assume the presence of two fluids, matter and radiation. In this case, the dynamical equations become and the functions m(Z) and n(Z) are given by We will assume that λ p = 0 as this would reduce the model back to GR with a cosmological constant term. The sign of the parameter λ p determines whether the phase space is naturally compact or not, so we will look at these two cases separately.
From Table II the existence of generic fixed points can be examined, but the Friedmann constraint (36) will need to be used to determine which are physically relevant. Both conditions for P 1 and P 2 are satisfied. For the point P n with n(Z) = 2, the condition is that either λ p < 0 or λ p ≥ 0 with λ 2 − 3λ p ≥ 0. In fact, one can also see that for λ p < 0 the algebraic equation n(Z) = 2 has one solution for 0 ≤ Z ≤ 1, except in the special case when λ is chosen such that n(Z) does not diverge for any Z in the physical range (see Appendix A for details, where in this pathological case the point P n is replaced by P m ). For λ p > 0 the equation n(Z) = 2 can either have zero, one, or two solutions within 0 ≤ Z ≤ 1 depending on whether λ is greater than, equal to, or less than − 3λ p .
The point P m does not exist -except in the special case mentioned above which is discussed in Appendix A -because m(Z) and n(Z) have the same denominators, meaning they both diverge for the same values of Z. Evaluating X * 2 at Z = 0 and Z = 1 reveals that point A has X 2 → ±∞, whilst point B has X 2 = 1. Lastly, the line L1 with m(Z) = −1 exists for λ p > 0 with at most one solution in the physical Z range, whereas L2 does not exist because n(Z) cannot diverge with m(Z) staying finite. Next we will look at the phase space plots for this model.
First, assume λ p < 0 and λ ≤ 0. In this case, the model is compact and X 1 and X 2 are bounded between [0, 1], which can be deduced from (36). There is also the additional constraint from (36) for the physical region to satisfy the inequality which must hold for any values of the parameters. The phase portraits for λ < 0 and λ = 0 are given in Fig. 5a and Fig. 5b. The fixed points B, P 2 and P n represent the radiation dominated repeller Ω r = 1, matter dominated saddle point Ω m = 1 and de Sitter attractor Ω Λ = 1, respectively. This can be deduced from the the definitions of X 2 and Z at the fixed points, and by using the Friedmann constraint to determine the value of X 1 (Ω m ). From Table III the effective equation of state parameter can be read off at the points B and P n as w eff = 1/3 and w eff = −1. For the matter saddle point P 2 one easily finds w eff = 0 as expected. All orbits start at B and end at P n , with some trajectories attracted towards the matter saddle P 2 . The qualitative features are very similar to those of GR with a positive cosmological constant, shown in Fig. 4. When λ p is negative but λ > 0 the phase space is still compact but no longer bounded between [0, 1] in the fluid variables X 1 and X 2 . The late time de Sitter point P n at (0, Z * ) remains the only global late-time attractor of the physical phase space, and the other two fixed points are the same. However, trajectories beginning at the pastattractor, point B, can instead follow trajectories in the the positive X 2 direction before terminating at the future-time attractor P n , shown in Fig. 5c. This is because equation (51), which defines the physical phase space, changes when the parameters of the model change.
For the case when λ p > 0 the situation differs. Firstly, the variables X 1 and X 2 are no longer bounded from above, meaning the phase space is not compact. It is instead constrained by the inequality along with 0 ≤ Z ≤ 1, as usual. The phase space is then infinite in the positive X 2 direction as Z tends to zero. Analogous to the 1D case (Fig. 3), the denominator of equations (47) and (48) vanishes when the Z coordinate takes the fixed value Z = Z c (λ p ) such that m(Z c ) = −1. This corresponds to the line L1 with Z coordinate Z c = (−2 3λ p + λ p )/(λ p − 12). If λ > − 3λ p the physical region is connected, whereas for λ ≤ − 3λ p the space bifurcates into two disconnected parts at the point point {0, Z = Z c }. However, for all λ, trajectories are always confined to regions either above or below the line L1 at Z = Z c . This will be shown clearly in the compactified phase portraits, Fig. 6.
In order to compactify the phase space we introduce the variablẽ such that 0 ≤X 2 ≤ 1. This is what we referred to earlier as our poor man's version of compactification. The new dynamical system in {X 2 , Z} can then be found by taking the derivative ofX 2 with respect to N and using equations (47)-(48). Similarly, the Friedmann constraint is rewritten withX 2 leading to the compactified phase spaces given in Fig. 6. Note that the fixed points in Table II and Table III are still valid in the original variables X 2 and Z, therefore they will also represent critical points in the compactified variables (at the newX 2 coordinates). The two critical points B and P 2 are the same as in the λ p < 0 case, being the unstable radiation dominated repeller Ω r = 1 and the matter dominated saddle Ω m = 1 respectively. The line L1 splits the phase space, as shown by Figs. 6a-6c, and it is quite manifest how the trajectories change direction near this line. This is due to a sign At pointÃ one finds that X 2 = Ω r → ∞ and X 1 = Ω m → ∞, whilst at P 1 we have X 2 = Ω r = 0 and X 1 = Ω m → ∞. At P 1 the effective equation of state parameter is w eff = −2, representing a phantom regime. The new pointÃ can be seen as an early time repeller as all trajectories in the bottom section of the phase space originate there. This is best interpreted as a Big Bang like state where both matter sources diverge. This is consistent with the standard GR framework where ρ m ∝ a −3 and ρ r ∝ a −4 which both diverge as a → 0. Figs. 6a and 6b show the phase space for λ > − 3λ p , with the phase space beginning to pinch off in Fig. 6b as λ decreases. The pointÃ atX 2 = 1 corresponds to the (asymptotic) point with X 2 → ∞. The points P 2 and B are all still present, and P 1 is now within the physical phase space, as mentioned above. The dashed (red) line L1 separates the phase space into two disconnected parts, with trajectories approaching the line from either side. Here the analysis is largely the same as the 1D case in the previous section, with all trajectories attracted to the singular line L1. The relative magnitudes of the parameters determine the shape of the phase space. For the case with λ < − 3λ p , the phase space splits into two disjoint parts along the line L1. Fig. 6c shows the limiting case with λ = − 3λ p . The late-time de Sitter points P n exist for λ ≤ − 3λ p and are solutions to the equation n(Z) = Z 2 (36 + 12λ − λ p ) − λ p + 2Z(λ p − 6λ) λ p − 2λ p Z + Z 2 (36 + λ p ) = 2 , which has two solutions in the physical Z-range when the inequality is strict, see Fig. 6d. In the particular case λ = − 3λ p , the two points merge into one, which coincides with the line L1 (Fig. 6c). The location of the 'first' P n is on this singular line, whilst for smaller values of λ the 'second' point P n appears. As mentioned in the stability analysis, the points P n are always stable, and they attract all trajectories in their respective sections of the phase space. As λ decreases relative to λ p , the P n 's move apart along the Z-axis. This is an interesting bifurcation rarely discussed before in cosmological models.

IV. CONCLUSIONS
We looked at cosmological dynamical systems arising in different modified theories of gravity. By focusing on theories with second order field equations -for example, those in f (G), f (T ) and f (Q) gravity -a model independent approach could be introduced. This is accomplished by introducing the standard cosmological matter variables along with a dynamical variable directly related to the Hubble function H, which itself is directly related to the geometric scalars G, T and Q which define the different theories. This is a unique approach that has not been studied before and is only applicable to cosmological equations that are of second order in the derivatives of the metric. For example, in f (R) gravity the Ricci scalar includes terms proportional to both H andḢ, meaning the variables chosen here are insufficient to close the system. In our case, the system will always be closed once a function f is specified because any function of G, T or Q can be written as a function of H. This can also include non-minimal couplings of matter provided the equations remain second order.
The dimensionality of the dynamical systems constructed using these variables is equal to the number of matter sources, making a very general analysis of two-fluid cosmologies possible. In the past, such models required a higher dimensional approach. In fact, even without specifying a given model f , we have shown it is possible to determine all of the critical points of the system (Table II) as well as most of their physical properties (Table III). For example, we generally find stable points corresponding to the de Sitter solutions with effective equation of state parameter w eff = −1 and deceleration parameter q = −1. This leaves little room for models that exhibit entirely different qualitative properties than those studied here, although we note that some of the points' properties depend on the free functions m(Z) and n(Z) (which led to a phantom equation of state at P 1 in the Beltrán model). We also note that more complicated functions f will lead to a more complicated modified Friedmann constraint equation. It is the latter that determines the physical phase space for each model, and this depends on the values of all parameters in f . It is also this constraint equation that determines which of the critical points are physically important, and this can take a slightly convoluted form, as in one of the models we discussed.
Our model independent approach can be used to discriminate between contending models and constrain their free parameters. To discriminate between the classes of geometric theories themselves, one must look beyond the background cosmological dynamics, which are equivalent. Any viable cosmological background solution, for instance in f (T ) gravity, can also be realised in one of the other geometric formulations. As is clear from the approach taken here, this work is limited to considering just those background dynamics, but it would be interesting to consider how they differ beyond this level.
Despite the large freedom in the choice of functions f and models that can be explored, we have shown that the qualitative dynamics remain simple enough to be captured in the general analysis performed here. With this in mind, one could look to see if observational constraints on physical quantities could be used to constrain the theory space of all possible allowed models. At the simplest level, this would mean only considering models exhibiting a late-time de Sitter point, along with the matter saddle and an early-time radiation repeller. A more detailed treatment would be to use the effective equation of state and deceleration parameters to rule out models that are in conflict with observations, or to give tighter constraints on the free parameters in those models.