Cosmological Constraints on $f(G)$ Dark Energy Models

Modified gravity theories with the Gauss-Bonnet term $G=R^2-4R^{\mu\nu}R_{\mu\nu}+R^{\mu\nu\rho\sigma}R_{\mu\nu\rho\sigma}$ have recently gained a lot of attention as a possible explanation of dark energy. We perform a thorough phase space analysis on the so-called $f(G)$ models, where $f(G)$ is some general function of the Gauss-Bonnet term, and derive conditions for the cosmological viability of $f(G)$ dark energy models. Following the $f(R)$ case, we show that these conditions can be nicely presented as geometrical constraints on the derivatives of $f(G)$. We find that for general $f(G)$ models there are two kinds of stable accelerated solutions, a de Sitter solution and a phantom-like solution. They co-exist with each other and which solution the universe evolves to depends on the initial conditions. Finally, several toy models of $f(G)$ dark energy are explored. Cosmologically viable trajectories that mimic the $\Lambda$CDM model in the radiation and matter dominated periods, but have distinctive signatures at late times, are obtained.


I. INTRODUCTION
The discovery that our universe recently started accelerating [1,2,3,4] has baffled particle physicists and cosmologists for a decade. Many suggestions about the origin of the cosmic acceleration, or the nature of dark energy, have been proposed (for a review, see [5]). While the cosmological constant (ΛCDM) introduced by Einstein is the simplest proposal, albeit requiring fine-tuning to avoid coincidence problems [6], other exotic negative-pressure fluids, often described in terms of scalar fields [7], have been proposed to address these issues. However, it is probably fair to say that no compelling and well-motivated solutions have yet been developed. As a result, attention has naturally turned to the possibility of dark energy originating from modifications of gravity itself [8,9,10,11,12,13].
Among the recent modifications of general relativity, perhaps the simplest approach is to replace the Einstein-Hilbert Lagrangian density with some general function of the Ricci scaler √ −gf (R) [9,10,14,15,16,17,18,19]. For these kind of models, generally it is easy to obtain a late time accelerated epoch. In [20,21], conditions for the cosmological viability of general f (R) models have been established. But as with the popular model f (R) = R+µ/R, they are also generally found to be inconsistent with solar system constraints [22,23].
A well motivated curvature invariant beyond the Ricci scalar is the the Gauss-Bonnet term which is inspired by string theory [24,25] and has gained special interest in cosmology [26,27,28,29,30,31]. It is also well known that the Gauss-Bonnet combination can reduce the number of spin-2 ghosts, which otherwise severely haunt the perturbation theory [32]. The 4 dimensional Gauss-Bonnet term is a topological invariant, thus has no dynamical effects if added into the Lagrangian linearly. To introduce extra dynamics, one may couple the Gauss-Bonnet term to a scalar field, as naturally appears in low-energy string effective actions [24]. For the exponential coupling with a scalar field potential, this model can produce a matter dominated period followed by an accelerated period [33,34]. However, it has been shown that it tends to conflict with solar system constraints [35]. Alternatively, we may consider the Lagrangian density as some general function of the Ricci scalar and the Gauss-Bonnet term √ −gF (R, G) [13,31], and as a natural simplification, F (R, G) = R/2 + f (G) has been investigated [30,36,37,38,39]. Uddin et al. [40] studied the existence and stability of power-law scaling solutions in f (G) models, finding that scaling solutions exist in the model f (G) = ±2 √ αG with α an arbitrary constant. In [41], the authors derived the stability conditions for a de Sitter solution and a standard matter/radiation solution in general f (G) models; the authors also constructed several cosmologically viable f (G) models, but found it difficult to simulate the models to high redshifts.
In this paper, we systematically survey the cosmological viability of general f (G) models. The layout of this paper is as follows. In Section II, we first derive field equations of general f (G) models in the flat Robertson-Walker background and recast them nicely into a 5 dimensional autonomous system. Analogous to the case of f (R) models [20], we construct the curve m(r), where m = Gf GG /f G and r = −Gf G /f with f G = df /dG, etc.. In Section III, the properties of the phase space and their cosmological implications are determined. Interestingly, instead of critical points, we obtain lines of critical points in the phase space. We then establish conditions for cosmologically viable f (G) dark energy models in the m-r plane. Consequently, when deciding the cosmological viability of a particular f (G) model, one's first resort may be to inspect its m(r) curve in the m-r plane. The particular results of [41] can be easily recovered from the more general conditions obtained with our approach. In Section IV, we investigate several specific f (G) models and examine their cosmological viability by numerical simulations. Cosmologically viable trajectories of the universe are obtained. We summarise the results in Section V.
We consider the following action where R is the Ricci scalar, f (G) is a general function of the Gauss-Bonnet term G and √ −gL r and √ −gL m are the radiation and matter Lagrangian densities respectively. We have chosen 8πG b = 1, where G b is a bare gravitational constant. Varying the action (2) with respect to g µν , we obtain the field equations for f (G) gravity where T µν is the energy momentum tensor of radiation (energy density ρ r ) and matter (energy density ρ m ). In the observationally favoured flat Robertson-Walker metric [4] we have R = 6(Ḣ + 2H 2 ) , and Eqs. (3) reduce to where H is the Hubble parameter and an overdot stands for a derivative with respect to t. Additionally, the densities ρ r and ρ m satisfy the usual continuity equationṡ Eqs. (6)(7)(8)(9) determine the dynamics of the f (G) gravity system (2) in a homogeneous and isotropic background. Drawing the analogy with 3H 2 = ρ and −2Ḣ = p + ρ, we can naturally define gravitationally induced forms of dark energy density ρ DE and pressure p DE as In this way, the continuity equation of dark energy holds, and its equation of state parameter becomes We also define the effective equation of state [20] and dimensionless energy densities To discuss the dynamics of a general f (G) model, it is convenient to rewrite the equations of motion as a dynamical system [42,43,44]. To achieve this we introduce the following dimensionless variables Then Eqs. (6)(7)(8)(9) can be recast as a first order dynamical system where N = ln(a/a i ) (a i is the initial value of the scalar factor) and Note that Eq. (6) is recast as In terms of x i , we can rewrite w DE and w ef f as From the r.h.s. of Eq. (24), we can factor out x 4 , which suggests that x 4 = 0 is an invariant submanifold [45] of the dynamical system, meaning the system can not go through the subspace x 4 = 0 and can only approach it asymptotically. We can also factor out x 5 from the r.h.s. of Eq. (25). However, x 5 = 0 is not an invariant submanifold, because when x 5 = 0 (i.e., G = 0), the first term (−x 3 x 2 5 /x 1 m =Ġ/24H 5 ) on the r.h.s. of Eq. (25) is not necessarily zero.
From Eq. (27) we can express G as a function of x 1 /x 2 , and then substituting it into Eq. (26), m can be expressed in terms of x 1 /x 2 . Therefore, given a form of f (G), we obtain a function m(r = x 1 /x 2 ) and the dynamical system (21-25) becomes autonomous. For example, the model f (G) = α(G p − β) q corresponds to the straight line m(r) = (1 − q)r/q + p − 1 in the m-r plane. A similar function of m(r) exists in the f (R) models and examining the cosmological viability of f (R) dark energy models according to the corresponding m(r) curve in the m-r plane has proven to be a fruitful approach [20]. On the other hand, if we know a form of m(r) that is cosmologically viable, one may use it to obtain f (G) from Eqs. (26) and (27).
From Eqs. (22) and (25), given a form of f (G), we may express G and H in terms of x 2 and x 5 , and substituting it into Eq. (21), we may then express x 1 in terms of x 2 and x 5 , so x 1 is actually not an independent variable. However, beyond the power-law form f (G) = αG n , which allows x 1 to be eliminated directly, even simple forms of f (G) yield a complicated dependence for x 1 on x 2 and x 5 , which is often impossible to obtain in closed form. Note that in [41], the authors recast the field equations into a 4 dimensional dynamical system (autonomous if the form of f (G) is given). However, phase space analysis, especially for general f (G) models, in that reformulation is not easy. Consequently we solve the original 5 dimensional system (21)(22)(23)(24)(25), which is analytically simple and numerically more reliable and efficient.

III. ANALYTICAL RESULTS OF GENERAL f (G) MODELS
Writing the field equations in terms of the new dimensionless variables x i makes it easy to look at the analytical properties of general f (G) models. The critical points (equilibria) are where the r.h.s. of the first order dynamical system (21-25) is zero. A new feature emerges here. Rather than having isolated critical points, we obtain continuous lines of critical points (equilibrium manifolds), which we call critical lines: L 1 , L 2 and L 4 are straight lines in the phase space, while L 3 is not. Since L 1 , L 2 and L 4 have different constant values of x 5 , 1, −1/2 and −1 respectively, they do not intersect. However, Generally, if a nonlinear system has a critical line, the Jacobian matrix of the linearised system at a critical point on the line has a zero eigenvalue with an associated eigenvector tangent to the critical line at the chosen point. This kind of nonlinear system is a special subclass of the non-hyperbolic system (whose linearised system has one or more eigenvalues with zero real parts.). The emergence of equilibrium manifolds may be due to some symmetry in the nonlinear system or that the nonlinear system can be reduced to a lower dimensional one. For the case of f (G) = αG n , as we shall see in Section IV, when the system is reduced to a lower 4 dimensional one with just x 2 -x 5 , L 1 shrinks to a point but the other 3 lines remain. The stability of a particular critical point on the line can be determined by the nonzero eigenvalues, because near this critical point there is essentially no dynamics along the critical line (i.e., along the direction of the eigenvector associated with the zero eigenvalue), so the dynamics near this critical point may be viewed in a reduced phase space obtained by suppressing the zero eigenvalue direction. It is also interesting to note that the dynamical system (21)(22)(23)(24)(25) does not have any parameters in it, but, as we shall see below, the stability of a critical point changes along the critical line. Thus this system might be considered as the case of bifurcation without parameters [46].
Note that for L 2 , L 3 and L 4 to exist, m should be equal to −1/2, as we have noted in the definitions of the critical lines above. Also, for L 2 , L 3 and L 4 , we have r = x 1 /x 2 = −1/2. Hence for these critical lines, is satisfied, i.e., they must be located at the point (−1/2, −1/2) in the m-r plane. Thus for L 2 , L 3 and L 4 to exist, r = −1/2 should solve the equation m(r) = −1/2. From the values of Ω DE , w DE , etc., we can see that L 2 and L 4 correspond to solutions in which dark energy scales with matter and radiation respectively, thus they can hopefully correspond to the standard matter and radiation epochs. Therefore, for a standard matter or radiation epoch to exist, the constraint (31) should be satisfied, i.e., the m(r) curve should pass through (−1/2, −1/2) in the m-r plane. Following the case in general f (R) models [20], we can derive the following equation from Eqs. (26) and (27) dr When r = 0 is satisfied and (m + r + 1)Ġ/HG does not diverge or m(r) = −r − 1 is satisfied and rĠ/HG does not diverge, we have dr/dN = 0. From the m-r plane's view, this means that, evolving along the curve m(r), typically the system can not go through any intersection points between the curve m(r) and the m axis (r = 0) or the particular straight line m c (r) = −r − 1. However, the curve m(r) can be "connected" at r = ±∞ or m = ±∞. So for example, the system may evolve from the r < 0 half plane to the r > 0 half plane by "tunnelling" through the infinite point of the m(r) curve at r = ±∞. This happens because basically r is not an essential dynamical variable of the dynamical system (21)(22)(23)(24)(25). For example, r will go through the infinite point if x 2 goes from 0 − to 0 + with x 1 remaining finite.
Let us now discuss the properties and cosmological consequences of an arbitrary point on each of the critical lines in turn. Note that supposing the system starts from somewhere near a particular critical line and the evolution of the system is initially attracted to the critical line, then the particular critical point on the line the system will finally evolve to or pass nearby depends on the initial conditions. We denote an arbitrary point on, for example, L 1 as P 1 .
This is a de Sitter point for any point along the line L 1 , i.e., for any value of x 20 (a given value of x 2 ). The eigenvalues of the linearised system are: where m 1 = m(r 1 = (1 − x 20 )/x 20 ). Hence this is a stable spiral in the subspace of the last two eigenvalues when a stable node when and an unstable point otherwise. In the m-r plane, the stable spiral condition (34) corresponds to the area enclosed by the r axis and the curve m dS * (r) = 8(r + 1)/25r, whilst the stable node condition (35) corresponds to the area enclosed by the curve m dS * (r) = 8(r + 1)/25r and the curve m dS (r) = (r + 1)/2r. Note that the point (−1/2, −1/2) is on the curve m dS (r) = (r + 1)/2r (the boundary of the stable area), so the stable area of L 1 is connected to L 2 , L 3 and L 4 in the m-r plane, as seen in Fig. 1. One can check that these conditions are consistent with those established for a stable de Sitter point in [41]. For example, from the condition (34), we can derive the condition for the existence of a stable spiral in [41] 0 where H 1 and G 1 are the Hubble parameter and the Gauss-Bonnet term at the de Sitter point respectively and we have used , scaling with matter point At this point, dark energy mimics the evolution of matter, and the density of dark energy scales with that of matter (Ω m /Ω DE = (6 − 5x 30 )/5x 30 ), whereas there is no radiation. The eigenvalues of the linearised system are: where m (−1/2) = dm(r)/dr| r=−1/2 . Along this line, the stability of the last two eigenvalues changes at the points x 30 = 0 or x 30 = 96/71, which we may call bifurcation points: the points with x 30 < 0 or x 30 > 96/71 are stable, while other points are unstable. In order for P 2 to be a standard matter era (Ω m 1), x 30 needs to be near 0, but x 30 = 0 is a singularity of the last two eigenvalues. It is interesting to note that the situation here is quite similar to that of the standard matter point in general f (R) models [20], where for P 5 in their notation to be a standard matter era m 5 needs to be near 0, but 0 is a singularity of two eigenvalues. For our current point P 2 , if x 30 → 0 + , the last two eigenvalues diverge as −∞ and +∞. However, if x 30 is some small (non-infinitesimal) positive value, the two eigenvalues are large and finite with opposite signs. If the system flows towards these kind of points, unless the initial conditions are extremely fine-tuned, the system will not remain around it for a reasonable long time. Since an elongated matter dominated epoch is needed for cosmic structure formation, we generally consider this case as unacceptable. On the other hand, if x 30 → 0 − , then the last two eigenvalues are complex with negative real parts, which means the two eigenvalues are stable. However, the imaginary parts diverge, and hence the system oscillates rapidly in the subspace associated with the two eigenvalues when it approaches the standard matter point. The frequency of the oscillation can be reduced if x 30 is slightly less than 0, while the amplitude of the oscillation will not be very large if the initial conditions are set near the critical point in the associated subspace. The system should finally leave the standard matter point and enter an accelerated point, so we shall set 3(m (−1/2) + 1) to be greater than 0 in order to make the matter point unstable, i.e., In the m-r plane, it means that there are unacceptable directions of the m(r) curve approaching the point (−1/2, −1/2), see Fig. 1. Note that the properties and their cosmological implications obtained above assume that the system is exactly at a critical point. Since the critical point here for the standard matter epoch is not supposed to be totally stable, the system will just pass nearby the critical point rather than fall onto it, so the properties and their cosmological implications when the system passing by the critical point should be just approximate to what are inferred from the critical point analysis.
x 50 −2 , 0, x 50 ), dark energy dominated point This is a dark energy dominated point with no radiation and no matter. It has a fascinating wealth of possibilities, mimicking radiation when x 50 = −1 and matter when x 50 = −1/2; it is a quintessence-like point (−1 < w DE < −1/3) if 0 < x 50 < 1 and a phantom-like point (w DE < −1) if x 50 > 1. The eigenvalues of the linearised system are: This is a stable point when either of the following two conditions is satisfied or Note that the second condition (43) is consistent with the condition (40), which means that if we require an unstable direction for the system to leave the standard matter point P 2 , then P 3 can be a stable phantom-like point, co-existing with the stable de Sitter point P 1 . Which point the universe finally evolves to depends on the initial conditions, as we shall see in Section IV. Note that although dark energy with w < −1 is problematic theoretically [47], it is still consistent with current data [4,48,49].
, scaling with radiation point At this point, dark energy mimics the evolution of radiation, and the density of dark energy scales with that of radiation (Ω r /Ω DE = (4 − 3x 30 )/3x 30 ), whereas there is no matter. The eigenvalues of the linearised system are: The stability structure in the subspace of the last two eigenvalues is similar to that of P 2 : the points with x 30 < 0 or x 30 > 64/47 are stable in the subspace (still unstable in the total phase space); in order for P 4 to be a standard radiation point (Ω r 1), x 30 needs to be near 0, but at the same time, x 30 = 0 is a singularity of the two eigenvalues; if x 30 is slightly greater than 0, P 2 is extremely unstable; if x 30 is slightly less than 0, P 2 is stable in the subspace but the system oscillates significantly when flowing to it. Note that the condition f GG > 0 from [41] for the existence of a standard radiation/matter point can also be recovered. For example, to avoid the violent instability of the standard radiation point, we require x 30 be slightly less than 0, so −f /3H 2 = x 20 = −x 30 /2 > 0, which means f < 0. Additionally, for the standard radiation point, we have r = −Gf G /f = −1/2 < 0 and m(−1/2) = Gf GG /f G = −1/2 < 0, so f GG > 0 holds at the standard radiation point.
Starting from the radiation dominated era, a cosmologically viable trajectory of the universe in the phase space would start somewhere near the standard radiation point P 4 (with x 30 slightly less than 0), then slowly pass nearby the standard matter point P 2 (with x 30 slightly less than 0) and finally land at the stable de Sitter point P 1 or the stable phantomlike point P 3 . In the m-r plane, with the trajectory fixed along the curve m(r) and subject to the constraint from Eq. (32), the universe would start near the point (−1/2, −1/2) where the standard radiation point P 4 and the standard matter point P 2 are located, slowly moving away from it, and finally fall onto some point in the stable area of P 1 or come back to evolve to P 3 that also resides at the point (−1/2, −1/2). Note that it might also be possible that the universe has an unstable accelelated epoch after the matter dominated period, whose cosmological viability is not the concern of this paper and is left for future work. Indeed, for phase spaces higher than 2 dimensions, more complicated dynamical phenomena such as chaos may appear, but again this is beyond the scope of the current paper.

IV. SPECIFIC MODELS AND NUMERICAL RESULTS
Having described possible cosmic trajectories for general f (G) dark energy models, we now turn our attention to the cosmological viability of a few specific toy models of f (G) whose m(r) curves can be analytically obtained. Note that if the analytical form of m(r) can not be obtained, one may get a numerical expression of it for a certain range of r. On the other hand, given a form of m(r), it may not be easy to obtain an analytical form for f (G), but again we can solve it numerically.
The Gauss-Bonnet term G evolves from negative values in the radiation and matter dominated epochs to positive values in the acceleration epoch, thus the lagrangian is not well defined if there is any function of the Gauss-Bonnet term such as √ G or ln(βG). In such cases, G in the f (G) should be understood as |G| in this section as well as in Fig. 1. We will start the system deep in the radiation dominated epoch, which requires we initially set r = x 1 /x 2 −1/2, x 4 1, x 5 −1 and x 1 + x 2 + x 3 0. Since the matter dominated point is located at (−1/2, −1/2), setting r = x 1 /x 2 −1/2 is also a requirement for the universe to have an elongated matter dominated period. We find that the dynamical systems of f (G) models are quite often stiff. The simulation figures in this section are produced using the Matlab stiff ordinary differential equation solver ode15s.
We first consider the simplest case f (G) = αG n . This model gives a point (r, m) = (−n, n − 1) rather than a curve in the m-r plane, see Fig. 1.
So we may consider this system in the 4 dimensional subspace x 2 -x 5 . In this subspace, while L 2 , L 3 , L 4 remain critical lines, L 1 shrinks to a critical point (1/(1−n), 0, 0, 1). In order to have the standard matter and radiation points, the point (−n, n − 1) should be located at (−1/2, −1/2), thus n = 1/2. Note that this result is consistent with what has been obtained in surveying the scaling solutions of f (G) models [40]. So in this case, all the 4 critical manifolds (points or lines) are at (−1/2, −1/2). Re-calculating the eigenvalues, we find that, for P 1 , they are −4, −3, 0, −3; for P 2 , P 3 and P 4 , the eigenvalues containing m (−1/2) disappear, compared to the 5 dimensional case. The 0 eigenvalue for P 1 here is not related to any critical line, thus the stability of this isolated non-hyperbolic point may be obtained by the centre manifold method [45]. Since the eigenvalue containing m (−1/2) for P 2 has vanished, this point is either totally stable (the system will stay in the matter dominated era for ever and not evolve to the dark energy dominated era.) or extremely unstable (the system can not stay in the matter dominated era for long enough to allow cosmic structure to form.). Therefore this model is ruled out unless we allow extreme fine-tuning of initial conditions.
The power-law model may be generalised to f (G) = α(G p − β) q . The corresponding curve m(r) and its first derivative m (r) for this model are respectively To have a standard matter point, i.e., to employ the relation (31), we find it leads to a constraint equation for the parameters 2pq = 1 .
The system needs to be able to leave the matter point, which leads to a further constraint Thus if these constraints are satisfied, cosmologically viable trajectories can in principle be obtained.
To be specific, we set p = 3/4 and q = 2/3 (see Fig. 1) and numerically solve this model to get cosmic trajectories. In Fig. 2, we obtain a trajectory that is ended with a stable de Sitter era. For this trajectory, while the behaviours of Ω r , Ω m , Ω DE and w ef f are quite similar to those of the ΛCDM model in the radiation and matter dominated epochs, w ef f oscillates rapidly and goes deep into the phantom-like region as the universe enters the de Sitter period. On the other hand, w DE at first mimics the background fluids (although oscillating significantly when leaving the radiation period), then at the end of the matter dominated epoch plunging below the w = −1 divide for a while, thereafter increasing above the divide, and finally oscillating rapidly into the de Sitter period. In the m-r plane, confined on the straight line m = 1/2 · r − 1/4, the universe slowly leaves the point (−1/2, −1/2) towards the r = 0 line, bounces back and forth and then goes into the stable de Sitter point with r < 0. The reason why w DE as well as w ef f oscillate rapidly when the system enters the de Sitter period is related to the fact that the line r = 0 actually corresponds to the bifurcation point of L 1 , similar to the bifurcation point of L 2 that is discussed in details in Section III. To see this, we express the last two eigenvalues of P 1 in terms of r, and we find that these eigenvalues have divergent imaginary parts if r → 0 − . Note that the system may also evolve to the stable phantom-like point P 3 with x 50 > 1 after the matter dominated period, see Fig. 3. In this case, w DE oscillates rapidly in the radiation and matter dominated epochs, then plunges well below the w = −1 divide at the end of the matter dominated epoch and thereafter shortly passes nearby the de Sitter point P 1 (with rapid oscillation) and the point P 3 with x 50 0 before going into the stable phantom-like point. In the m-r plane's view, the system slowly leaves the point (−1/2, −1/2) towards the r = 0 line and then bounces back and falls onto the point (−1/2, −1/2).

V. CONCLUSIONS
In this paper we have derived conditions for cosmologically viable f (G) models. By recasting the field equations of f (G) gravity with a flat Robertson-Walker metric into a 5 dimensional autonomous system, we have seen that the equilibria of the system are not isolated but form so-called equilibrium manifolds (in the present case they are 1 dimensional lines, called critical lines.). We have shown that the emergence of the critical lines is not totally due to the fact that the 5 dimensional autonomous system can in principle be reduced to a 4 dimensional one. Discovering 4 critical lines in the phase space for a general f (G) model, the curved line L 3 is found to connect to the other 3 straight lines, L 1 , L 2 and L 4 . We have also found that for all 4 lines the stability of a critical point changes when moving along the critical line, which may be referred to as bifurcation without parameters. For the critical line L 2 (or L 4 ) in which the dark energy density scales with that of matter (or radiation), matter (or radiation) becomes dominant only if the parameter of the critical line x 30 is around 0. At the same time, the stability of the critical line changes rapidly around x 30 = 0, and we require x 30 < 0 to avoid a violent instability (i.e., to avoid extreme fine-tuning of the initial conditions).
Similar to the case of f (R), for a specific f (G) model, we can construct the m(r) curve in the m-r plane. In this plane, L 2 , L 3 and L 4 all lie at the point (−1/2, −1/2). Thus for a standard matter point P 2 (on L 2 ) to exist, we require the m(r) curve satisfy the relation (31): The condition (40): is also employed to let the universe finally leave the matter dominated point, which corresponds to some forbidden directions around the point (−1/2, −1/2) in the m-r plane. The stable condition for the de Sitter point P 1 (on L 1 ) corresponds to the area enclosed by the curve m dS (r) = (r + 1)/2r and the r axis in the m-r plane, in which the area enclosed by the curve m dS * (r) = 8(r + 1)/25r and the r axis corresponds to a stable spiral (in a 2 dimensional subspace) with the rest corresponding to a stable node. P 3 (on L 3 ) is a dark energy dominated point. It has a wealth of possibilities. Besides mimicking radiation when x 50 = −1 and matter when x 50 = −1/2, it can be a stable quintessence-like point or a stable phantom-like point. But if we require a standard matter epoch (requiring the condition (40) to hold unless we allow extreme fine-tuning of initial conditions), it can only be stable as a phantom-like point. Which stable accelerated point the universe finally evolves to, the de Sitter point P 1 or the phantom-like point P 3 , depends on the initial conditions. We have also studied several toy models of f (G) dark energy whose m(r) curves can be obtained analytically. Note that if the analytical form of m(r) can not be obtained, one may obtain it numerically. On the other hand, if a form of m(r) can be obtained from observations or other kinds of analyses, it may not be easy to derive the corresponding f (G), but still we might find it numerically. For the model f (G) = α(G p − β) q in particular, cosmologically viable trajectories have been obtained for the case with a de Sitter epoch as the final stage as well as for the case with a phantom-like epoch as the final stage. These trajectories mimic the ΛCDM scenario in the radiation and matter dominated periods (although w DE may oscillate for a certain period) but have distinctive signatures at late times. Most significantly, w DE plunges below the w = −1 divide in the late matter dominated epoch (pole-like behaviour, for similar behaviours of w DE in other scenarios see [50]) and then oscillates rapidly, also going below the w = −1 divide, when evolving to (for the case with a de Sitter epoch as the final stage) or passing nearby (for the case with a phantom-like epoch as the final stage) the de Sitter point P 1 . For other potentially viable models, cosmologically viable trajectories may still exist, but we find it is non-trivial to obtain them. We leave this task for future work.