Dynamical behavior in f (T, T_G) cosmology

Since the discovery of the universe late-times acceleration, a large amount of research has been devoted to its explanation. In principle, one can follow two main directions to achieve this. The first way is to modify the content of the universe, introducing the dark energy concept, with its simpler candidates being a canonical scalar field, a phantom field, or the combination of both fields in a unified model dubbed quintom (for reviews on dark energy see [1, 2] and references therein). The second direction that one can follow is to modify the gravitational sector itself (for a review see [3] and references therein), acquiring a modified cosmological dynamics. However, note that apart from the interpretation, one can transform from one approach to the other, since the crucial issue is just the number of degrees of freedom beyond General Relativity and standard model particles (see [4] for a review on such a unified point of view). Finally, note that the above scenarios, apart from late-times implications, can be also used for the description of the inflationary stage [5].

In the majority of modified gravitational theories, one suitably extends the curvature-based Einstein–Hilbert action of General Relativity. However, an interesting class of gravitational modification arises when one modifies the action of the equivalent formulation of General Relativity, based on torsion. In particular, it is known that Einstein himself constructed the so-called 'Teleparallel Equivalent of General Relativity' (TEGR) [6–9] using the curvature-less Weitzenböck connection instead of the torsion-less Levi-Civita one. The corresponding Lagrangian, namely the torsion scalar T, is constructed by contractions of the torsion tensor, in a similar way that the usual Einstein–Hilbert Lagrangian R is constructed by contractions of the curvature (Riemann) tensor. Thus, inspired by the f(R) modifications of the Einstein–Hilbert Lagrangian [10, 11], one can construct the f(T) modified gravity by extending T to an arbitrary function [12–14]. Note that although TEGR coincides with General Relativity at the level of equations, f(T) does not coincide with f(R); i.e., they represent different modification classes. Thus, the cosmological implications of f(T) gravity are new and interesting [14–43].

However, in curvature gravity, apart from the simple f(R) modification, one can construct more complicated extensions using higher-curvature corrections such as the Gauss–Bonnet term G [44–46] or functions of it [47, 48], Lovelock combinations [49–51], and Weyl combinations [52–54]. Inspired by these, in recent work [55], the $f(T,{{T}_{G}})$ gravitational modification was constructed. It is based on the old quadratic torsion scalar T as well as on the new quartic torsion scalar T_G, i.e., the teleparallel equivalent of the Gauss–Bonnet term. Obviously, $f(T,{{T}_{G}})$ theories cannot arise from the f(T) ones; additionally, they are different from the $f(R,G)$ class of curvature modified gravity. Thus, $f(T,{{T}_{G}})$ is a novel class of gravitational modification.

The cosmological applications of $f(T,{{T}_{G}})$ gravity prove to be very interesting [55]. Therefore, it is both interesting and necessary to perform a dynamical analysis, examining in a systematic way the allowed cosmological behaviors, focusing on the late-times stable solutions. The phase-space and stability analysis is a very powerful tool, since it reveals the global features of a given cosmological scenario independently of the initial conditions and the specific evolution of the universe. In the present investigation, we perform such a detailed phase-space analysis, and we extract the late-times, asymptotic solutions, calculating also the corresponding observable quantities such as the deceleration parameter, the effective dark energy equation-of-state parameter, and the various density parameters.

The plan of the work is the following: in section 2 we briefly review the scenario of $f(T,{{T}_{G}})$ gravity, and in section 3 we present its application in cosmology. In section 4 we perform the detailed dynamical analysis for the simplest non-trivial model of $f(T,{{T}_{G}})$ gravity. In section 5 we discuss the cosmological implications and the physical behavior of the scenario. Finally, in section 6 we summarize our results.

In this section we briefly review the $f(T,{{T}_{G}})$ gravitational modification following [55]. In the whole manuscript we use the following notation: Greek indices run over the coordinate space-time, while Latin indices run over its tangent space.

In this framework, the dynamical variable is the vierbein field ${{e}_{a}}({{x}^{\mu }})$. In terms of coordinates, it can be expressed in components as ${{e}_{a}}=e_{a}^{\mu }{{\partial }_{\mu }}$, while the dual vierbein is defined as ${{e}^{a}}=e_{\mu }^{a}{\rm d}{{x}^{\mu }}$. Concerning the other field, i.e., the connection 1-forms $\omega _{}^{\;\;\;ba}({{x}^{\mu }})$ (which defines the parallel transportation), one uses the Weitzenböck one, which in all coordinate frames is defined as:

$$\begin{eqnarray}&&\omega _{\mu \nu }^{\lambda }=e_{a}^{\lambda }e_{\mu ,\nu }^{a}.\end{eqnarray} \tag{ 2.1 }$$

Due to its inhomogeneous transformation law, it has tangent-space components $\omega _{bc}^{a}=0$, assuring the property of vanishing non-metricity. Additionally, for an orthonormal vierbein, the metric tensor is given by the relation:

$$\begin{eqnarray}&&{{g}_{\mu \nu }}={{\eta }_{ab}}\;e_{\mu }^{a}\;e_{\nu }^{b},\end{eqnarray} \tag{ 2.2 }$$

where ${{\eta }_{ab}}={\rm diag}(-1,1,1,1)$ and indices $a,b,...$ are raised/lowered with η_ab.

One can now define the torsion tensor as:

$$\begin{eqnarray}&&T_{\mu \nu }^{\lambda }=e_{a}^{\lambda }\left( {{\partial }_{\nu }}e_{\;\mu }^{a}-{{\partial }_{\mu }}e_{\;\nu }^{a} \right),\end{eqnarray} \tag{ 2.3 }$$

while the Riemann tensor is zero by construction, due to the teleparallelism condition which is imposed with the use of the Weitzenböck connection. Moreover, the contorsion tensor, which equals the difference between the Weitzenböck and Levi-Civita connections, is defined as:

$$\begin{eqnarray}&&\mathcal{K}_{\ \ \ \ \rho }^{\mu \nu }=-\frac{1}{2}\left( T_{\ \ \ \ \rho }^{\mu \nu }-T_{\ \ \ \ \rho }^{\nu \mu }-T_{\rho }^{\ \ \mu \nu } \right).\end{eqnarray} \tag{ 2.4 }$$

Since in this formulation all the information concerning the gravitational field is included in the torsion tensor $T_{}^{\;\;\;\mu \nu \lambda }$, one can use it in order to construct torsion invariants. The simplest invariants that one can build are quadratic in the torsion tensor. In particular, the combination

$$\begin{eqnarray}&&T=\frac{1}{4}{{T}^{\mu \nu \lambda }}{{T}_{\mu \nu \lambda }}+\frac{1}{2}{{T}^{\mu \nu \lambda }}{{T}_{\lambda \nu \mu }}-T_{\nu }^{\nu \mu }T_{\lambda \mu }^{\lambda },\end{eqnarray} \tag{ 2.5 }$$

which can in general be defined in an arbitrary dimension D, is the 'torsion scalar', and if it is used as a Lagrangian and varied in terms of the vierbein, it gives rise to the Einstein field equations. That is why the gravitational theory characterized by the action:

$$\begin{eqnarray}&&S=-\frac{1}{2{{\kappa }^{2}}}\int {{{\rm d}}^{4}}x\;e\;T,\end{eqnarray} \tag{ 2.6 }$$

with $e={\rm det} (e_{\mu }^{a})=\sqrt{|g|}$ and ${{\kappa }^{2}}\equiv 8\pi G$ the gravitational constant, is called TEGR. In these lines, one can be based on T in order to construct modified gravitational theories extending the TEGR action to [14–43]:

$$\begin{eqnarray}&&S=\frac{1}{2{{\kappa }^{2}}}\int {{{\rm d}}^{4}}x\;e\;f(T).\end{eqnarray} \tag{ 2.7 }$$

We stress here that although the field equations of TEGR are identical with those of General Relativity, f(T) modification gives rise to different equations than f(R) modification.

However, one can use the torsion tensor in order to construct higher-order torsion invariants, in a similar way that one uses the Riemann tensor in order to construct higher-order curvature invariants. In particular, in [55] the invariant:

$$\begin{eqnarray}\begin{array}{rcl} {{T}_{G}} & = & \left( \mathcal{K}_{\varphi \pi }^{\kappa }\mathcal{K}_{\,\rho }^{\varphi \lambda }\mathcal{K}_{\,\chi \sigma }^{\mu }\mathcal{K}_{\tau }^{\chi \nu }-2\mathcal{K}_{\,\pi }^{\kappa \lambda }\mathcal{K}_{\,\,\varphi \rho }^{\mu }\mathcal{K}_{\chi \sigma }^{\varphi }\mathcal{K}_{\,\tau }^{\chi \nu } \right. \\ {} & {} & \left. +2\mathcal{K}_{\,\pi }^{\kappa \lambda }\mathcal{K}_{\,\varphi \rho }^{\mu }\mathcal{K}_{\,\chi }^{\varphi \nu }\mathcal{K}_{\,\sigma \tau }^{\chi }+2\mathcal{K}_{\,\pi }^{\kappa \lambda }\mathcal{K}_{\,\varphi \rho }^{\mu }\mathcal{K}_{\sigma ,\tau }^{\varphi \nu } \right)\delta _{\kappa \lambda \mu \nu }^{\pi \rho \sigma \tau } \\ \end{array}\end{eqnarray} \tag{ 2.8 }$$

was constructed in an arbitrary dimension D, where the generalized $\delta _{\kappa \lambda \mu \nu }^{\pi \rho \sigma \tau }$ is the determinant of the Kronecker deltas. This invariant is just the Teleparallel Equivalent of the Gauss–Bonnet combination $G={{R}^{2}}-4{{R}_{\mu \nu }}{{R}^{\mu \nu }}+{{R}_{\mu \nu \kappa \lambda }}{{R}^{\mu \nu \kappa \lambda }}$, and in four dimensions it reduces to a topological term. Thus, inspired by the f(G) extensions of General Relativity [47, 48], one can consider general functions $f({{T}_{G}})$ in the action too.

Taking the above into account, one can propose a new class of gravitational modifications as [55]:

$$\begin{eqnarray}&&S=\frac{1}{2{{\kappa }^{2}}}\int {{{\rm d}}^{4}}x\;e\;f(T,{{T}_{G}}),\end{eqnarray} \tag{ 2.9 }$$

which is also valid in higher dimensions. Since T_G is quartic in the torsion tensor, $f(T,{{T}_{G}})$ gravity is more general than the f(T) class. Additionally, $f(T,{{T}_{G}})$ gravity is obviously different from the $f(R,G)$ one [47, 48, 56–58]. Note that the usual Einstein–Gauss–Bonnet theory for $D\gt 4$ arises in the special case $f(T,{{T}_{G}})=-T+\alpha {{T}_{G}}$ (with α the Gauss–Bonnet coupling), while TEGR (that is GR) is obtained for $f(T,{{T}_{G}})=-T$.

In order to investigate the cosmological implication of the above action (2.9), we consider a spatially flat cosmological ansatz:

$$\begin{eqnarray}&&{\rm d}{{s}^{2}}=-{\rm d}{{t}^{2}}+{{a}^{2}}(t){{\delta }_{\hat{i}\hat{j}}}{\rm d}{{x}^{{\hat{i}}}}{\rm d}{{x}^{{\hat{j}}}},\end{eqnarray} \tag{ 3.1 }$$

where a(t) is the scale factor. This metric arises from the diagonal vierbein:

$$\begin{eqnarray}&&e_{\mu }^{a}={\rm diag}(1,a(t),a(t),a(t))\end{eqnarray} \tag{ 3.2 }$$

through (2.2), while the dual vierbein is $e_{a}^{\mu }={\rm diag}(1,{{a}^{-1}}(t),{{a}^{-1}}(t),{{a}^{-1}}(t))$, and its determinant $e=a{{(t)}^{3}}$. Thus, inserting the vierbein (3.2) into relations (2.5) and (2.8), we find:

$$\begin{eqnarray}&&T=6\;{{H}^{2}}\end{eqnarray} \tag{ 3.3 }$$

$$\begin{eqnarray}&&{{T}_{G}}=24\;{{H}^{2}}\left( \dot{H}+{{H}^{2}} \right),\end{eqnarray} \tag{ 3.4 }$$

where $H=\frac{{\dot{a}}}{a}$ is the Hubble parameter and dots denote differentiation with respect to t.

Finally, in order to acquire a realistic cosmology, we additionally consider a matter action S_m, corresponding to an energy-momentum tensor ${{\Theta }^{\mu \nu }}$, focusing on the case of a perfect fluid of energy density ρ_m and pressure p_m.

As shown in [55], the variation of the total action $S+{{S}_{m}}$ gives, in the case of FRW geometry, the following Friedmann equations:

$$\begin{eqnarray}&&f-12\;{{H}^{2}}{{f}_{T}}-{{T}_{G}}{{f}_{{{T}_{G}}}}+24\;{{H}^{3}}\mathop{{{f}_{{{T}_{G}}}}}\limits^{.}=2{{\kappa }^{2}}{{\rho }_{m}}\end{eqnarray} \tag{ 3.5 }$$

$$\begin{eqnarray}&&f-4\left( 3\;{{H}^{2}}+\dot{H} \right){{f}_{T}}-4\;H\mathop{{{f}_{T}}}\limits^{.}-{{T}_{G}}{{f}_{{{T}_{G}}}}+\frac{2}{3\;H}{{T}_{G}}\mathop{{{f}_{{{T}_{G}}}}}\limits^{.}+8\;{{H}^{2}}\ddot{{{f}_{{{T}_{G}}}}}=-2{{\kappa }^{2}}{{p}_{m}},\end{eqnarray} \tag{ 3.6 }$$

where $\mathop{{{f}_{T}}}\limits^{.}={{f}_{TT}}\dot{T}+{{f}_{T{{T}_{G}}}}{{\dot{T}}_{G}}$, $\mathop{{{f}_{{{T}_{G}}}}}\limits^{.}={{f}_{T{{T}_{G}}}}\dot{T}+{{f}_{{{T}_{G}}{{T}_{G}}}}{{\dot{T}}_{G}}$, $\ddot{{{f}_{{{T}_{G}}}}}={{f}_{TT{{T}_{G}}}}{{\dot{T}}^{2}}+2{{f}_{T{{T}_{G}}{{T}_{G}}}}\dot{T}{{\dot{T}}_{G}}+{{f}_{{{T}_{G}}{{T}_{G}}{{T}_{G}}}}\dot{T}_{G}^{2}+{{f}_{T{{T}_{G}}}}\ddot{T}+{{f}_{{{T}_{G}}{{T}_{G}}}}{{\ddot{T}}_{G}}$, with f_TT, ${{f}_{T{{T}_{G}}}}$,..., denoting multiple partial differentiations of f with respect to T, T_G. Here, the involved time-derivatives of $\dot{T}$, $\ddot{T}$, ${{\dot{T}}_{G}}$, ${{\ddot{T}}_{G}}$ are straightforwardly obtained using (3.3) and (3.4).

Therefore, we can rewrite the Friedmann equations (3.5) and (3.6) in the usual form:

$$\begin{eqnarray}&&{{H}^{2}}=\frac{{{\kappa }^{2}}}{3}\left( {{\rho }_{m}}+{{\rho }_{{\rm DE}}} \right)\end{eqnarray} \tag{ 3.7 }$$

$$\begin{eqnarray}&&\dot{H}=-\frac{{{\kappa }^{2}}}{2}\left( {{\rho }_{m}}+{{p}_{m}}+{{\rho }_{{\rm DE}}}+{{p}_{{\rm DE}}} \right),\end{eqnarray} \tag{ 3.8 }$$

defining the energy density and pressure of the effective dark energy sector as:

$$\begin{eqnarray}&&{{\rho }_{{\rm DE}}}\equiv \frac{1}{2{{\kappa }^{2}}}\left( 6\;{{H}^{2}}-f+12\;{{H}^{2}}{{f}_{T}}+{{T}_{G}}{{f}_{{{T}_{G}}}}-24\;{{H}^{3}}\mathop{{{f}_{{{T}_{G}}}}}\limits^{.} \right)\end{eqnarray} \tag{ 3.9 }$$

$$\begin{eqnarray}&&\begin{array}{l} {{p}_{{\rm DE}}}\equiv \frac{1}{2{{\kappa }^{2}}}\left[ -2(2\dot{H}+3\;{{H}^{2}})+f-4\left( \dot{H}+3\;{{H}^{2}} \right){{f}_{T}}-4\;H\mathop{{{f}_{T}}}\limits^{.} \right. \\ \qquad \left. -{{T}_{G}}{{f}_{{{T}_{G}}}}+\frac{2}{3\;H}{{T}_{G}}\mathop{{{f}_{{{T}_{G}}}}}\limits^{.}+8\;{{H}^{2}}\ddot{{{f}_{{{T}_{G}}}}} \right]. \\ \end{array}\end{eqnarray} \tag{ 3.10 }$$

The standard matter ρ_m is conserved independently; i.e., ${{\dot{\rho }}_{m}}+3\;H({{\rho }_{m}}+{{p}_{m}})=0$. One can easily verify that the dark energy density and pressure satisfy the usual evolution equation:

$$\begin{eqnarray}&&{{\dot{\rho }}_{{\rm DE}}}+3\;H({{\rho }_{{\rm DE}}}+{{p}_{{\rm DE}}})=0,\end{eqnarray} \tag{ 3.11 }$$

and we can also define the dark energy equation-of-state parameter as usual:

$$\begin{eqnarray}&&{{w}_{{\rm DE}}}\equiv \frac{{{p}_{{\rm DE}}}}{{{\rho }_{{\rm DE}}}}.\end{eqnarray} \tag{ 3.12 }$$

In order to perform the stability analysis of a given cosmological scenario, one first transforms it to its autonomous form ${\bf X}^{\prime} ={\bf f}({\bf X})$ [59–65], where ${\bf X}$ are some auxiliary variables presented as a column vector, and primes denote derivatives with respect to $N={\rm ln} a$. Then, one extracts the critical points ${{{\bf X}}_{c}}$ by imposing the condition ${\bf X}^{\prime} ={\bf 0}$. In order to determine their stability properties, one expands around them with ${\bf U}$ the column vector of the perturbations of the variables. Therefore, for each critical point, the perturbation equations are expanded to first order as ${\bf U}^{\prime} ={\bf Q}\cdot {\bf U}$, with the matrix ${\bf Q}$ containing the coefficients of the perturbation equations. The eigenvalues of ${\bf Q}$ determine the type and stability of the specific critical point.

In order to perform the above analysis, we need to specify the $f(T,{{T}_{G}})$ form. In usual f(T) gravity, one starts adding corrections of T-powers. However, in the scenario at hand, since T_G contains quartic torsion terms, it is of the same order with T². Therefore, T and $\sqrt{{{T}^{2}}+{{\alpha }_{2}}{{T}_{G}}}$ are of the same order, and thus, one should use both in a modified theory. Hence, the simplest non-trivial model, which does not introduce a new mass scale into the problem and differs from General Relativity, is the one based on:

$$\begin{eqnarray}&&f(T,{{T}_{G}})=-T+{{\alpha }_{1}}\sqrt{{{T}^{2}}+{{\alpha }_{2}}{{T}_{G}}}.\end{eqnarray} \tag{ 4.1 }$$

The couplings ${{\alpha }_{1}},{{\alpha }_{2}}$ are dimensionless and the model is expected to play an important role at late times. Indeed, this model, although simple, can lead to interesting cosmological behavior, revealing the advantages, the capabilities, and the new features of $f(T,{{T}_{G}})$ cosmology. We mention here that when ${{\alpha }_{2}}=0$, this scenario reduces to TEGR, i.e., to General Relativity, with just a rescaled Newtonʼs constant, whose dynamical analysis has been performed in detail in the literature [61–63]. Thus, in the following we restrict our analysis to the case ${{\alpha }_{2}}\ne 0$.

In this case, the cosmological equations are the Friedmann equations (3.7) and (3.8), with the effective dark energy density and pressure (3.9) and (3.10) becoming:

$$\begin{eqnarray}&&{{\kappa }^{2}}{{\rho }_{{\rm DE}}}=\frac{\sqrt{3}{{\alpha }_{1}}{{H}^{2}}\left\{ \alpha _{2}^{2}\ddot{H}+9{{\alpha }_{2}}H\dot{H}+\left[ (3-2{{\alpha }_{2}}){{\alpha }_{2}}+9 \right]{{H}^{3}} \right\}}{{{D}^{3/2}}}\end{eqnarray} \tag{ 4.2 }$$

$$\begin{eqnarray}\begin{array}{rcl} {{\kappa }^{2}}{{p}_{{\rm DE}}} & = & \frac{\begin{array}{ccccccccccccccc} {{\alpha }_{1}}\left\{ (2{{\alpha }_{2}}+3)\left[ {{\alpha }_{2}}(10{{\alpha }_{2}}-51)-18 \right]{{H}^{4}} \right.\left. +{{\alpha }_{2}}\left[ 4{{\alpha }_{2}}(5{{\alpha }_{2}}-21)-90 \right]{{H}^{2}}\dot{H}-54\alpha _{2}^{2}{{{\dot{H}}}^{2}} \right\}H\dot{H} \\ \end{array}}{\sqrt{3}{{D}^{5/2}}} \\ {} & {} & -\frac{{{\alpha }_{1}}\alpha _{2}^{2}H\ddot{H}}{\sqrt{3}{{D}^{3/2}}}+\frac{\sqrt{3}{{\alpha }_{1}}\alpha _{2}^{3}\;H{{\ddot{H}}^{2}}}{{{D}^{5/2}}}-\frac{2{{\alpha }_{1}}\alpha _{2}^{2}\ddot{H}\left[ 2({{\alpha }_{2}}-3){{H}^{2}}\dot{H}+2{{\alpha }_{2}}{{{\dot{H}}}^{2}}+(6{{\alpha }_{2}}+9){{H}^{4}} \right]}{\sqrt{3}{{D}^{5/2}}} \\ {} & {} & +\frac{\sqrt{3}{{\alpha }_{1}}({{\alpha }_{2}}-3){{\left( 2{{\alpha }_{2}}+3 \right)}^{2}}\;{{H}^{7}}}{{{D}^{5/2}}}, \\ \end{array}\end{eqnarray} \tag{ 4.3 }$$

where $D=3\;{{H}^{2}}+2{{\alpha }_{2}}(\dot{H}+{{H}^{2}})$. In order to perform the dynamical analysis of this cosmological scenario, we introduce the following auxiliary variables:

$$\begin{eqnarray}&&x=\sqrt{\frac{D}{3\;{{H}^{2}}}}=\sqrt{1+\frac{2{{\alpha }_{2}}}{3}\left( 1+\frac{{\dot{H}}}{{{H}^{2}}} \right)}\end{eqnarray} \tag{ 4.4 }$$

$$\begin{eqnarray}&&{{\Omega }_{m}}=\frac{{{\kappa }^{2}}{{\rho }_{m}}}{3\;{{H}^{2}}}.\end{eqnarray} \tag{ 4.5 }$$

Thus, the cosmological system is transformed to the following autonomous form:

$$\begin{eqnarray}&&x^{\prime} =-\frac{x\left[ 3{{\alpha }_{1}}{{x}^{2}}-6(1-{{\Omega }_{m}})x+{{\alpha }_{1}}(3-4{{\alpha }_{2}}) \right]}{2{{\alpha }_{1}}{{\alpha }_{2}}}\end{eqnarray} \tag{ 4.6 }$$

$$\begin{eqnarray}&&\Omega _{m}^{\prime }=-\frac{{{\Omega }_{m}}\left( 3{{x}^{2}}+{{\alpha }_{2}}+3{{\alpha }_{2}}{{w}_{m}}-3 \right)}{{{\alpha }_{2}}}\;,\end{eqnarray} \tag{ 4.7 }$$

where primes denote differentiation with respect to the new time variable N, so $f^{\prime} ={{H}^{-1}}\dot{f}$. The above dynamical system is defined in the phase space $\left\{ (x,{{\Omega }_{m}})|x\in [0,\infty ),{{\Omega }_{m}}\in [0,\infty ] \right\}$.

One can now express the various observables in terms of the above auxiliary variables Ω_m and x. (Note that Ω_m is an observable itself, i.e., the matter density parameter.) In particular, the deceleration parameter $q\equiv -1-\dot{H}/{{H}^{2}}$ is given by:

$$\begin{eqnarray}&&q=\frac{3\left( 1-{{x}^{2}} \right)}{2{{\alpha }_{2}}}.\end{eqnarray} \tag{ 4.8 }$$

Similarly, the dark energy density parameter straightaway reads:

$$\begin{eqnarray}&&{{\Omega }_{{\rm DE}}}\equiv \frac{{{\kappa }^{2}}{{\rho }_{{\rm DE}}}}{3\;{{H}^{2}}}=1-{{\Omega }_{m}}.\end{eqnarray} \tag{ 4.9 }$$

The dark energy equation-of-state parameter w_DE is given by the relation $2q=1+3({{w}_{m}}{{\Omega }_{m}}+{{w}_{{\rm DE}}}{{\Omega }_{{\rm DE}}})$, and therefore:

$$\begin{eqnarray}&&{{w}_{{\rm DE}}}=\frac{3{{x}^{2}}+{{\alpha }_{2}}+3{{\alpha }_{2}}{{w}_{m}}{{\Omega }_{m}}-3}{3{{\alpha }_{2}}({{\Omega }_{m}}-1)},\end{eqnarray} \tag{ 4.10 }$$

where ${{w}_{m}}\equiv \frac{{{p}_{m}}}{{{\rho }_{m}}}$ is the matter equation-of-state parameter. In the following, without loss of generality, we assume dust matter (w_m = 0), but the extension to general w_m is straightforward.

4.1. Finite phase space analysis

We now proceed to the detailed phase-space analysis. The real and physically interesting (i.e., corresponding to an expanding universe) critical points of the autonomous system (4.6)–(4.7), obtained by setting the left hand sides of these equations to zero, are presented in table 1 . In the same table we provide their existence conditions. Their stability is extracted by examining the sign of the real part of the eigenvalues of the 2 × 2 matrix ${\bf Q}$ of the corresponding linearized perturbation equations. This procedure is shown in appendix A, and in table 1 we summarize the stability results. Furthermore, for each critical point, we calculate the values of the deceleration parameter q and the dark energy equation-of-state parameter w_DE given by (4.8) and (4.10), and we present the results in table 2. Finally, in the same table we summarize the physical description of the solutions, which we analyze in the next section.

Table 1.  The real and physically interesting critical points of the autonomous system (4.6)–(4.7). Existence and stability conditions. We use the notations ${{\Omega }_{m1}}=\frac{{{\alpha }_{1}}\sqrt{9-3{{\alpha }_{2}}}(6-5{{\alpha }_{2}})+6({{\alpha }_{2}}-3)}{6({{\alpha }_{2}}-3)}$, ${{x}_{2}}=\frac{3-\sqrt{3\alpha _{1}^{2}(4{{\alpha }_{2}}-3)+9}}{3{{\alpha }_{1}}}$, and ${{x}_{3}}=\frac{3+\sqrt{3\alpha _{1}^{2}(4{{\alpha }_{2}}-3)+9}}{3{{\alpha }_{1}}}$.

Cr. P.	x	Ωm	Existence	Stability
P1	$\sqrt{1-\frac{{{\alpha }_{2}}}{3}}$	Ωm1	$\frac{6}{5}\lt {{\alpha }_{2}}\lt 3,{{\alpha }_{1}}\geqslant -2\sqrt{\frac{3(3-{{\alpha }_{2}})}{{{(-6+5{{\alpha }_{2}})}^{2}}}}$ or	Stable spiral for ${{\alpha }_{2}}\lt 3$ and
			${{\alpha }_{2}}=\frac{6}{5}$ or	$-32\sqrt{3}\sqrt{\frac{{{(3-{{\alpha }_{2}})}^{3}}}{{{\left( 71\alpha _{2}^{2}-336{{\alpha }_{2}}+288 \right)}^{2}}}}\lt {{\alpha }_{1}}\lt 0$ or
			${{\alpha }_{2}}\lt \frac{6}{5},{{\alpha }_{1}}\leqslant 2\sqrt{\frac{3(3-{{\alpha }_{2}})}{{{(-6+5{{\alpha }_{2}})}^{2}}}}$	${{\alpha }_{1}}\lt 0,{{\alpha }_{2}}\leqslant \frac{1}{71}(168-36\sqrt{6})\approx 1.124$.
				Saddle otherwise (hyperbolic cases).
P2	x2	0	${{\alpha }_{2}}\lt \frac{3}{4},\;0\lt {{\alpha }_{1}}\leqslant \sqrt{\frac{3}{3-4{{\alpha }_{2}}}}$ or	Stable node for ${{\alpha }_{2}}\lt 0,\;0\lt {{\alpha }_{1}}\lt 2\sqrt{\frac{3(3-{{\alpha }_{2}})}{{{(5{{\alpha }_{2}}-6)}^{2}}}}$
			${{\alpha }_{1}}\ne 0,\;{{\alpha }_{2}}=\frac{3}{4}$ or	or $\frac{6}{5}\lt {{\alpha }_{2}}\leqslant 3,\;{{\alpha }_{1}}\lt -2\sqrt{\frac{3(3-{{\alpha }_{2}})}{{{(5{{\alpha }_{2}}-6)}^{2}}}}$
			${{\alpha }_{2}}\gt \frac{3}{4},\;{{\alpha }_{1}}\lt 0$	or ${{\alpha }_{2}}\gt 3,\;{{\alpha }_{1}}\lt 0$.
				Unstable node for $0\lt {{\alpha }_{2}}\lt \frac{3}{4},\;0\lt {{\alpha }_{1}}\lt \sqrt{\frac{3}{3-4{{\alpha }_{2}}}}.$
				Saddle otherwise (hyperbolic cases).
P3	x3	0	${{\alpha }_{2}}\lt \frac{3}{4},\;0\lt {{\alpha }_{1}}\leqslant \sqrt{\frac{3}{3-4{{\alpha }_{2}}}}$ or	Stable node for ${{\alpha }_{1}}\gt 0,{{\alpha }_{2}}\geqslant \frac{6}{5}.$
			${{\alpha }_{2}}\geqslant \frac{3}{4},\;{{\alpha }_{1}}\gt 0$	Unstable node for ${{\alpha }_{2}}\lt 0,\;0\lt {{\alpha }_{1}}\lt \frac{\sqrt{3}}{\sqrt{3-4{{\alpha }_{2}}}}.$
				Saddle otherwise (hyperbolic cases).
P4	0	0	Always	Unstable node for $\frac{3}{4}\lt {{\alpha }_{2}}\lt 3$.
				Saddle otherwise (hyperbolic cases).

Table 2.  The real and physically interesting critical points of the autonomous system (4.6)–(4.7), and the corresponding values of the dark energy density parameter ${{\Omega }_{{\rm DE}}}$, the deceleration parameter q, and the dark energy equation-of-state parameter w_DE. We use the notation ${{\Omega }_{m1}}=\frac{{{\alpha }_{1}}\sqrt{9-3{{\alpha }_{2}}}(6-5{{\alpha }_{2}})+6({{\alpha }_{2}}-3)}{6({{\alpha }_{2}}-3)}$, ${{q}_{2}}=\frac{\sqrt{3\alpha _{1}^{2}(4{{\alpha }_{2}}-3)+9}}{\alpha _{1}^{2}{{\alpha }_{2}}}-\frac{3}{\alpha _{1}^{2}{{\alpha }_{2}}}+\frac{3}{{{\alpha }_{2}}}-2$, ${{q}_{3}}=-\frac{\sqrt{3\alpha _{1}^{2}(4{{\alpha }_{2}}-3)+9}}{\alpha _{1}^{2}{{\alpha }_{2}}}-\frac{3}{\alpha _{1}^{2}{{\alpha }_{2}}}+\frac{3}{{{\alpha }_{2}}}-2$, ${{w}_{{\rm DE}2}}=\frac{2\sqrt{3\alpha _{1}^{2}(4{{\alpha }_{2}}-3)+9}}{3\alpha _{1}^{2}{{\alpha }_{2}}}-\frac{2}{\alpha _{1}^{2}{{\alpha }_{2}}}+\frac{2}{{{\alpha }_{2}}}-\frac{5}{3}$, and ${{w}_{{\rm DE}3}}=-\frac{2\sqrt{3\alpha _{1}^{2}(4{{\alpha }_{2}}-3)+9}}{3\alpha _{1}^{2}{{\alpha }_{2}}}-\frac{2}{\alpha _{1}^{2}{{\alpha }_{2}}}+\frac{2}{{{\alpha }_{2}}}-\frac{5}{3}$. In the last column we summarize their physical description.

Cr. P.	${{\Omega }_{{\rm DE}}}$	q	wDE	Properties of solutions
P1	$1-{{\Omega }_{m1}}$	$\frac{1}{2}$	0	Dark Energy–Dark Matter scaling solution
P2	1	q2	wDE2	Decelerating solution for
				${{\alpha }_{2}}\lt 0,\ \frac{3}{3-2{{\alpha }_{2}}}\lt {{\alpha }_{1}}\leqslant \sqrt{\frac{3}{3-4{{\alpha }_{2}}}}$ or
				$0\lt {{\alpha }_{2}}\lt \frac{3}{4},\ 0\lt {{\alpha }_{1}}\leqslant \sqrt{\frac{3}{3-4{{\alpha }_{2}}}}$ or
				${{\alpha }_{1}}\ne 0,\ {{\alpha }_{2}}=\frac{3}{4}$ or
				$\frac{3}{4}\lt {{\alpha }_{2}}\leqslant \frac{3}{2},\ {{\alpha }_{1}}\lt 0$ or ${{\alpha }_{2}}\gt \frac{3}{2},\ \frac{3}{3-2{{\alpha }_{2}}}\lt {{\alpha }_{1}}\lt 0$.
				Quintessence solution for
				${{\alpha }_{2}}\leqslant -\frac{3}{2},\ \;0\lt {{\alpha }_{1}}\lt \frac{3}{3-2{{\alpha }_{2}}}$ or
				$-\frac{3}{2}\lt {{\alpha }_{2}}\lt 0,\ -\sqrt{\frac{3(2{{\alpha }_{2}}+3)}{{{({{\alpha }_{2}}-3)}^{2}}}}\lt {{\alpha }_{1}}\lt \frac{3}{3-2{{\alpha }_{2}}}$ or
				$\frac{3}{2}\lt {{\alpha }_{2}}\leqslant 3,\ {{\alpha }_{1}}\lt \frac{3}{3-2{{\alpha }_{2}}}$ or
				${{\alpha }_{2}}\gt 3,\ -\sqrt{\frac{3(2{{\alpha }_{2}}+3)}{{{({{\alpha }_{2}}-3)}^{2}}}}\lt {{\alpha }_{1}}\lt \frac{3}{3-2{{\alpha }_{2}}}$.
				De Sitter solution for
				$-\frac{3}{2}\lt {{\alpha }_{2}}\lt 0,\ {{\alpha }_{1}}=\sqrt{\frac{3(2{{\alpha }_{2}}+3)}{{{({{\alpha }_{2}}-3)}^{2}}}}$ or ${{\alpha }_{2}}\gt 3,\ {{\alpha }_{1}}=-\sqrt{\frac{3(2{{\alpha }_{2}}+3)}{{{({{\alpha }_{2}}-3)}^{2}}}}$.
				Phantom solution for
				$-\frac{3}{2}\lt {{\alpha }_{2}}\lt 0,\ 0\lt {{\alpha }_{1}}\lt \sqrt{\frac{3(2{{\alpha }_{2}}+3)}{{{({{\alpha }_{2}}-3)}^{2}}}}$ or ${{\alpha }_{2}}\gt 3,\ {{\alpha }_{1}}\lt -\sqrt{\frac{3(2{{\alpha }_{2}}+3)}{{{({{\alpha }_{2}}-3)}^{2}}}}$.
P3	1	q3	wDE3	Decelerating solution for
				${{\alpha }_{2}}\lt 0,\ 0\lt {{\alpha }_{1}}\leqslant \sqrt{\frac{3}{3-4{{\alpha }_{2}}}}$ or
				$0\lt {{\alpha }_{2}}\lt \frac{3}{4},\ \frac{3}{3-2{{\alpha }_{2}}}\lt {{\alpha }_{1}}\leqslant \sqrt{\frac{3}{3-4{{\alpha }_{2}}}}$ or
				$\frac{3}{4}\leqslant {{\alpha }_{2}}\lt \frac{3}{2},\ {{\alpha }_{1}}\gt \frac{3}{3-2{{\alpha }_{2}}}$.
				Quintessence solution for
				$0\lt {{\alpha }_{2}}\lt \frac{3}{2},\ \sqrt{\frac{3(2{{\alpha }_{2}}+3)}{{{({{\alpha }_{2}}-3)}^{2}}}}\lt {{\alpha }_{1}}\lt -\frac{3}{2{{\alpha }_{2}}-3}$ or
				$\frac{3}{2}\leqslant {{\alpha }_{2}}\lt 3,\ {{\alpha }_{1}}\gt \sqrt{\frac{3(2{{\alpha }_{2}}+3)}{{{({{\alpha }_{2}}-3)}^{2}}}}$.
				De Sitter solution for $0\lt {{\alpha }_{2}}\lt 3,\ {{\alpha }_{1}}=\sqrt{\frac{3(2{{\alpha }_{2}}+3)}{{{({{\alpha }_{2}}-3)}^{2}}}}$.
				Phantom solution for
				$0\lt {{\alpha }_{2}}\lt 3,\ 0\lt {{\alpha }_{1}}\lt \sqrt{\frac{3(2{{\alpha }_{2}}+3)}{{{({{\alpha }_{2}}-3)}^{2}}}}$ or ${{\alpha }_{2}}\geqslant 3,\ {{\alpha }_{1}}\gt 0$.
P4	1	$\frac{3}{2{{\alpha }_{2}}}$	$\frac{1}{{{\alpha }_{2}}}-\frac{1}{3}$	Decelerating solution for ${{\alpha }_{2}}\gt 0$.
				Quintessence DE dominated solution for ${{\alpha }_{2}}\lt -\frac{3}{2}$.
				De Sitter solution for ${{\alpha }_{2}}=-\frac{3}{2}$.
				Phantom solution for $-\frac{3}{2}\lt {{\alpha }_{2}}\lt 0$.

4.2. Phase space analysis at infinity

Due to the fact that the dynamical system (4.6)–(4.7) is non-compact, there could be non-trivial dynamical features could exist in the asymptotic regime too. Therefore, in order to complete the phase space analysis, we must extend our investigation with the analysis at infinity using the Poincaré projection method [66, 67].

We introduce the new coordinates $(r,\theta )$, defined by:

$$\begin{eqnarray}&&x=\frac{r}{1-r}{\rm cos} \theta \end{eqnarray} \tag{ 4.11 }$$

$$\begin{eqnarray}&&{{\Omega }_{m}}=\frac{r}{1-r}{\rm sin} \theta ,\end{eqnarray} \tag{ 4.12 }$$

with $\theta \in \left[ 0,\frac{\pi }{2} \right]$ and $r\in \left[ 0,1 \right)$. Thus, the critical points at infinity, i.e., $x\to +\infty$ or ${{\Omega }_{m}}\to +\infty$ (i.e., ${{R}^{2}}\equiv {{x}^{2}}+\Omega _{m}^{2}\to +\infty$), correspond to $r\to {{1}^{-}}$. Moreover, the region of the plane $(r,\theta )$ that corresponds to $0\leqslant x,\;0\leqslant {{\Omega }_{m}}\leqslant 1$ is given by:

$$\begin{eqnarray}&&\begin{array}{l} \left\{ (r,\theta ):0\leqslant r\leqslant \frac{1}{2},0\leqslant \theta \leqslant \frac{\pi }{2} \right\} \\ \quad \times \cup \left\{ (r,\theta ):\frac{1}{2}\lt r\lt 1,0\leqslant \theta \leqslant {\rm arcsin} \left( \frac{1-r}{r} \right) \right\}. \\ \end{array}\end{eqnarray} \tag{ 4.13 }$$

Using relations (4.11) and (4.12) and substituting into (4.8) and (4.10), we obtain the deceleration and equation-of-state parameters as a function of the new variables, namely:

$$\begin{eqnarray}&&q=\frac{3\left( 1-2r+{{r}^{2}}{{{\rm sin} }^{2}}\theta \right)}{2{{\alpha }_{2}}{{(1-r)}^{2}}}\end{eqnarray} \tag{ 4.14 }$$

$$\begin{eqnarray}&&{{w}_{{\rm DE}}}=\frac{{{\alpha }_{2}}{{(1-r)}^{2}}-3\left( 1-2r+{{r}^{2}}{{{\rm sin} }^{2}}\theta \right)}{3{{\alpha }_{2}}(1-r)\left[ r({\rm sin} \theta +1)-1 \right]}\;,\end{eqnarray} \tag{ 4.15 }$$

while ${{\Omega }_{{\rm DE}}}$ is just $1-{{\Omega }_{m}}$, i.e.:

$$\begin{eqnarray}&&{{\Omega }_{{\rm DE}}}=\frac{1-r(1+{\rm sin} \theta )}{1-r}.\end{eqnarray} \tag{ 4.16 }$$

Applying the procedure described in appendix B, we conclude that there are three critical point at infinity. These critical points, along with their stability conditions, are presented in table 3. In the same table we include the corresponding values of the observables ${{\Omega }_{{\rm DE}}}$, q, and w_DE, calculated using (4.14), (4.15), and (4.16). These points correspond to Big Rip, sudden, or other forms of singularities [68–73], depending on whether the singularity is reached at finite or infinite time and on their observable features.

Table 3.  The real critical points of the autonomous system (4.6)–(4.7) at infinity, stability conditions, and the corresponding values of the dark energy density parameter ${{\Omega }_{{\rm DE}}}$, of the deceleration parameter q, and of the dark energy equation-of-state parameter w_DE. All points correspond to a form of future, past, or intermediate singularity, depending on the parameters [68–73].

Cr. P.	θ	Stability	${{\Omega }_{{\rm DE}}}$	q	wDE
Q1	0	saddle point	1	$-{\rm sgn} ({{\alpha }_{2}})\infty$	$-{\rm sgn} ({{\alpha }_{2}})\infty$
Q2	${\rm arctan} \left( \frac{{{\alpha }_{1}}}{2} \right)$	unstable for ${{\alpha }_{2}}\gt 0$	$-\infty$	$-{\rm sgn} ({{\alpha }_{2}})\infty$	${\rm sgn} ({{\alpha }_{2}})\infty$
		stable for ${{\alpha }_{2}}\lt 0$
Q3	$\frac{\pi }{2}$	see numerical elaboration	$-\infty$	$\frac{3}{2{{\alpha }_{2}}}$	0

In the previous section we performed the complete phase-space analysis of the physically interesting model with $f(T,{{T}_{G}})=-T+{{\alpha }_{1}}\sqrt{{{T}^{2}}+{{\alpha }_{2}}{{T}_{G}}}$, both at the finite region and at infinity. Thus, in the present section we discuss the corresponding cosmological behavior. As usual, the features of the solutions can be easily deduced by the values of the observables. In particular, $q\lt 0$ ($q\gt 0$) corresponds to acceleration (deceleration), $q=-1$ to the de Sitter solution; ${{w}_{{\rm DE}}}\gt -1$ (${{w}_{{\rm DE}}}\lt -1$) corresponds to quintessence-like (phantom-like) behavior, and ${{\Omega }_{{\rm DE}}}=1$ implies a dark-energy dominated universe.

Point P₁ is stable for the conditions presented in table 1, and thus it can attract the universe at late times. Since the dark energy and matter density parameters are of the same order, this point represents a dark energy–dark matter scaling solution, alleviating the coincidence problem. (Note that in order to handle the coincidence problem, one should provide an explanation of why the present Ω_m and ${{\Omega }_{{\rm DE}}}$ are of the same order, although they follow different evolution behaviors.) However, it has the disadvantage that w_DE is 0, and the universe is not accelerating, as expected [74]. Although this picture is not favored by observations, it may simply imply that the universe today has not yet reached its asymptotic regime.

Point P₂ is stable for the conditions presented in table 1, and therefore it can be the late-times state of the universe. It corresponds to a dark-energy dominated universe that can be accelerating. Interestingly enough, depending on the model parameters, the dark energy equation-of-state parameter can lie in the quintessence regime, it can be equal to the cosmological constant value $-1$, or it can even lie in the phantom regime. These features are a great advantage of the scenario at hand, since they are compatible with observations; moreover, they are obtained only due to the novel features of $f(T,{{T}_{G}})$ gravity, without the explicit inclusion of a cosmological constant or a scalar field, either a canonical or a phantom one.

Point P₃ is stable for the conditions presented in table 1, and therefore it can attract the universe at late times. It has similar features with P₂, but for different parameter regions. Namely, it corresponds to a dark-energy dominated universe that can be accelerating, where the dark energy equation-of-state parameter can lie in the quintessence or phantom regime, or it can be exactly $-1$. These features also make this point a good candidate for the description of Nature.

Point P₄ corresponds to a dark-energy dominated universe that can be accelerating, where the dark energy equation-of-state parameter can lie in the quintessence or phantom regime, or it can be exactly $-1$. However, P₄ is not stable, and thus it cannot attract the universe at late times.

Finally, the present scenario possesses three critical points at infinity, two of which can be stable. They correspond to Big Rip, sudden, or other forms of singularities, depending on the parameter choice. We mention that as the universe moves toward these stable points, the matter density parameter Ω_m will be larger than 1. Although this is not theoretically excluded, growth-index observations could indeed exclude these regions (as it happens in f(R) gravity [10]), and thus the corresponding parameter range that leads the universe to their basin of attraction should be excluded in the model at hand. Such a detailed investigation has not been performed in torsion-based gravity, and therefore it has to be done from the beginning. However, since it lies outside the scope of the present work, it is left for a future project.

In order to present the aforementioned behavior more transparently, we first evolve the autonomous system (4.6)–(4.7) numerically for the parameter choices ${{\alpha }_{1}}=-\sqrt{33}$ and ${{\alpha }_{2}}=4$, assuming the matter to be dust (w_m = 0). The corresponding phase-space behavior is depicted in figure 1. For completeness we also present the projection in the 'Poincaré plane' $(r,\theta )$, where we depict the behavior at both the finite and the infinite region. In this case the universe at late times is attracted by the dark-energy dominated de Sitter attractor P₂, where the effective dark energy behaves like a cosmological constant. At infinity, there is not any stable point, and thus the universe cannot result in any form of singularity.

Zoom In Zoom Out Reset image size

In figure 2 we present the phase-space behavior of the autonomous system (4.6)–(4.7) for the choice ${{\alpha }_{1}}=\frac{1}{2}$ and ${{\alpha }_{2}}=-\frac{1}{2}$ (assuming w_m = 0), and its projection on the 'Poincaré plane' $(r,\theta )$. In this case the attractor at the finite region is the phantom solution P₂. Additionally, the attractor in the infinite region is Q₂, i.e., a future singularity.

Zoom In Zoom Out Reset image size

Finally, in figure 3 some orbits and the corresponding Poincaré projections are present for the choice ${{\alpha }_{1}}=3$ and ${{\alpha }_{2}}=\frac{3}{2}$, with w_m = 0. In this case, the universe is attracted by the quintessence solution P₃. Furthermore, in the infinite region the attractor is Q₃, i.e., a future singularity⁴ .

Zoom In Zoom Out Reset image size

In the present work we studied the dynamical behavior of the recently proposed scenario of $f(T,{{T}_{G}})$ cosmology [55]. This class of modified gravity is based on the quadratic torsion scalar T, which is the Lagrangian of the teleparallel equivalent of General Relativity, and also based on the new quartic torsion scalar T_G, which is the teleparallel equivalent of the Gauss–Bonnet term. Obviously, $f(T,{{T}_{G}})$ theories are more general and cannot be spanned by the simple f(T) ones; additionally, they are different from the $f(R,G)$ class of curvature modified gravity too.

Without loss of generality, as a simple but non-trivial example, capable of revealing the advantages and the new features of the theory, we considered a model where T and T_G corrections are of the same order and thus expected to play an important role at late times. We performed for a spatially flat universe the complete and detailed phase-space behavior, both in the finite and infinite regions, calculating also the values of basic observables such as the various density parameters, the deceleration parameter, and the dark energy equation-of-state parameter.

This scenario exhibits interesting cosmological behaviors. In particular, depending on the model parameters, the universe can result in a dark-energy dominated accelerating solution. Additionally, the dark energy equation-of-state parameter can lie in the quintessence regime, it can be equal to the cosmological constant value $-1$, or it can even lie in the phantom regime. It can also result in a dark energy–dark matter scaling solution, and thus it can alleviate the coincidence problem. Finally, under certain parameter choices, the universe can result in Big Rip, sudden, or other forms of singularities, as is usual in many modified gravitational theories. Definitely, before the scenario at hand can be considered as a good candidate for the description of Nature, a detailed confrontation with observations should be performed. In particular, one should use data from local gravity experiments (Solar System observations), as well as type Ia Supernovae (SNIa), Baryon Acoustic Oscillations (BAO), and Cosmic Microwave Background (CMB) radiation data, in order to impose constraints on the model. These necessary investigations lie beyond the scope of the present work and are left for a future project.

GL was supported by COMISIÓN NACIONAL DE CIENCIAS Y TECNOLOGÍA through Proyecto FONDECYT DE POSTDOCTORADO 2014 grant 3140244 and by DI-PUCV grant 123.730/2013. The research of ENS is implemented within the framework of the Action 'Supporting Postdoctoral Researchers' of the Operational Program 'Education and Lifelong Learning' (Actions Beneficiary: General Secretariat for Research and Technology) and is co-financed by the European Social Fund (ESF) and the Greek State.

For the critical points $({{x}_{c}},{{\Omega }_{mc}})$ of the autonomous system (4.6)–(4.7), presented in table 1, the coefficients of the perturbation equations form a 2 × 2 matrix ${\bf Q}$, which reads:

$$\begin{eqnarray*}\begin{array}{rcl} {{{\bf Q}}_{11}} & = & \frac{{{\alpha }_{1}}(4{{\alpha }_{2}}-3)-9{{\alpha }_{1}}{{x}^{2}}-12x({{\Omega }_{m}}-1)}{2{{\alpha }_{1}}{{\alpha }_{2}}} \\ {{{\bf Q}}_{12}} & = & -\frac{3{{x}^{2}}}{{{\alpha }_{1}}{{\alpha }_{2}}} \\ {{{\bf Q}}_{21}} & = & -\frac{6x{{\Omega }_{m}}}{{{\alpha }_{2}}} \\ {{{\bf Q}}_{22}} & = & -\frac{{{\alpha }_{2}}+3{{x}^{2}}-3}{{{\alpha }_{2}}}. \\ \end{array}\end{eqnarray*}$$

Thus, we can straightforwardly see that using the explicit critical points shown in table 1, the matrix ${\bf Q}$ acquires a simple form that allows for an easy calculation of its eigenvalues. Hence, by examining the sign of the real parts of these eigenvalues, we can classify the corresponding critical point. In particular, if all eigenvalues of a critical point have positive real parts, then this point is unstable; if they all have negative real parts, then it is stable; and if they change sign, then the point is a saddle one. In the following we present the results for each separate point.

Point P₁ has the coordinates:

$$\begin{eqnarray*}&&{{P}_{1}}:(x,{{\Omega }_{m}})=\left( \sqrt{1-\frac{{{\alpha }_{2}}}{3}},\frac{{{\alpha }_{1}}\sqrt{9-3{{\alpha }_{2}}}(6-5{{\alpha }_{2}})+6({{\alpha }_{2}}-3)}{6({{\alpha }_{2}}-3)} \right);\end{eqnarray*}$$

i.e., it exists for either ${{\alpha }_{2}}=\frac{6}{5}$ or $\frac{6}{5}\lt {{\alpha }_{2}}\lt 3,{{\alpha }_{1}}\geqslant -2\sqrt{\frac{3(3-{{\alpha }_{2}})}{{{(5{{\alpha }_{2}}-6)}^{2}}}}$, or ${{\alpha }_{2}}\lt \frac{5}{6},{{\alpha }_{1}}\leqslant 2\sqrt{\frac{3(3-{{\alpha }_{2}})}{{{(5{{\alpha }_{2}}-6)}^{2}}}}$. The eigenvalues of the corresponding linearization matrix are:

$$\begin{eqnarray}&&\begin{array}{l} \left\{ -\frac{\sqrt{{{\alpha }_{1}}\left[ (336-71{{\alpha }_{2}}){{\alpha }_{2}}-288 \right]+32\sqrt{3}{{\left( 3-{{\alpha }_{2}} \right)}^{3/2}}}}{4\sqrt{{{\alpha }_{1}}}{{\alpha }_{2}}}-\frac{3}{4}, \right. \\ \left. \frac{\sqrt{{{\alpha }_{1}}\left[ (336-71{{\alpha }_{2}}){{\alpha }_{2}}-288 \right]+32\sqrt{3}{{\left( 3-{{\alpha }_{2}} \right)}^{3/2}}}}{4\sqrt{{{\alpha }_{1}}}{{\alpha }_{2}}}-\frac{3}{4} \right\}. \\ \end{array}\end{eqnarray} \tag{ A.1 }$$

Therefore, P₁ is a stable spiral for:

$$\begin{eqnarray*}&&{{\alpha }_{2}}\lt 3,\ \ -32\sqrt{3}\sqrt{\frac{{{\left( 3-{{\alpha }_{2}} \right)}^{3}}}{{{\left( 71\alpha _{2}^{2}-336{{\alpha }_{2}}+288 \right)}^{2}}}}\lt {{\alpha }_{1}}\lt 0\end{eqnarray*}$$

or

$$\begin{eqnarray*}&&{{\alpha }_{1}}\lt 0,{{\alpha }_{2}}\leqslant \frac{1}{71}\left( 168-36\sqrt{6} \right)\lesssim 1.12421.\end{eqnarray*}$$

Otherwise it is a saddle. (We have excluded the parameter values that leads to non-hyperbolic critical points.)

Point P₂ has the coordinates

$$\begin{eqnarray*}&&{{P}_{2}}:(x,{{\Omega }_{m}})=\left( \frac{3-\sqrt{3\alpha _{1}^{2}(4{{\alpha }_{2}}-3)+9}}{3{{\alpha }_{1}}},0 \right);\end{eqnarray*}$$

i.e., it exists for either ${{\alpha }_{2}}\lt \frac{3}{4},0\lt {{\alpha }_{1}}\leqslant \sqrt{\frac{3}{3-4{{\alpha }_{2}}}}$, or ${{\alpha }_{1}}\ne 0,{{\alpha }_{2}}=\frac{3}{4}$, or ${{\alpha }_{2}}\gt \frac{3}{4},{{\alpha }_{1}}\lt 0$. The eigenvalues of the linearization matrix are:

$$\begin{eqnarray}&&\begin{array}{l} \left\{ \frac{2\sqrt{3\alpha _{1}^{2}(4{{\alpha }_{2}}-3)+9}}{\alpha _{1}^{2}{{\alpha }_{2}}}-\frac{6}{\alpha _{1}^{2}{{\alpha }_{2}}}+\frac{6}{{{\alpha }_{2}}}-5, \right. \\ \qquad \left. \frac{\sqrt{3\alpha _{1}^{2}(4{{\alpha }_{2}}-3)+9}}{\alpha _{1}^{2}{{\alpha }_{2}}}-\frac{3}{\alpha _{1}^{2}{{\alpha }_{2}}}+\frac{3}{{{\alpha }_{2}}}-4 \right\}. \\ \end{array}\end{eqnarray} \tag{ A.2 }$$

Hence, it is a stable node for either:

$$\begin{eqnarray*}&&{{\alpha }_{2}}\lt 0,\ \ 0\lt {{\alpha }_{1}}\lt 2\sqrt{\frac{3(3-{{\alpha }_{2}})}{{{\left( 5{{\alpha }_{2}}-6 \right)}^{2}}}},\end{eqnarray*}$$

or

$$\begin{eqnarray*}&&\frac{6}{5}\lt {{\alpha }_{2}}\leqslant 3,\ \ {{\alpha }_{1}}\lt -2\sqrt{\frac{3(3-{{\alpha }_{2}})}{{{\left( 5{{\alpha }_{2}}-6 \right)}^{2}}}},\end{eqnarray*}$$

or

$$\begin{eqnarray*}&&{{\alpha }_{2}}\gt 3,\ \ {{\alpha }_{1}}\lt 0.\end{eqnarray*}$$

Additionally, it is an unstable node for:

$$\begin{eqnarray*}&&0\lt {{\alpha }_{2}}\lt \frac{3}{4},\ \ 0\lt {{\alpha }_{1}}\lt \sqrt{\frac{3}{3-4{{\alpha }_{2}}}}.\end{eqnarray*}$$

Finally, for the remaining parameter range in the hyperbolic domain, the point behaves as a saddle.

Point P₃ has the coordinates:

$$\begin{eqnarray*}&&{{P}_{3}}:(x,{{\Omega }_{m}})=\left( \frac{3+\sqrt{3\alpha _{1}^{2}(4{{\alpha }_{2}}-3)+9}}{3{{\alpha }_{1}}},0 \right);\end{eqnarray*}$$

i.e., it exists for ${{\alpha }_{2}}\lt \frac{3}{4},\;0\lt {{\alpha }_{1}}\leqslant \sqrt{\frac{3}{3-4{{\alpha }_{2}}}}$, or ${{\alpha }_{2}}\geqslant \frac{3}{4},{{\alpha }_{1}}\gt 0$. The eigenvalues of the corresponding linearization matrix are:

$$\begin{eqnarray}&&\begin{array}{l} \left\{ -\frac{2\sqrt{3}\sqrt{4\alpha _{1}^{2}{{\alpha }_{2}}-3\alpha _{1}^{2}+3}}{\alpha _{1}^{2}{{\alpha }_{2}}}-\frac{6}{\alpha _{1}^{2}{{\alpha }_{2}}}+\frac{6}{{{\alpha }_{2}}}-5, \right. \\ \qquad \left. -\frac{\sqrt{3}\sqrt{4\alpha _{1}^{2}{{\alpha }_{2}}-3\alpha _{1}^{2}+3}}{\alpha _{1}^{2}{{\alpha }_{2}}}-\frac{3}{\alpha _{1}^{2}{{\alpha }_{2}}}+\frac{3}{{{\alpha }_{2}}}-4 \right\}. \\ \end{array}\end{eqnarray} \tag{ A.3 }$$

That is, it is a stable node for:

$$\begin{eqnarray*}&&{{\alpha }_{1}}\gt 0,\ \ {{\alpha }_{2}}\geqslant \frac{6}{5},\end{eqnarray*}$$

and it is an unstable node for

$$\begin{eqnarray*}&&{{\alpha }_{2}}\lt 0,\ \ 0\lt {{\alpha }_{1}}\lt \frac{\sqrt{3}}{\sqrt{3-4{{\alpha }_{2}}}};\end{eqnarray*}$$

otherwise it is a saddle with the exclusion of the parameter values that leads to a non-hyperbolic critical point.

Point P₄ has the coordinates ${{P}_{4}}:(x,{{\Omega }_{m}})=(0,0)$, and it exists always. The eigenvalues of the linearization matrix read:

$$\begin{eqnarray*}&&\left\{ 2-\frac{3}{2{{\alpha }_{2}}},\;\frac{3}{{{\alpha }_{2}}}-1 \right\}.\end{eqnarray*}$$

Therefore, it is an unstable node for $\frac{3}{4}\lt {{\alpha }_{2}}\lt 3$, and it is non-hyperbolic for ${{\alpha }_{2}}\in \left\{ \frac{3}{4},3 \right\}$; otherwise it is a saddle.

The above results are summarized in table A.1 .

Table A.1.  The real and physically interesting critical points at the finite region of the autonomous system (4.6)–(4.7), their existence conditions, and the corresponding eigenvalues ${{\nu }_{1}},{{\nu }_{2}}$ of the matrix ${\bf Q}$ of the perturbation equations. We denote ${{\Omega }_{m1}}=\frac{{{\alpha }_{1}}\sqrt{9-3{{\alpha }_{2}}}(6-5{{\alpha }_{2}})+6({{\alpha }_{2}}-3)}{6({{\alpha }_{2}}-3)}$.

Cr. P.	x	Ωm	Existence	ν1	ν2
P1	$\sqrt{1-\frac{{{\alpha }_{2}}}{3}}$	Ωm1	${{\alpha }_{2}}=\frac{6}{5}$ or	$-\frac{\sqrt{{{\alpha }_{1}}\left[ (336-71{{\alpha }_{2}}){{\alpha }_{2}}-288 \right]+32\sqrt{3}{{\left( 3-{{\alpha }_{2}} \right)}^{3/2}}}}{4\sqrt{{{\alpha }_{1}}}{{\alpha }_{2}}}-\frac{3}{4}$	$\frac{\sqrt{{{\alpha }_{1}}\left[ (336-71{{\alpha }_{2}}){{\alpha }_{2}}-288 \right]+32\sqrt{3}{{\left( 3-{{\alpha }_{2}} \right)}^{3/2}}}}{4\sqrt{{{\alpha }_{1}}}{{\alpha }_{2}}}-\frac{3}{4}$
			$\frac{6}{5}\lt {{\alpha }_{2}}\lt 3,-2\sqrt{\frac{3(3-{{\alpha }_{2}})}{{{\left( -6+5{{\alpha }_{2}} \right)}^{2}}}}\leqslant {{\alpha }_{1}}\leqslant 0$
			or ${{\alpha }_{2}}=\frac{6}{5}$
			or ${{\alpha }_{2}}\lt \frac{6}{5},0\leqslant {{\alpha }_{1}}\leqslant 2\sqrt{\frac{3(3-{{\alpha }_{2}})}{{{\left( -6+5{{\alpha }_{2}} \right)}^{2}}}}$
P2	$\frac{3-\sqrt{3\alpha _{1}^{2}(4{{\alpha }_{2}}-3)+9}}{3{{\alpha }_{1}}}$	0	${{\alpha }_{2}}\lt \frac{3}{4},0\lt {{\alpha }_{1}}\leqslant \sqrt{\frac{3}{3-4{{\alpha }_{2}}}}$ or	$\frac{2\sqrt{3\alpha _{1}^{2}(4{{\alpha }_{2}}-3)+9}}{\alpha _{1}^{2}{{\alpha }_{2}}}-\frac{6}{\alpha _{1}^{2}{{\alpha }_{2}}}+\frac{6}{{{\alpha }_{2}}}-5$	$\frac{\sqrt{3\alpha _{1}^{2}(4{{\alpha }_{2}}-3)+9}}{\alpha _{1}^{2}{{\alpha }_{2}}}-\frac{3}{\alpha _{1}^{2}{{\alpha }_{2}}}+\frac{3}{{{\alpha }_{2}}}-4$
			${{\alpha }_{1}}\ne 0,{{\alpha }_{2}}=\frac{3}{4}$ or
			${{\alpha }_{2}}\gt \frac{3}{4},{{\alpha }_{1}}\lt 0$
P3	$\frac{3+\sqrt{3\alpha _{1}^{2}(4{{\alpha }_{2}}-3)+9}}{3{{\alpha }_{1}}}$	0	${{\alpha }_{2}}\lt \frac{3}{4},0\lt {{\alpha }_{1}}\leqslant \sqrt{\frac{3}{3-4{{\alpha }_{2}}}}$ or	$-\frac{2\sqrt{3}\sqrt{4\alpha _{1}^{2}{{\alpha }_{2}}-3\alpha _{1}^{2}+3}}{\alpha _{1}^{2}{{\alpha }_{2}}}-\frac{6}{\alpha _{1}^{2}{{\alpha }_{2}}}+\frac{6}{{{\alpha }_{2}}}-5$	$-\frac{\sqrt{3}\sqrt{4\alpha _{1}^{2}{{\alpha }_{2}}-3\alpha _{1}^{2}+3}}{\alpha _{1}^{2}{{\alpha }_{2}}}-\frac{3}{\alpha _{1}^{2}{{\alpha }_{2}}}+\frac{3}{{{\alpha }_{2}}}-4$
			${{\alpha }_{2}}\geqslant \frac{3}{4},{{\alpha }_{1}}\gt 0$
P4	0	0	Always	$2-\frac{3}{2{{\alpha }_{2}}}$	$\frac{3}{{{\alpha }_{2}}}-1$

We introduce the new coordinates $(r,\theta )$, defined by:

$$\begin{eqnarray}\begin{array}{rcl} x & = & \frac{r}{1-r}{\rm cos} \theta \\ {{\Omega }_{m}} & = & \frac{r}{1-r}{\rm sin} \theta , \\ \end{array}\end{eqnarray} \tag{ B.1 }$$

with $\theta \in \left[ 0,\frac{\pi }{2} \right]$ and $r\in \left[ 0,1 \right)$. The limit $r\to {{1}^{-}}$ corresponds to ${{R}^{2}}\equiv {{x}^{2}}+\Omega _{m}^{2}\to \infty .$ Note that the physical region of the plane $(r,\theta )$ that corresponds to $0\leqslant x,\;0\leqslant {{\Omega }_{m}}\leqslant 1$, is given by:

$$\begin{eqnarray}&&\begin{array}{l} \left\{ (r,\theta ):0\leqslant r\leqslant \frac{1}{2},0\leqslant \theta \leqslant \frac{\pi }{2} \right\}\cup \\ \qquad \left\{ (r,\theta ):\frac{1}{2}\lt r\lt 1,0\leqslant \theta \leqslant {\rm arcsin} \left( \frac{1-r}{r} \right) \right\}. \\ \end{array}\end{eqnarray} \tag{ B.2 }$$

The leading terms of the equations for $r^{\prime}$ and $\theta ^{\prime}$ as $r\to {{1}^{-}}$ are:

$$\begin{eqnarray}&&r^{\prime} \to \frac{3{{{\rm cos} }^{2}}(\theta )[{{\alpha }_{1}}({\rm cos} (2\theta )-3)-2{\rm sin} (2\theta )]}{4{{\alpha }_{1}}{{\alpha }_{2}}(1-r)}\end{eqnarray} \tag{ B.3 }$$

$$\begin{eqnarray}&&\theta ^{\prime} \to -\frac{3{\rm sin} (\theta ){{{\rm cos} }^{2}}(\theta )\left[ {{\alpha }_{1}}{\rm cos} (\theta )-2{\rm sin} (\theta ) \right]}{2{{\alpha }_{1}}{{\alpha }_{2}}{{\left( 1-r \right)}^{2}}}.\end{eqnarray} \tag{ B.4 }$$

Hence, the fixed points at infinity (i.e., for $r\to {{1}^{-}}$) are obtained by setting $\theta ^{\prime} =0$ and solving for θ.

Let us denote a generic fixed point by $\theta ={{\theta }^{*}}$. The stability of this point is studied by analyzing first the stability of the angular coordinates from equation (B.4) and then deducing, from the sign of the equation (B.3), the stability on the radial direction⁵ . A fixed point $\theta ={{\theta }^{*}}$ is said to be stable if both:

$$\begin{eqnarray}&&\frac{{\rm d}\theta ^{\prime} }{{\rm d}\theta }\left| _{\theta ={{\theta }^{*}}} \right.\lt 0,r^{\prime} \left| _{\theta ={{\theta }^{*}}} \right.\gt 0.\end{eqnarray} \tag{ B.5 }$$

The first condition implies stability of the angular coordinate θ. The second condition implies that the r-values increase before reaching the limit value r = 1 (i.e., before the boundary 'at infinity' is reached) at the fixed point $\theta ={{\theta }^{*}}$. Similarly, a fixed point $\theta ={{\theta }^{*}}$ is said to be unstable if both:

$$\begin{eqnarray}&&\frac{{\rm d}\theta ^{\prime} }{{\rm d}\theta }\left| _{\theta ={{\theta }^{*}}} \right.\gt 0,r^{\prime} \left| _{\theta ={{\theta }^{*}}} \right.\lt 0.\end{eqnarray} \tag{ B.6 }$$

Finally, it is a saddle point if either:

$$\begin{eqnarray}&&\frac{{\rm d}\theta ^{\prime} }{{\rm d}\theta }\left| _{\theta ={{\theta }^{*}}} \right.\gt 0,r^{\prime} \left| _{\theta ={{\theta }^{*}}} \right.\gt 0\end{eqnarray} \tag{ B.7 }$$

or

$$\begin{eqnarray}&&\frac{{\rm d}\theta ^{\prime} }{{\rm d}\theta }\left| _{\theta ={{\theta }^{*}}} \right.\lt 0,r^{\prime} \left| _{\theta ={{\theta }^{*}}} \right.\lt 0.\end{eqnarray} \tag{ B.8 }$$

In summary, the fixed points of the autonomous system at hand at infinity are the following:

• ${{Q}_{1}}:{{\theta }^{*}}=0,{{r}^{*}}=1,x=\infty ,{{\Omega }_{m}}=0$.

Since from the definition (4.4) we have $x=\sqrt{1+\frac{2{{\alpha }_{2}}}{3}(1+\frac{{\dot{H}}}{{{H}^{2}}})}$, we deduce that the corresponding cosmological solution satisfies ${\rm sign}({{\alpha }_{2}}\dot{H})\frac{|\dot{H}|}{{{H}^{2}}}\to \infty$ or $H\to 0$ ($\dot{H}$ is bounded, with ${\rm sign}({{\alpha }_{2}}\dot{H})\gt 0$). Since $q={{w}_{{\rm DE}}}=-{\rm sgn} ({{\alpha }_{2}})\infty$, the point represents a super-accelerated phantom solution for ${{\alpha }_{2}}\gt 0$, where eventually the universe ends in Big Rip, sudden, or other forms of singularities (depending on whether the singularity is reached at finite or infinite time, what its features are, etc) [68–73]. For ${{\alpha }_{2}}\lt 0$, it is a decelerating solution where the universe asymptotically stops expanding. Since $\left( \frac{{\rm d}\theta ^{\prime} }{{\rm d}\theta },r^{\prime} \right){{|}_{\theta =0}}=\left( -\frac{3}{2{{\alpha }_{2}}},-\frac{3}{2{{\alpha }_{2}}} \right)$, we conclude that Q₁ is always a saddle point.

• ${{Q}_{2}}:{{\theta }^{*}}={\rm arctan} (\frac{{{\alpha }_{1}}}{2}),{{\alpha }_{1}}\ne 0,{{r}^{*}}=1,\frac{{{\Omega }_{m}}}{x}\to \frac{{{\alpha }_{1}}}{2},x\to \infty ,{{\Omega }_{m}}\to \infty$. Since $0\leqslant \theta \leqslant \frac{\pi }{2},$ then ${{\alpha }_{1}}\gt 0$. Since $\left( \frac{{\rm d}\theta ^{\prime} }{{\rm d}\theta },r^{\prime} \right){{|}_{\theta =\frac{\pi }{2}}}=\left( \frac{6}{(1+\alpha _{1}^{2}){{\alpha }_{2}}},-\frac{12}{(1+\alpha _{1}^{2}){{\alpha }_{2}}} \right)$, we deduce that Q₂ is unstable for ${{\alpha }_{2}}\gt 0$ or stable for ${{\alpha }_{2}}\lt 0$, as confirmed in figure 2. Since at this point ${{\Omega }_{m}}$ diverges, it corresponds to some form of future (respectively past) singularity for ${{\alpha }_{2}}\lt 0$ (respectively ${{\alpha }_{2}}\gt 0$) [68–73]. Its detailed classification for the various parameter regions lies beyond the scope of the present work.
• ${{Q}_{3}}:{{\theta }^{*}}=\frac{\pi }{2},{{r}^{*}}=1,x=0,{{\Omega }_{m}}=\infty$. Since $\left( \frac{{\rm d}\theta ^{\prime} }{{\rm d}\theta },r^{\prime} \right){{|}_{\theta =\frac{\pi }{2}}}=\left( 0,0 \right)$, we cannot rely on the linearization to examine the stability, and therefore we need to resort to numerical examination (see figures 1, 2). Since at this point ${{\Omega }_{m}}$ diverges, it corresponds to some form of future, past, or intermediate singularity [68–73].

The above results are summarized in table B.1 .

Table B.1.  The real critical points of the autonomous system (4.6)–(4.7) at infinity, their existence conditions, the corresponding values of $\frac{{\rm d}\theta ^{\prime} }{{\rm d}\theta }$ and $r^{\prime}$, and the resulting stability conditions.

Cr. P.	${{\theta }^{*}}$	${{\left. \frac{{\rm d}\theta ^{\prime} }{{\rm d}\theta } \right|}_{\theta ={{\theta }^{*}}}}$	${{\left. r^{\prime} \right|}_{\theta ={{\theta }^{*}}}}$	Stability
Q1	0	$-\frac{3}{2{{\alpha }_{2}}}$	$-\frac{3}{2{{\alpha }_{2}}}$	saddle point
Q2	${\rm arctan} \left( \frac{{{\alpha }_{1}}}{2} \right)$	$\frac{6}{(1+\alpha _{1}^{2}){{\alpha }_{2}}}$	$-\frac{12}{(1+\alpha _{1}^{2}){{\alpha }_{2}}}$	unstable for ${{\alpha }_{2}}\gt 0$
				stable for ${{\alpha }_{2}}\lt 0$
Q3	$\frac{\pi }{2}$	0	0	see numerical elaboration