Friedmann’s Universe Controlled by Gauss-Bonnet Modified Gravity

The accepted idea that the expansion of the universe is accelerating needs, for compatibility to general relativity, the introduction of some unusual forms of matter. However, several authors have proposed that instead of making weird hypothesis on some yet unobservable species of matter, one should follow the original idea of the first Einstein’s paper on cosmology and consider that in the cosmic scene one has to modify the equations that control the gravitational metric. This possibility led us to re-examine the evolution of the topological invariant containing two duals in a dynamical universe, the so called Gauss-Bonnet topological invariant. The particular interest on this invariant is due to the fact that in a homogeneous and isotropic universe this invariant drives the cosmic acceleration. In a decelerating scenario and as a necessary previous condition of an ulterior acceleration this invariant must have an extremum identified to its maximum value. We will examine the conditions for this to occur, and a description of the universe with epochs of accelerated and decelerated expansion.

Modified theories of gravity have been proposed in order to account for astronomical data, like the cosmic late acceleration. Two of such propositions that have been gaining notoriety are the f(R) theory and the Gauss-Bonnet (GB) topological invariant modified gravity [1−4]. In GB modified gravity it makes use of the topological invariant of the same name, which is only an invariant in 4D space-time, commonly expressed as $ {\cal{G}} $.

This topological invariant $ {\cal{G}} $ by itself does not affect the dynamics of General Relativity in 4D when added to the Einstein-Hilbert action, instead, there are three basic ways for implementing it. First, by coupling $ {\cal{G}} $ to a scalar field ϕ [5−7]. Second, by considering a nontrivial function of the topological invariant, $ f({\cal{G}}) $ modified gravity. And third, by implementing the topological invariant to a Einstein-Hilbert action of higher dimension D, in which $ {\cal{G}} $ is no longer a topological invariant, and taking the limit $ D\rightarrow4 $, yielding nontrivial modified dynamics [8]. Some groups had used observational data in order to find constraints to different approaches to GB modified gravity [9−13]. Great reviews of modified theories of gravity can be found at [14−18].

Studies in $ f({\cal{G}}) $ gravity have shown its capabilities of mimicking the λCDM cosmology and describe radiation/matter epochs followed by cosmic acceleration [19]. The main critique for $ f({\cal{G}}) $ is due to the appearance of higher order derivative terms in the equations, which may lead to instabilities in the theory in first order perturbation models [20, 21]. This requires the establishment of criteria and conditions to exclude possible unphysical solutions (see [5] for example).

In this present work we use the $ f({\cal{G}}) $ Gauss-Bonnet Modified Gravity approach, by applying the simplest non-trivial function of $ {\cal{G}} $, we peer into its consequences on a homogeneous and isotropic Friedmann's geometry, first in an empty space-time, then in the case with source matter as cosmic dust. Section II. gives a very brief definition of topological invariants in 4D. Section III. presents the standard $ f({\cal{G}}) $ gravity; subsection III.a. presents a relationship between $ f({\cal{G}}) $ and the Fierz representation of a spin-2 particle. Section IV. applies the theory to the Friedmann's universe; subsection IV.A. deals with the case for $ {\cal{G}} = constant $, which yields a number of important constraints to filter out non-physical solutions; subsection IV.B. considers the simplest non-trivial function of $ {\cal{G}} $, that is $ f({\cal{G}}) = \sigma {\cal{G}}^2 $; first we assume an empty space-time, which produces a 2D dynamical system, and a complete qualitative description of numerical solutions near and around the critical points is presented; and second, we assume cosmic dust as source, producing a 3D dynamical system, and present a set of qualitative behavior of numerical solutions. Section V, the conclusion.

There are two topological invariants in a 4-d Riemannian geometry, that is

where

and $ {\cal{G}} $ is the Gauss-bonnet topological invariant

The dual is defined for any anti-symmetric tensor $ F_{\mu\nu} $ as

where $ \epsilon_{\mu\nu\alpha\beta} $ is the completely anti-symmetric Levi-Cività symbol.

We note that the invariant $ {\cal{G}} $ can be written in terms of the curvature tensor (without the dual operation) by the identity

This $ {\cal{G}} $ satisfies the identity,

For an arbitrary function $ f({\cal{G}}) $, the variation of the action

where $ L_m $ is the matter Lagrangian, leads to the following set of dynamic equations,

where $ Z_{\mu\nu} $ is defined as

where $ f' = \dfrac{df}{d{\cal{G}}} $.

Note that contraction of the last two indices of $ {H}^{\mu\nu\alpha} $ implies the identity

Let us show that such quantity $ {H}^{\mu\nu\beta} $ defined by equation (12) may be associated to a spin-2 field. We start by remembering that there are two basic representations to deal with a spin-2 field that we call the Einstein and the Fierz representations (both were introduced by Fierz):

● Einstein representation: $ \varphi_{\mu\nu}, $

● Fierz representation: $ F_{\alpha\beta\mu}. $

The tensor $ \varphi_{\mu\nu} $ is symmetric and it has 10 independent components. The tensor $ F_{\alpha\beta\mu} $ has the symmetries

which states that this tensor has 24 independent components. The second property states that the dual has no trace, that is,

or, equivalently,

which reduces this number to 20. The identity

makes this number to 14, and finally, the property of vanishing trace,

reduces the number of independent components to 10.

It is a direct exercise to show that the tensor $ {H}^{\mu\nu\alpha} $, defined in equation (12),

satisfies all the conditions above, showing that there is a spin-2 tensor hidden in the above dynamics, constructed with the invariant $ {\cal{G}} $.

The tensor $ {H}^{\mu\nu\beta} $ has four properties:

● Anti-symmetric in the first two indices: $ {H}^{\mu\nu\beta} = -{H}^{\nu\mu\beta} $

● Cyclic identity: $ {H}_{\mu\nu\beta}+{H}_{\beta\mu\nu}+{H}_{\nu\beta\mu} = 0 \Leftrightarrow {H^{\overset{\ast}{\mu\alpha}}}{_\alpha} = 0 $

● Null divergent in the last indice: $ {H^{\mu\nu\beta}}{_{;\beta}} = 0 $

● Contraction of the last two indices implies $ :{H^{\mu\nu}}{_\nu} = -X_{\nu} \left({R}^{\mu\nu}-\frac{1}{2}R{g}^{\mu\nu}\right) $

So the quantity $ {H}^{\mu\nu\beta} $ satisfies all these conditions, showing that it could be associated to a spin-2 tensor hidden in the above dynamics. This is only a formal property of the quantity $ H^{\mu\nu\alpha} $ related to derivatives of second order of the metric tensor.

Let us consider the homogeneous and isotropic Friedmann's metric,

where $ a(t) $ is the scale factor. Using the definition of the Hubble parameter $ H\equiv \dot{a}/a $, the Gauss-Bonnet topological invariant yields

note that the sign of $ {\cal{G}} $ indicates the sign of the acceleration of the scale factor.

From (10), the new term $ Z_{\mu\nu} $ can be treated as a geometric source, and in case of the metric (20), when $ {\cal{G}} $ is only a function of time, it can be associated to a perfect fluid [22], that is,

where $ Z_{\mu\nu} $ can be written under the form

Thus, the corresponding geometric density $ \rho_z $ and geometric pressure $ p_z $ are given by,

and

The modified dynamical equations from (22) are straightforward,

where the energy-matter source satisfies the conservation law,

Observe that, differentiating (26) with respect to time, and applying equations (21) and (28), results in (27), that is, (27) is not an independent equation and can be interpreted as a cinematic relation with respect to the expansion parameter of the model.

In order to observe qualitative and quantitative properties of the model we make use of the theory of dynamical systems. Choosing the variables $ {\cal{G}} $, H and $ \rho_m $, and through equations (21), (24), (26) and (28), the problem can be enunciated in a general dynamical system form:

where (31) refers to one class of matter-energy source, and when other independent sources of matter-energy are added to the system there are the corresponding source density dynamical variables and equations. In addition, it is necessary to define the equation of state relating $ \rho_m $ and $ p_m $, usually of the form $ p = \omega \rho $, for a given value of ω.

One class of functions for the topological invariant $ {\cal{G}} $ to consider is of the type $ f({\cal{G}}) = \sigma{\cal{G}}^n $, where σ is a constant parameter and $ n = 2, 3, ... $, in order to give rise to nontrivial effects. The choice of a particular $ f({\cal{G}}) $ only affects equation (29),

It is apparent that for $ n>2 $ the dynamical system presents a singularity when $ {\cal{G}} = 0 $ in the phase space, which may provide a barrier for the system to transition between positive and negative values of $ {\cal{G}} $. Then, $ n = 2 $ provides the simplest case in order for $ {\cal{G}} $ to have the ability to change sign, and consequently, give rise to a model with epochs of accelerated and decelerated rates of expansion. We will see later that this system can have solutions where the Gauss-Bonnet topological invariant $ {\cal{G}} $ is cyclically dampened, changing it's sign periodically, which could be used to describe a universe with stages of acceleration and deceleration, but first let us consider the case where $ {\cal{G}} $ is constant for an interval or indefinitely.

The most simple behavior we could get from $ {\cal{G}} $ is a constant α, let

so $ \dot{{\cal{G}}} = 0 $ and $ \ddot{{\cal{G}}} = 0 $, then the geometric density and pressure (24) and (25) are constants too given by

From the definition of $ {\cal{G}} $ in (21) the problem is to solve the integral,

There are three classes of solutions for (36), depending if α is positive, negative or zero. And there is a fourth class of solutions by assuming H as a constant.

The implicit solution of (36) for $ \alpha>0 $ is

where k is a constant of integration that simply shifts the solution in time. Assuming $ k = 0 $ and taking the limit for $ H\rightarrow \pm\infty $ on the left side of (37), it gives the result

This implies that the Hubble parameter diverges to $ \pm\infty $ in a finite interval $ [-L,L] $. The general solution is expressed in Figure 1a.

Figure 1.  (color online) Solution of $ H(t) $ for $ {\cal{G}} $ as a positive constant.

From equation (26) the source density can be expressed as

this means that if $ \rho_z $ is positive, then the condition $ \rho_m\geq 0 $ is only satisfied in the intervals $ -L<t\leq -m $ and $ m\leq t<L $, where $ \pm m $ is the time t for which the relation $ 3H^2 -\rho_z = 0 $ is satisfied. The value of m depends on the choice of the function $ f({\cal{G}}) $ and α. The conclusion is that, for $ {\cal{G}} $ as a constant α and $ \rho_z $ positive, the evolution of the Hubble parameter H is confined in one of the intervals $ [-L,-m] $ or $ [L,m] $ in order for the source density $ \rho_m $ to be greater or equal to zero, so H can be either positive or negative in this configuration. See Figure 1b for the graph of $ 3H^2 $ and its respective allowed intervals.

The implicit solution of (36) for $ \alpha<0 $ is

Where k is a constant of integration. For $ k = 0 $, the implicit solution above implies two possible solutions for $ H(t) $, one is $ H_1 (t) $ which is always positive, it diverges to $ +\infty $ at a finite past time $ t = -L $, and for the infinite future ($ t\rightarrow\infty $) it stabilizes at the value $ H_f = \left(\dfrac{|\alpha|}{24}\right)^{1/4} $. The other solution is $ H_2 (t) $ which is always negative, it diverges to $ -\infty $ in a finite future $ t = L $, and for the infinite past ($ t\rightarrow -\infty $) it approaches the value $ -H_f $. The solutions $ H_1 $ and $ H_2 $ are the inverse mirror images of each other, obtained by the transformation $ H_2(t) = -H_1 (-t) $, see Figure 2a. The values of $ \pm L $ is obtained by the limits $ H_1\rightarrow +\infty $ and $ H_2\rightarrow -\infty $, given by

Figure 2.  (color online) Solution of $ H(t) $ for $ {\cal{G}} $ as a negative constant.

The functions $ 3H_1^2 $ and $ 3H_2^2 $ have the same minimum value at $ 3H_f^2 $. Using the condition for the positivity of matter source ($ \rho_m \geq 0 $) given by (26), that is, $ 3H^2\geq \rho_z $, it leads to two possibilities: (a) If $ 3H_f^2 \geq\rho_z $ then the solutions $ H_1 $ and $ H_2 $ are valid over their respective time intervals ($ [-L,+\infty] $ and $ [-\infty,L] $); (b) If $ 3H_f^2 \leq\rho_z $ then $ H_1 $ is only valid in the interval $ [-L,-m] $, and the same for $ H_2 $ in the interval $ [m,L] $, where $ \mp m $ is the time for which $ 3H_1^2 = \rho_z $ and $ 3H_2^2 = \rho_z $ respectively, see Figure 2b. The value of $ \pm m $ will depend on the specific choice of the function $ f({\cal{G}}) $ and the value of α.

In the case for $ {\cal{G}} $ as a constant with value zero, it is possible to obtain an explicit solution for the Hubble parameter. Solving (36) for $ \alpha = 0 $ gives the solution for H,

where $ K_1 $ is a constant of integration. As for the scale factor $ a(t) $

with $ K_2 $ a constant of integration.

So if $ {\cal{G}} $ is constant and zero, then the scale factor increases or decreases linearly over time, and the Hubble parameter presents a singularity at $ t = -K_1 $. The condition for the positivity of the matter source ($ \rho_m\geq 0 $), can be obtained explicitly from equation (26),

From (44) we see that if the geometric density $ \rho_{z} $ is less or equal to zero, then $ \rho_m $ is always positive; but if $ \rho_{z} $ is positive, then $ \rho_m $ will be positive only in the interval

which depends on the the initial condition for $ H(t) $.

And finally, there is one last class of solutions obtained by assuming the Hubble parameter as a constant $ H_c $. For $ H_c\neq 0 $, and using the equation of state $ p_m = \omega \rho_m $, the dynamical equation (26) and the conservation equation (28), that is,

can only be satisfied if $ \rho_m = 0 $, since (46) implies $ \rho_m $ as a constant but (47) makes $ \dot{\rho_m}\neq 0 $ if $ \rho_m\neq 0 $. From (46) and (21) we obtain the relationships between $ H_c $, α and $ \rho_z $,

where α must be negative for the quantities to be real numbers. So a constant Hubble parameter can only describe an empty space-time.

It is worth noting that for $ H_c = 0 $ then from (21) it implies $ \alpha = 0 $, and equations (46) and (47) are both satisfied, with the matter source assuming the condition,

where (34) was used. Here we have a description of a static universe with homogeneous matter density $ \rho_m $ as a constant directly controlled by the function $ f(0) $. Since the Hubble parameter is zero, this means the scale factor is a constant $ a_c $ in the metric. The positivity of the matter source imposes the condition,

In summary, when the topological invariant $ {\cal{G}} $ is assumed to be constant, there emerges specific conditions in order for the matter source $ \rho_m $ to be strictly positive, even in some cases the topological invariant $ {\cal{G}} $ cannot be constant indefinitely but must vary over time. These specific conditions are highly dependent on the choice of the function $ f({\cal{G}}) $. This may be useful in the future in order to filter only solutions with $ \rho_m \geq 0 $. The separation of solutions by the sign of $ {\cal{G}} $ is useful because it indicates if the system is under accelerated, decelerated or null acceleration of expansion, according to (21). When the Hubble parameter is assumed constant it is clear the direct influence of the function $ f({\cal{G}}) $ over the system, either dictating the expansion rate in an empty universe, or by controlling the value of the homogeneous source matter density $ \rho_m $ in a static universe.

We now consider the first non trivial term in the Taylor expansion of $ f({\cal{G}}) $, that is, the function $ f({\cal{G}}) = \sigma {\cal{G}}^2 $, with σ as a constant. In this case the geometric density and pressure assume the form,

and we have to solve the set of differential equations,

Let us consider the case with an empty space-time. Under these conditions equation (54) is identically zero. Applying the dynamical system expression of the problem, given by (29) and (30), we have

This 2D dynamical system represents the whole problem in an empty space-time. Its numerical solutions will prove useful in giving a map to look for certain types of structures and behaviors in the complete problem with energy-matter sources, which is a dynamical system of higher dimension.

For $ \dot{{\cal{G}}} = 0 $ and $ \dot{H} = 0 $, the equilibrium points in the phase space $ ({\cal{G}}(t) \times H(t)) $ are $ {\cal{G}}_0 = -\frac{-1}{\sigma^{2/3}}\left(\frac{3}{2}\right)^{\frac{1}{3}} $ and $ H_0 = \pm\dfrac{1}{\left({\sigma 96}\right)^{1/6}} $, these imply the constant σ must be positive in order for the equilibrium points to be real. We already have seen in subsection IV.A.4 that solutions for the Hubble parameter as a constant $ H_0 $ are viable in an empty space-time, which means that in these equilibrium points the scale factor behaves as an exponential manner, given by $ a(t)\sim \exp{(H_0 t)} $. Observe that the line $ H = 0 $ in the phase space is a singularity region, and for the origin $ H = {\cal{G}} = 0 $ of the phase space, the dynamical system is undefined, but from (53) and (55) we see that the equations are identically satisfied. Then, in order to determine the behavior of solutions near the origin, specific methods of analyses must be applied.

The next step is to determine the behavior of solutions near the equilibrium points. For this, we analyze the Jacobian matrix of the system in the equilibrium points, $ a = \left(\dfrac{\partial\dot{{\cal{G}}}}{\partial {\cal{G}}}\right)_{({\cal{G}}_0,H_0)} $, $ b = \left(\dfrac{\partial\dot{{\cal{G}}}}{\partial H}\right)_{({\cal{G}}_0,H_0)} $, $ c = \left(\dfrac{\partial\dot{H}}{\partial {\cal{G}}}\right)_{({\cal{G}}_0,H_0)} $ and $ d = \left(\dfrac{\partial\dot{H}}{\partial H}\right)_{({\cal{G}}_0,H_0)} $, and compute the eigenvalues of this matrix,

It turns out that the two eigenvalues are real and of opposite signs for both equilibrium points, so these equilibrium points are saddle points, that is, trajectories near the equilibrium points are attracted or repelled depending on the direction they approach these points (see [23]).

Now to gain more information on the solutions for this problem we can numerically integrate the solutions for given initial conditions. To numerically integrate the solutions we assumed $ \sigma\equiv 1 $, and using the Runge–Kutta method ([24]), integrating for different initial conditions inside the interval $ [{\cal{G}} = -2...1, H = -1...1] $ which contains the two equilibrium saddle points, gives the behavior map in Figure 3.

Figure 3.  (color online) Integral solutions for the 2D dynamical system in the phase space $ [{\cal{G}}(t)\times H(t)] $. The points $ P1 $ and $ P2 $ are the corresponding equilibrium points of the system. Observe that the $ {\cal{G}} $ axis is a singularity for the system and we cannot infer their behavior numerically.

Numerical solutions are only possible if not near the region of singularity of the system, but through approximation it is possible to determine the behavior of solutions near the singularity region (for $ H \approx 0 $). From (56) and (57) it follows,

For $ |H|<<1 $ and $ |{\cal{G}}|>>1 $, equation (59) simplifies to $ \dfrac{dH}{d{\cal{G}}}\approx 0 $, with solution $ H(t)\approx cte $, that is, in the phase space the $ H(t) $ is constant near the $ {\cal{G}} $ axis when not close to the origin. Let us focus now on what happens to solutions near the origin of the phase space, since $ {\cal{G}}\approx 0 $ and $ H\approx 0 $, the $ H^4 $ term can be dropped out compared to $ {\cal{G}} $ in equation (59), yielding

and the solutions to this system are given by,

and

where a and b are constants that depend on the initial conditions. We note that the value of $ H(t) $ is bounded. The plots of (61) and (62), for different values for a and b, are given in Figures 4a and 4b. Although the dynamical system is undefined at the origin of the phase space, we see that its solutions near it are well defined, the origin acts as a convergence and divergence of solutions, depending on the region it is approached.

Figure 4.  (color online) Solutions near the origin of the phase space $ {\cal{G}} \times H $.

The last step, to understand the whole qualitative picture for the behavior of this simplified problem, is to determine how the numerical solutions given in Figure 3 interact with the solutions near the origin of the phase space (61) and (62). We can take a direct approach, first consider the solution (61) with an arbitrary value for a, because of its square root, both $ {\cal{G}} $ and H have two extreme values in this solution; so choosing one of the extreme values and considering them as initial conditions for the 2D dynamical system (56) and (57), we can numerically integrate to obtain an integral solution. Two things can happen to this integral solution, first it can grow away from the origin of the phase space (see Figure 3), or it can go back to the origin in two ways: it can go back to the original solution (61), in which case the integral solution is cyclic closed; or it can return but for a different value for the constant a in (61), in which case the integral solution is cyclic open, or have a vortex like behavior. The result is, depending on the initial conditions chosen, the integral solution grows away from the origin or it cycles back in a vortex like behavior nearing the origin of the phase space, Figure 5 shows some examples of the first result, and Figure 6 shows an example of the second result.

Figure 5.  (color online) Examples of integral solutions that pass one time through the origin of the phase space $ {\cal{G}} \times H $.

Figure 6.  (color online) One example of an integral solution presenting an open cyclic behavior (vortex like) near the origin, each time passing through the origin with an smaller arch of travel, that is, the origin acts as an attraction pole.

The Figures 3, 5 and 6 give us a summary of the global behavior for the 2D dynamical system (56) and (57). Focusing for H positive, the phase space can be divided into four quadrants, see Figure 7. The first quadrant (I) gives us initial conditions where the Hubble parameter H decreases to a minimum positive value and increases afterwards; the second quadrant (II) have initial conditions with integral solutions that converge to the origin of the phase space and are connected to the fourth quadrant (IV), then for initial conditions in (II) the Hubble parameter decreases to zero and increases afterwards in region (IV); and in region (III) the initial conditions result in integral solutions that are open cyclic around the origin of the phase space, where the origin acts as an attraction pole.

Figure 7.  (color online) Global behavior of the 2D dynamical system(56) and (57) for H positive. Quadrants (II) and (IV) are connected through the origin, and quadrant (III) have open cyclic solutions attracted to the origin.

So now, having a global qualitative behavior for the empty space-time problem, this gives us a map of behavior to look for in the real case with energy-matter sources. The complete problem is given by the set of equations (53) through (55), for an empty space-time this set yields a 2D dynamical system, and adding matter-energy source yields a higher dimension dynamical system. Figures 3 and 7 provide a rough idea of how the integral solutions behave in these higher dimensional cases.

The behavior of the Hubble parameter, and consequently of the scale factor, is highly dependent on the initial conditions, as seen in Figure 7. On quadrant II it can describe a universe initially in an epoch of accelerated expansion, followed by a period of deceleration, and finally one last epoch of steady increase in acceleration and rate of expansion in region IV. Region I offers only the behavior of an ever increasing rate of expansion and acceleration of the system. Solutions on quadrants I and II/IV tend to converge to a single path in the phase space. While region III describe an empty Friedmann's universe in which the topological invariant $ {\cal{G}} $ changes sign periodically, from (21) it implies that the acceleration of the scale factor $ a(t) $ changes sign periodically, but the Hubble parameter H is always positive or negative, in other words, it can describe an empty universe that is always expanding but periodically accelerating and decelerating, and moreover, the intensity of this acceleration and deceleration decreases over time. The points of equilibrium $ {\cal{G}}_0 $ and $ H_0 $ are saddle points, this implies that perturbations in these points may lead to two outcomes, either the system is attracted to the convergence of solutions of quadrants I and II/IV, or is attracted to the periodical solutions of quadrant III. The point $ {\cal{G}} = 0 $ and $ H = 0 $ is a special point of convergence and divergence of solutions, this means that small perturbations to the system near this origin may lead to radical differences in the evolution of the system, hinting to a possibly chaotic nature of behavior. In future work we will present a study when perturbations are applied to the system.

From (26), (27) and (30) we can infer that $ \rho_z = 3H^2 $, and the geometrical equation of state relating the geometrical density and geometrical pressure is given by,

Since the geometrical density $ \rho_z $ is always positive, then the topological invariant $ {\cal{G}} $ is the factor responsible for the sign of the effective pressure in the model. The geometrical density $ \rho_z $ will be associated to a positive effective pressure when $ {\cal{G}} $ satisfies the condition,

and the negative pressure occurs in the condition,

In particular, if $ {\cal{G}}<0 $ then $ p_z<0 $ is satisfied automatically.

In the phase space $ {\cal{G}} \times H $, the positive and negative sign of the effective pressure is divided by the curve $ {\cal{G}} = 12H^4 $ (see Figure 8). If we compare with the global behavior map given in Figure 7, then, if $ {\cal{G}} $ is in the cyclic region III, the effective pressure changes sign periodically, while in region II the pressure may change sign from negative to positive, and in region IV it becomes negative again, and finally, in region I the pressure is always negative.

Figure 8.  (color online) Regions in the phase space $ {\cal{G}} \times H$ where the effective pressure is positive (region with blue plus), and negative pressure (region with negative red). The solid curve $ {\cal{G}}=12H^4$ is where the pressure is zero.

Let us consider now a Friedmann's universe with cosmic dust, in this case $ p_m = 0 $. The set of differential equations are

This set of equations may be expressed as an autonomous 3D dynamical system in terms of the independent variables $ {\cal{G}} $, H and $ \rho_m $ given by (29), (30) and (31), that is,

The equilibrium points of the system obtained by setting $ \dot{{\cal{G}}} = 0 $, $ \dot{H} = 0 $ and $ \dot{\rho_m} = 0 $ are,

which are the same equilibrium points as in the case of an empty universe. The eigenvalues of the Jacobian matrix of the system, applied at the equilibrium points, have mixed signs, this implies they behave as saddle points for the integral solutions, just like in the empty space case. Besides these two equilibrium points, the solution $ H = 0 $ satisfies equations (66) through (68), with the conditions $ {\cal{G}} = 0 $ and $ \rho_m = 0 $, as it was shown in subsection IV.A.4, with $ \rho_m = -\rho_z = \dfrac{1}{2}f(0) = 0 $.

Numerical solutions to this system present behaviors like those observed in the empty space case, that is, for $ H(t) $ and σ positive, solutions in which the Hubble parameter decreases to a minimum positive or zero value and increases afterwards, or open cyclic behavior in the phase space. Since now we're dealing with a 3D dynamical system, a compilation of integral solutions in one 3D graph on the phase space becomes convoluted. A sample of two integral solutions, representing the two described behaviors, with given initial conditions, are given below:

In Figures 9a, 9b, 10a and 10b we have a case of integral solution with cyclic behavior. This shows that the cosmic dust has a dampening effect on the evolution of the system, the integral solution does not necessarily need to pass through the origin in the $ {\cal{G}} \times H $ phase plane (as seen in Figure 9a); the Hubble parameter H starts with an initial value and then decreases periodically over time, like a ball bouncing from a surface a few times, approaching a stable value (as observed in Figure 9b); the topological invariant $ {\cal{G}} $ changes its sign a few times until a point in time when it approaches a stable value (as viewed in Figure 10a), that is, it describes a Friedmann's universe that goes through epochs of accelerated and decelerated expansion, approaching a stable value with zero acceleration; and finally, the matter density decreases as expected.

Figure 9.  (color online) Numerical solution of the dynamical system in the phase space $ {\cal{G}} \times H \times \rho_m $. Initial conditions $ {\cal{G}}(0) = 1 $, $ H(0) = 1 $, $ \rho_m (0) = 15 $.

Figure 10.  (color online) Numerical solution of the dynamical system in the phase space $ {\cal{G}} \times H \times \rho_m $. Initial conditions $ {\cal{G}}(0) = 1 $, $ H(0) = 1 $, $ \rho_m (0) = 15 $.

Figures 11a, 11b, 12a and 12b show a case where the Hubble parameter has only one minimum positive value, and the topological invariant changes sign only one time.

Figure 11.  (color online) Numerical solution of the dynamical system in the phase space $ {\cal{G}} \times H \times \rho_m $. Initial conditions $ {\cal{G}}(0) = 1 $, $ H(0) = 1 $, $ \rho(0) = 1 $.

Figure 12.  (color online) Numerical solution of the dynamical system in the phase space $ {\cal{G}} \times H \times \rho_m$. Initial conditions $ {\cal{G}}(0)=1$, $ H(0)=1$, $\rho(0)=1$.

One important qualitative result here is that, given a positive value for σ, and fixed initial conditions for $ {\cal{G}}(0) $ and $ H(0) = 0 $, the initial value of the cosmic dust density $ \rho_m(0) $ will dictate which type of behavior the system will develop, that is, if the Hubble parameter will have only one minimum value, or have a cyclic behavior with epochs of acceleration and deceleration. Then $ \rho_m $ plays an important role on the evolution of the system.

The dampening effect due to the cosmic dust in the system is a consequence of its direct action on $ \dot{{\cal{G}}} $, as it can be seen in (69). When the system is about to change from $ {\cal{G}} $ positive to $ {\cal{G}} $ negative, this implies $ \dot{{\cal{G}}} < 0 $, and we see from (69) and (70) that, the presence of the cosmic dust increases the negative value of $ \dot{{\cal{G}}} $ but does not affect directly the $ \dot{H} $. The consequence of this modification to the dynamical system is that, the cosmic dust accelerates the change of sign of the topological invariant from positive to negative before the Hubble parameter can reach zero. For example, consider the cyclic solution in Figure 9a, when it is approaching the change of sign of $ {\cal{G}} $ from positive to negative, this implies the condition over the Hubble parameter,

and when the topological invariant is exactly zero we have the specific condition,

The above condition shows a mechanism for which the cosmic dust $ \rho_m $ limits the maximum value that the Hubble parameter H can have when the topological invariant $ {\cal{G}} $ change sign from positive to negative, that is, when the model change from an epoch of decelerated expansion to an epoch of accelerated expansion. As the $ \rho_m $ tends to zero (see Figure 10b), the condition (78) implies that $ H\rightarrow 0 $, which is the behavior of the empty spacetime solution observed, for example, in Figure 6.

The rank-3 tensor representation of a spin-2 field appears directly in the modified equations of $ f({\cal{G}}) $ gravity as a covariant divergent. In a homogeneous and isotropic Friedmann's universe, the trivial case for $ {\cal{G}} $ as a constant gives rise to specific conditions for $ \rho_m\geq 0 $ to be true, that is, depending on the sign of $ {\cal{G}} $, the geometric density $ \rho_z $ must satisfy certain conditions and the Hubble parameter $ H(t) $ may only vary over delimited time intervals. In other words, in some cases the topological invariant $ {\cal{G}} $ can only be constant over a finite interval of time in order for the density mass $ \rho_m $ to be strictly positive, this gives one necessary, but not a sufficient tool to filter out non-physical solutions.

The first non-trivial case we can consider, $ f({\cal{G}}) = \sigma G^2 $, gives rise to different behaviors for the Friedmann's universe. In the empty space-time case, it can describe a universe that is always in accelerated expansion, or that is decelerated and accelerates at a later phase; or describe a universe that is expanding in phases of decreasing acceleration and deceleration. When cosmic dust is added to the system, its initial value will dictate how this universe will evolve, if it will expand with increasing dampened phases of acceleration and deceleration, or if it will expand with ever increasing acceleration. The Hubble parameter H as a non zero constant is a viable solution for an empty space-time, while $ H = 0 $ is a possible solution when a constant source matter $ \rho_m $ is added.

The next step would to add a gas of photons to the system and study its influence. And to compare the numerical solutions with real observational data to constrain the value of σ and verify if it is reasonable. In our dynamical systems $ {\cal{G}} $ was unbounded, it would be worthwhile to bound $ {\cal{G}} $ to a finite interval using a Born-Infeld type function, summing or multiplying the original function, and verify how this boundary for $ {\cal{G}} $, and its variation, would affect the evolution of the system. Another point of interest would be to study perturbations to determine how stable and unstable the equilibrium points are, and in particular there is the origin of the phase space studied that may present chaotic behavior in the evolution of the system.

M. Novello thanks Fundação de Amparo e Pesquisa do Estado do Rio de Janeiro (FAPERJ) for a fellowship

Fábio H. M. dos Anjos thanks Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) fot the support