Dynamical system analysis of generalized energy-momentum-squared gravity

TL;DR

This paper analyzes the dynamics of a modified gravity theory, energy-momentum squared gravity, using dynamical systems to explore its phase space and cosmological implications, including accelerated expansion and matter eras.

Contribution

It introduces a dynamical system analysis for the $f(R, T^2)$ gravity model, revealing diverse cosmological behaviors and phase space structures depending on model parameters.

Findings

Phase space complexity varies with parameters $m$ and $n$.

Models can reproduce matter, dark energy, and accelerated expansion epochs.

Richer phase space allows for multiple cosmological scenarios.

Abstract

In this work we have investigated the dynamics of a recent modification to the general theory of relativity, the energy-momentum squared gravity model $f(R,\mathbf{T^2})$, where $R$ represents the scalar curvature and $\mathbf{T^2}$ the square of the energy-momentum tensor. By using dynamical system analysis for various types of gravity functions $f(R,\mathbf{T^2})$, we have studied the structure of the phase space and the physical implications of the energy-momentum squared coupling. In the first case of functional where $f(R,\mathbf{T^2})=f_0 R^n(\mathbf{T^2})^m$, with $f_0$ constant, we have shown that the phase space structure has a reduced complexity, with a high sensitivity to the values of the $m$ and $n$ parameters. Depending on the values of the $m$ and $n$ parameters, the model exhibits various cosmological epochs, corresponding to matter eras, solutions associated with an accelerated expansion, or decelerated periods. The second model studied corresponds to the $f(R,\mathbf{T^2})=\alpha R^n+\beta\mathbf{(T^2)}^m$ form with $\alpha, \beta$ constant parameters. In this case, a richer phase space structure is obtained which can recover different cosmological scenarios, associated to matter eras, de--Sitter solutions, and dark energy epochs. Hence, this model represents an interesting cosmological model which can explain the current evolution of the Universe and the emergence of the accelerated expansion as a geometrical consequence.

Tables

Table 1. Table 1: The structure of the phase space for the f ​ ( R , 𝐓 𝟐 ) = f 0 ​ R n ​ ( 𝐓 𝟐 ) m 𝑓 𝑅 superscript 𝐓 2 subscript 𝑓 0 superscript 𝑅 𝑛 superscript superscript 𝐓 2 𝑚 f(R,\mathbf{T^{2}})=f_{0}R^{n}(\mathbf{T^{2}})^{m} model. We have defined the following quantity: τ = 4 ​ m 2 ​ ( n ​ ( 49 ​ n − 78 ) + 9 ) + 4 ​ m ​ ( n ​ ( 7 ​ n ​ ( 8 ​ n − 19 ) + 78 ) − 9 ) + ( n ​ ( 8 ​ n − 13 ) + 3 ) 2 𝜏 4 superscript 𝑚 2 𝑛 49 𝑛 78 9 4 𝑚 𝑛 7 𝑛 8 𝑛 19 78 9 superscript 𝑛 8 𝑛 13 3 2 \tau=\sqrt{4m^{2}(n(49n-78)+9)+4m(n(7n(8n-19)+78)-9)+(n(8n-13)+3)^{2}} .

Cr.P.	$x_2$	$x_4$	$x_5$
$P_1$	$\frac{m\left(n+5\right)+4n-5}{\left(m+2\right)n^2+\left(2m-3\right)n-m+1}$	0	$\frac{-3mn\left(n+2\right)+3m+n\left(11-7n\right)-3}{3m\left(n\left(n+2\right)-1\right)+6n^2-9n+3}$
$P_2$	$\frac{52m^{2}n-12m^2+56m{n}^2-52mn+12m+16n^3-28n^2+13n-3+\left(1-2m\right)\tau }{4n^2\left(2m+n-1\right)\left(2m+2n-1\right)}$	$\frac{-3+m\left(6-26n\right)+\left(13-8n\right)n+\tau }{4n^2}$	$-\frac{3+6m\left(4m-3\right)-13n+26mn+8n^2+\tau }{6\left(2m+n-1\right)\left(2m+2n-1\right)}$
$P_3$	$\frac{52m^{2}n-12m^2+56m{n}^2-52mn+12m+16n^3-28n^2+13n-3-\left(1-2m\right)\tau }{4n^2\left(2m+n-1\right)\left(2m+2n-1\right)}$	$\frac{-3+m\left(6-26n\right)+\left(13-8n\right)n-\tau }{4n^2}$	$\frac{-26mn+6m\left(3-4m\right)-8n^2+13n-3+\tau }{6\left(2m+n-1\right)\left(2m+2n-1\right)}$

Table 2. Table 2: Critical points, existence condition, effective state parameter and acceleration for the f ​ ( R , 𝐓 𝟐 ) = α ​ R n + β ​ ( 𝐓 𝟐 ) m 𝑓 𝑅 superscript 𝐓 2 𝛼 superscript 𝑅 𝑛 𝛽 superscript superscript 𝐓 2 𝑚 f(R,\mathbf{T^{2}})=\alpha R^{n}+\beta(\mathbf{T^{2}})^{m} model considering a dust matter w = 0 𝑤 0 w=0 .

Cr.P.	$x_3$	$x_4$	$x_5$	Existence	$w_{\mathrm{eff}}$	Acceleration
$P_1$	0	0	0	Always	$1/3$	Never
$P_2$	0	0	$-1$	Always	$1/3$	Never
$P_3$	$\frac{n\left(4n-5\right)}{2n^2-3n+1}$	0	0	Always	$\frac{-6n^2+7n+1}{6n^2-9n+3}$
$P_4$	$\frac{n\left(4n-5\right)}{2n^2-3n+1}$	0	$-1$	Always	$\frac{-6n^2+7n+1}{6n^2-9n+3}$
$P_5$	$\frac{4n-3}{2n}$	$\frac{-8n^2+13n-3}{2n^2}$	0	$\frac{13-\sqrt{73}}{16}\le n\le \frac{3}{10}\vee 1\le n\le \frac{13+\sqrt{73}}{16}$	$\frac{1}{n}-1$	Never
$P_6$	2	$-\frac{2m\left(n-2\right)}{n}$	$-1$	$0\le -\frac{2m\left(n-2\right)}{n}\le 1$	$-1$	Always
$P_7$	0	$10m$	$-1$	$0\le m\le \frac{1}{10}$	$1/3$	Never

Table 3. Table 3: Stability of the critical points for the f ​ ( R , 𝐓 𝟐 ) = α ​ R n + β ​ ( 𝐓 𝟐 ) m 𝑓 𝑅 superscript 𝐓 2 𝛼 superscript 𝑅 𝑛 𝛽 superscript superscript 𝐓 2 𝑚 f(R,\mathbf{T^{2}})=\alpha R^{n}+\beta(\mathbf{T^{2}})^{m} model with a pressureless matter ( w = 0 𝑤 0 w=0 ). The (approximated) real numbers appearing in the inequalities above come from the numerical analysis of the sign of the eigenvalues of the fixed point

Cr.P.	Stable	Unstable	Saddle
$P_1$	Never
$P_2$	Never
$P_3$
	$\vee$		$\vee$
	$\vee$ $n>\frac{1}{16}\left(\sqrt{73}+13\right)\wedge m>\frac{1}{2}$		$\vee$
$P_4$
$P_5$		Never
	$\vee m>1/2\wedge \left(1.327⪅n⪅1.346\right)$
$P_6$		Never
$P_7$		Never

Equations

$$S=\frac{1}{2{\kappa }^2}\int d^{4}x\sqrt{-g}f(R,{\mathbf{T}}^2)+S_{\mathrm{m}},$$

$$R_{\mu \nu }{f}_R+g_{\mu \nu }\square f_R-{\nabla }_{\mu }{\nabla }_{\nu }{f}_R-\frac{1}{2}{g}_{\mu \nu }f={\kappa }^2{T}_{\mu \nu }-f_{{\mathbf{T}}^2}{\Theta }_{\mu \nu },$$

$$\Theta_{\mu\nu}=\frac{\delta(\mathbf{T^{2}})}{\delta g^{\mu\nu}}=\frac{\delta(T^{\alpha\beta}T_{\alpha\beta})}{\delta g^{\mu\nu}}=-2L_{\rm m}\Big{(}T_{\mu\nu}-\frac{1}{2}g_{\mu\nu}T\Big{)}-T\,T_{\mu\nu}+2T^{\alpha}_{\mu}T_{\nu\alpha}-4T^{\alpha\beta}\frac{\partial^{2}L_{\rm m}}{\partial g^{\mu\nu}\partial g^{\alpha\beta}}\,,$$

$${\kappa }^2{\nabla }^{\mu }{T}_{\mu \nu }=-\frac{1}{2}{g}_{\mu \nu }{\nabla }^{\mu }f+{\nabla }^{\mu }(f_{{\mathbf{T}}^2}{\Theta }_{\mu \nu }).$$

$$d{s}^2=-d{t}^2+a^2(t){\delta }_{ik}d{x}^{i}d{x}^k,$$

$$\Theta_{\mu\nu}=-\Big{(}\rho^{2}+4p\rho+3p^{2}\Big{)}u_{\mu}u_{\nu}\,.$$

$$\displaystyle-3f_{R}\Big{(}\dot{H}+H^{2}\Big{)}+\frac{f}{2}+3H\dot{f_{R}}$$

$$-f_R(\dot{H}+3H^2)+\frac{1}{2}f+\ddot{f_R}+2H\dot{f_R}$$

$${\mathbf{\Theta }}^2:={\Theta }_{\mu \nu }{\Theta }^{\mu \nu }={\rho }^2+4p\rho +3p^2$$

$${\kappa }^2(\dot{\rho }+3H(\rho +p))$$

$$p=w\rho ,$$

$${\mathbf{\Theta }}^2=(1+4w+3w^2){\rho }^2,$$

$$\dot{\rho }+3H(w+1)\rho$$

$$x_1=\frac{{\dot{f}}_R}{f_{R}H},x_2=\frac{f}{6H^2{f}_R},x_3=\frac{R}{6H^2},x_4=\frac{\rho }{3H^2{f}_R},x_5=\left(3w^2+4w+1\right)\rho f_{{\mathbf{T}}^2}.$$

$$x_3+x_4+x_4{x}_5-x_1-x_2=1.$$

$$\frac{d{x}_1}{dN}$$

$$\frac{d{x}_2}{dN}$$

$$\frac{d{x}_3}{dN}$$

$$\frac{d{x}_4}{dN}$$

$$\frac{d{x}_5}{dN}$$

$$\mho =\frac{\dot{R}}{6H^3},\Xi =\frac{\dot{f}}{6f_R{H}^3},\Psi =\frac{\dot{H}}{H^2},\Delta =f_{R}H{\dot{f}}_{{\mathbf{T}}^2},\Pi =\frac{\dot{\rho }}{3f_R{H}^3}.$$

$$\displaystyle\Delta=-\frac{\Pi\Big{(}(w(3w+8)+1)x_{5}+(w+1)(3w+1)\Big{)}+3(w+1)(3w+1)x_{4}\left(x_{5}+w+1\right)}{3(w+1)^{2}(3w+1)^{2}x_{4}{}^{2}}\,.$$

$$\Psi =x_3-2,\Xi =\frac{\left(3w^2+1\right)\Pi x_5}{3w^2+4w+1}+\mho ,$$

$$\frac{d{x}_2}{dN}$$

$$\frac{d{x}_3}{dN}$$

$$\frac{d{x}_4}{dN}$$

$$\frac{d{x}_5}{dN}$$

$$x_3=\Gamma +3w{x}_4+2x_1+3x_2-1.$$

$$w_{\mathrm{eff}}=-1-\frac{2}{3}\left(x_3-2\right).$$

$$\Delta=\frac{2mnx_{2}^{2}}{3\left(3w^{2}+1\right)x_{3}^{2}x_{4}^{3}}\Big{(}2(m-1)\Pi x_{3}+nx_{4}\mho\Big{)}\,,$$

$$x_3=n{x}_2.$$

$$x_2(n-1)+x_4+x_4{x}_5-1=\frac{(n-1)\mho }{n{x}_2}+\frac{2m\Pi }{x_4},$$

$$\Delta$$

$$\Pi$$

$$\Xi$$

Full text

I Introduction

The discovery of the accelerated expansion of the Universe Perlmutter et al. (1999); Spergel et al. (2003) at the turn of the last century has come as the most unexpected and surprising result for the scientific society. The reason behind the reaction is quite obvious as gravity being an attractive force will compel the Universe and all the matter present inside it to contract. So the expansion of the Universe would gradually slow down and finally reach a situation where it stops totally. After this gravity will pull it back and the Universe would undergo a contraction. But the observational evidences are saying a totally different story. This is where our standard knowledge of physics is falling short and we are compelled to search for some new physics which will help us to explain the phenomenon that our observations are showing.

Over the last two decades, the scientific society have left no stones un-turned in its quest for a suitable physical theory that will explain the accelerated expansion Copeland et al. (2006) of the Universe. The whole effort can be broadly classified into two categories. The first one targets the nature of matter that fills the Universe. This theory tells us that the Universe is filled with a mysterious negative pressure component termed as “dark energy” which provides an antigravitating stress to sustain not only an expanding Universe but also fuels it to reach a level of accelerated expansion. The most common way of doing this is by introducing a cosmological constant $\Lambda$ in the Einstein’s equations of general relativity. The idea seems to be consistent but the concept is still in its infancy and is plagued by a lot of shortcomings. The most prominent one being its invisible nature and hence the term “mysterious”. Moreover there is the cosmological constant problem to deal with Martin (2012) which is associated to the lack of consistency between the observed values of vacuum energy density (given by a very small value of the cosmological constant) and the theoretically predicted large value given by the quantum field theory.

The second approach aims at modifying the geometry of spacetime, i.e. the Einstein’s gravity in the general theory of relativity (GR) at large distances, specifically beyond our Solar System to produce accelerating cosmological solutions Carroll et al. (2004); Cognola et al. (2008); Clifton et al. (2012). This has given rise to the concept of ”modified gravity”, with numerous theories available in the literature. Extensive reviews in modified gravity theories can be found in the Refs. Nojiri et al. (2017); Nojiri and Odintsov (2006a); Capozziello et al. (2019). Many of the theories of modified gravity aims at modifying the linear function of scalar curvature, $R$ responsible for the Einstein tensor in the Einstein equations of GR. So it is obvious that the alterations are brought about in such a way so as to generalize the gravitational Lagrangian which takes a special form $\mathcal{L}_{\rm GR}=R$ in case of GR. An extensively studied theory in this context is the $f(R)$ gravity where the gravitational Lagrangian $\mathcal{L}_{\rm GR}=R$ is replaced by an analytic function of $R$, i.e., $\mathcal{L}_{f(R)}=f(R)$. Via this generalization, we can explore the non-linear effects of the scalar curvature $R$ in the evolution of the universe by choosing a suitable function for $f(R)$. In this specific case, extensive reviews on this theory are available in the Refs. De Felice and Tsujikawa (2010); Sotiriou and Faraoni (2010).

From a theoretical point of view, the viability of $f(R)$ dark energy models have been studied in Amendola et al. (2007a). In this paper the authors ruled out the $f(R)$ theories where a power of $R$ is dominant at large or small $R$. The effects of a non-minimal curvature-matter coupling was studied in Azizi and Yaraie (2014), and constraints of the $f(R)$ dark energy models were derived in Amendola et al. (2007b). In Ref. Sotiriou (2006), the author studied the interplay between scalar–tensor theory and $f(R)$ theories of gravity considering the Palatini formalism. A specific scheme for $f(R)$ reconstruction was developed recently Nojiri and Odintsov (2006b), while large scale structure of $f(R)$ gravity was investigated in Ref. Song et al. (2007). Moreover, various papers have investigated different aspects of the latter theory by considering various techniques Shah and Samanta (2019); Odintsov and Oikonomou (2017); Bamba et al. (2014), and a survey of the generic f(R) models in various formulations is carried out in Capozziello et al. (2010). Further generalizations to the $f(R)$ modified theories of gravity have been proposed by introducing some couplings between the geometrical quantities and the matter sector. One interesting model is the one where the Lagrangian is constructed by considering a generic function of the Ricci scalar $R$ and of the trace of the stress-energy tensor $T$. Such modifications gave rise to $f(R,T)$ theories Harko et al. (2011a); Katırcı and Kavuk (2014); Nagpal et al. (2019). Moreover since scalar fields play a fundamental role in cosmology, $f(R,T^{\phi})$ theories were proposed by Harko et al in Ref. Harko et al. (2011a), where $T^{\phi}$ is the trace of the stress energy of the scalar field. A different type of coupling between geometry and matter was proposed Haghani et al. (2013), the generic $f(R,T,R_{\mu\nu}T^{\mu\nu})$ gravity theory. This is a more generic gravity theory in which matter is nonminimally coupled to geometry, the Lagrangian corresponding to the gravitational field has a general dependence of the Ricci scalar, the trace of the matter energy-momentum tensor, and the contraction between the Ricci tensor and the matter energy-momentum tensor. Further, in Bertolami et al. (2007); Harko (2010), a model with a non-minimally coupling between the density Lagrangian matter and the curvature $R$ was introduced. In this model, the Lagrangian corresponds to $f_{1}(R)+(1+\lambda f_{2}(R))L_{\rm m}$ where $\lambda$ is a constant and $L_{\rm m}$ is the matter Lagrangian density. Within these models in which the matter field is non-minimally coupled to gravity any particle is subject to an extra force appearing in a direction which is orthogonal to the four-velocity Bertolami et al. (2007). The latter model proposed in Bertolami et al. (2007) was extended to the case of generic couplings to both matter and gravity in a recent paper Harko (2008), considering a specific Lagrangian given by $f_{1}(R)+G(L_{\rm m})f_{2}(R)$. Furthermore, for the non-minimal geometry coupling models the Palatini formulation has been proposed in Harko et al. (2011b). In this context, a further specific extension related to the latter gravity theories was proposed in Harko and Lobo (2010) by embedding into the Lagrangian an arbitrary function of the $f(R,L_{\rm m})$ type.

More comprehensive ideas and reviews on modified gravity theories from different points of view has been considered in Nojiri et al. (2017); De Felice and Tsujikawa (2010); D’Agostino and Luongo (2018). The advances in the recent cosmology using a dynamical system approach in dark energy and modified gravity theories have appeared recently Bahamonde et al. (2018). In continuation of the above generalization procedure for the $f(R,L_{m})$ theory we can also choose to modify the corresponding Lagrangian by including some analytic function of $T_{\mu\nu}T^{\mu\nu}$, where $T_{\mu\nu}$ is the stress energy-momentum tensor of the matter component. Hence, such a choice of the corresponding Lagrangian will give rise to $f(R,T_{\mu\nu}T^{\mu\nu})$ theories of gravity. It should be kept in mind that in such a scenario, we are not introducing new forms of non-linear fluid stresses Barrow (1990); Nojiri and Odintsov (2005) like the scalar field, bulk viscosity Barrow (1988), or Chaplygin gas Bento et al. (2002); Gorini et al. (2003). In 2014, Katirci and Kavuk  Katırcı and Kavuk (2014) proposed such a theory for the first time, a covariant generalization of GR which allows the existence of a term proportional to $T_{\mu\nu}T^{\mu\nu}$ in the action functional. Further studies on this theory was carried out  Board and Barrow (2017); Roshan and Shojai (2016), where specific models of this gravity theory have been considered. Roshan et al  Roshan and Shojai (2016) analyzed the possibility of a bounce at early times within the energy momentum squared gravity (EMSG) model, with the specific functional given by $f(R,T^{2})=R+\eta T^{2}$, where $\eta$ is a constant. Board and Barrow Board and Barrow (2017) studied the cosmology of the energy momentum powered gravity (EMPG) model which is a generalization of the EMSG theory, where the model is characterized by $f(R,T^{2})=R+\eta(T^{2})^{n}$, with $\eta$ and $n$ constant parameters. Non exotic matter wormholes are studied in the framework of EMSG in Moraes and Sahoo (2018), and possible constraints from neutron stars were discussed in Akarsu et al. (2018a). Furthermore, recent studies Nari and Roshan (2018); Akarsu et al. (2018b); Keskin (2018) have considered various cosmological applications of the energy momentum squared gravity theory. It has been shown  Liu et al. (2016) that the quantum fluctuations associated to the metric tensor can produce additional cross terms between the Ricci and the energy-momentum tensor. In the framework of energy momentum squared gravity theories, the late time acceleration of the Universe have been investigated in Akarsu et al. (2018c), considering the case of a pressure-less fluid. In this specific case the authors Akarsu et al. (2018c) constrained different parameters of the corresponding model by relying on various values of the Hubble parameter.

From the above handful of literature it is clear that $f(R,T_{\mu\nu}T^{\mu\nu})$ gravity along with its EMSG and EMPG specializations need more attention and hence motivations to study such theories are quite high. As both the concepts of dark energy and modified gravity are till date inadequate to properly explain the observations, it is obvious that there is a lot of room for improvement in both the sectors and there is an open invitation to work and improve both the territories. This work is one such novel attempt to improve upon our existing knowledge of modified gravity theories. The present work will focus on some generalizations of the EMSG model of modified gravity and will try to explore its features via a dynamical system analysis.

The paper is organized as follows: Section II deals with cosmology of energy momentum squared gravity theories. In section III we give a detailed analysis of the dynamical system for two specific case of gravity functions. Finally the paper ends with the corresponding concluding remarks in section IV.

II Energy-momentum Squared Cosmology

The action of our model can be written as Board and Barrow (2017)

[TABLE]

where $f$ is a function depending on the square of the energy-momentum tensor $\mathbf{T^{2}}=T^{\mu\nu}T_{\mu\nu}$ and the scalar curvature $R$. Here, $\kappa^{2}=8\pi G$ and $S_{\rm m}$ represents the action corresponding to the matter component.

If we vary the action with respect to the metric we arrive at the following field equations

[TABLE]

where $\Box=\nabla_{\mu}\nabla^{\mu}$, $f_{R}=\partial f/\partial R$, $f_{\mathbf{T^{2}}}=\partial f/\partial\mathbf{T^{2}}$ and

[TABLE]

where $T$ is the trace of the energy-momentum tensor. By taking covariant derivatives with respect to the field equation (2), one finds the following conservation equation

[TABLE]

As one can see from the above equation that in general, the conservation equation does not hold for this theory. If one chooses $f(R,\mathbf{T}^{2})=2\alpha\log(\mathbf{T}^{2})$, one gets the same result reported in Akarsu et al. (2019).

In the following, we will concentrate on the flat FLRW cosmology for this model whose metric is described by

[TABLE]

with $\delta_{ik}$ being the Kronecker symbol and $a(t)$ the scale factor. Let us now consider that the matter content is described by a standard perfect fluid with $T_{\mu\nu}=(\rho+p)u_{\mu}u_{\nu}+pg_{\mu\nu}$ with $u_{\mu}$ being the 4-velocity and $\rho$ and $p$ are the energy density and the pressure of the fluid respectively. Using these, the energy-momentum tensor gives us $T^{2}=\rho^{2}+3p^{2}$. Further, let us assume $L_{\rm m}=p$ which allows us to rewrite $\Theta_{\mu\nu}$ defined in eqn. (3) as a quantity which does not depend on the function $f$, namely Board and Barrow (2017)

[TABLE]

The modified FLRW equations which corresponds to this particular action are given by

[TABLE]

where dots denote differentiation with respect to the cosmic time $t$, $H=\dot{a}/a$ is the Hubble parameter, and

[TABLE]

was defined. The conservation equation (4) reads as follows

[TABLE]

Clearly, the standard conservation equation does not hold in $f(R,\mathbf{T^{2}})$ cosmology for an arbitrary function. If one chooses $f(R,\mathbf{T^{2}})=f(R)$, all the terms on the RHS of the above equation are zero and the standard conservation equation is recovered. In the following, a standard barotropic equation of state will be assumed:

[TABLE]

where $w$ is the equation of state parameter. Using this relation one gets that

[TABLE]

and then the conservation equation (10) becomes

[TABLE]

In the next section, the dynamical system of the general model will be found and then some specific forms of the function will be adopted to analyse its cosmological properties by using dynamical system techniques. Hereafter, we will use geometric units such that $\kappa^{2}=1$.

III Dynamical system

In this section we will derive the general form of the dynamical system for the modified FLRW equations described by Eqs. (7)-(8). Let us first introduce the following dimensionless variables

[TABLE]

Using these quantities one finds that the first Friedmann equation given by (7) becomes

[TABLE]

The dynamical system for these five dimensionless variables become

[TABLE]

where $N=\log(a)$ and we have further defined the following parameters

[TABLE]

On the other hand, the conservation equation (13) can be written as follows ($w\neq-1$ and $w\neq-1/3$)

[TABLE]

Then, one can further notice the following extra relationships

[TABLE]

and then rewrite the dynamical system (16)-(20) as follows

[TABLE]

where we have also used the Friedmann constraint (15) and the conservation equation (22). Furthermore, in terms of dimensionless and auxiliary intermediate variables, the Friedmann acceleration equation (8) reduces to

[TABLE]

To close the dynamical system, one needs to impose a specific form of the function $f$. By doing this, one needs to derive how the parameters $\Pi$ and $\mho$ are either constants or depend on the dimensionless variables $x_{2},x_{3},x_{4}$ and $x_{5}$. Lastly, it can be seen that this particular choice of dimensionless variables does not necessarily imply that the effective matter density parameter needs to satisfy the usual standard existence conditions, i.e. $0\leq\Omega_{\rm m}^{\rm eff}=x_{4}\leq 1$ due to the appearance of the variation of the gravity functional with respect to the scalar curvature. However, in our analysis we shall consider the necessity of the standard existence conditions as a basic requirement for the validity of the corresponding critical points due to the complexity of the phase space.

Let us finish this section by noting that the effective state parameter is only related to the dimensionless parameter $x_{3}$ as follows

[TABLE]

This expression will be used to understand the nature of the critical points.

III.1 Case 1: $f(R,\mathbf{T^{2}})=f_{0}R^{n}(\mathbf{T^{2}})^{m}$

If we specify the following functional $f(R,\mathbf{T^{2}})=f_{0}R^{n}(\mathbf{T^{2}})^{m}$, where $m,n$ and $f_{0}$ are constant parameters, then we can close the dynamical system. By using the definition of $\Delta$ for this specific functional, we can deduce the following intermediate equation:

[TABLE]

which represents an interrelation between the definition of $\Delta$ in the conservation equation and the intermediate variables $\mho$ and $\Pi$. Moreover, in this specific case, one can note that we have an interrelation between the second and third dimensionless variables, namely,

[TABLE]

Furthermore, for this specific model, one can write the following:

[TABLE]

where we have used the constraint (15). Thus, we have three equations (22), (30) and (32) for the variables $\mho,\Delta$ and $\Pi$ that can be solved in terms of the dimensionless variables $x_{2},x_{4}$ and $x_{5}$, yielding

[TABLE]

where for simplicity, we have defined the quantity

[TABLE]

Therefore, the dynamical system (24)-(27) is reduced to a 3 dimensional one, given by

[TABLE]

where again, for simplicity, we have defined the function

[TABLE]

In the following, we will only focus on the dust case $(w=0)$. In this case, the structure of the phase space consists of three critical points which are valid from a physical point of view described in Table 1. Our analysis considered the requirement of the standard existence conditions which implies that the critical points are located in the real space, while the effective matter density parameter is non–negative and satisfies $0\leq\Omega_{\rm m}^{\rm eff}\leq 1$. The latter requirement represents an extra condition Carloni et al. (2016) added to the existence conditions due to the complexity of the phase space, by generalizing the usual conditions considered in the minimal coupling case. As can be observed from the Table 1, the structure of the phase space is sensitive to the values of the $m$ and $n$ parameters. Furthermore, the stability criteria for each critical point is analyzed in detail by determining the corresponding eigenvalues for each type of solution, constraining the possible values of the $m$ and $n$ parameters from a dynamical point of view. Due to the high complexity of the corresponding eigenvalues for the $P_{2}$ and $P_{3}$ critical points, the analysis relies on numerical evaluations which are described in figures. Each critical points have the following features:

• •

Point $P_{1}$: The first critical point $P_{1}$ represents a solution dominated by the geometrical dark energy component, with physical effects from the shape of the energy momentum squared gravity functional $f(R,\mathbf{T^{2}})$ and its variation with respect to the matter energy momentum tensor $\mathbf{T^{2}}$. In this case the matter density parameter $\Omega_{\rm m}^{\rm eff}$ is equal to zero, corresponding to an cosmological era dominated completely by the geometrical dark energy component. This critical point always exist, the only constraint come from the requirement that the solution has a non–zero denominator, e.g. resulting in: $(m+2)n^{2}+(2m-3)n-m+1\neq 0$ and also $3m(n(n+2)-1)+6n^{2}-9n+3\neq 0$. Using the definition of the effective equation of state in Eq. (29) one can find that at this critical point, the dynamical features depends heavily on the values the $m,n$ parameters, considering also the Eq. (31):

[TABLE]

The stability of this solution is also affected by the $m$ and $n$ parameters, as can be seen from its corresponding eigenvalues:

[TABLE]

The characterization of this critical point is displayed in the Fig. 1 where we have shown the physical effects of the $m$ and $n$ parameters, obtaining some possible constraints which corresponds to the stable case. The effective equation of state is sensitive to the variation of the $m$ and $n$ constants and the solution can manifest a large spectrum of dynamical eras, starting from a super–accelerated expansion to matter, stiff and super–stiff cosmological epochs. Due to the large spectrum of physical epochs present, we have displayed in the right panel of Fig. 1 the variation of the effective equation of state only in regions which corresponds to intervals of interest for modern cosmology. However, for some values of the $m$ and $n$ parameters, this point can represent either a super-accelerating late-time attractor or a saddle matter dominated point.

• •

Point $P_{2}$: The second critical point $P_{2}$ represents a solution characterized by the physical effects coming from the stress energy momentum coupling, the shape of the functional $f$, and matter energy density. In this case, the effective matter density is equal to

[TABLE]

showing a sensitivity to the values of the $m$ and $n$ constants. Analyzing the stability criteria in this case, the expression of the corresponding eigenvalues are too cumbersome to be written in the manuscript. Hence, we shall analyze the physical implications of the $P_{2}$ solution by relying only on numerical evaluations. The numerical description of the second critical point is shown in Fig. 2 where it can be seen that this critical point represents an epoch characterized by an accelerated or super–accelerated expansion. Note that in this case we have taken into account the standard existence condition for the effective matter density parameter which implies $0\leq\Omega_{\rm m}^{\rm eff}\leq 1$. Hence, depending on the values of the $m$ and $n$ parameters, we can obtain a stable solution characterized by an accelerated expansion, in agreement with the current evolution of the known Universe.

• •

Point $P_{3}$: If we take into account the standard existence conditions, then the last critical point $P_{3}$ represents a solution which is characterized by a decelerated expansion of the Universe. As in the previous case, the expressions of the corresponding eigenvalues are too complex to be written here. Hence, in our analysis we rely only on numerical evaluations, determining the dynamical behavior in some limited regions. From the stability analysis, we have observed that this particular solution cannot be stable or pure unstable, and corresponds to a saddle dynamical behavior. Taking into account the existence conditions, we have displayed in Fig. 3 the variation of the effective equation of state in this case for a limited region in the $m-n$ space. This solution corresponds to an epoch with a decelerated expansion, the evolution of the universe is highly sensitive to the values of the $m$ and $n$ parameters. Hence, this solution can explain the dust and radiation cosmological epochs in the evolution of the Universe. We have observed that for some regions the effective matter density parameter can be very close to zero and the effective equation of state can mimic a radiation behavior closely in spite of the absence of the radiation fluid.

The phase space analysis at infinity is performed in the Appendix A.1, where we show that at infinity only some of the critical points are physically viable due to the divergences of the effective equation of state.

Finally for this model we can note that the $f(R,\mathbf{T^{2}})=f_{0}R^{n}(\mathbf{T^{2}})^{m}$ gravity type represents an interesting cosmological model whose viability is very sensitive to the values of the $m$ and $n$ parameters. From a dynamical analysis it can recover the accelerated expansion era and can solve the dark energy problem without introducing a cosmological constant $\Lambda$. Moreover, depending on the values of $m$ and $n$ parameters this model can recover the matter dominated epoch and act towards a stiff fluid solution. As an example, Fig. 4 shows a model for the specific case where $m=2$ and $n=1.0001$, which gives $f(R,\mathbf{T^{2}})=f_{0}R^{1.0001}(\mathbf{T^{2}})^{2}$. This figure depicts the evolution of the effective equation of state for this model. One can notice that this case roughly describes the main epochs of our Universe, starting from a radiation dominated era with $w_{\rm eff}=1/3$, the passing to a matter dominated era with $w_{\rm eff}=0$ for a small interval, then facing an accelerating behaviour to finalizing in a super-accelerating era, the attractor of the corresponding model. As can be noted, the first model curiously exhibits the crossing of the phantom divide line boundary.

III.2 Case 2: $f(R,\mathbf{T^{2}})=\alpha R^{n}+\beta\mathbf{(T^{2})}^{m}$

In this section we will analyze another power-law model considering two power-law terms, the first one coming from the curvature scalar and the second one from the squared of the energy momentum tensor, implicitly given by $f(R,\mathbf{T^{2}})=\alpha R^{n}+\beta\mathbf{(T^{2})}^{m}$ where $\alpha,\beta,m$ and $n$ are constants. We will further assume that $m\neq 0$, $n\neq 0$ and $n\neq 1$. For this model, it is possible to find that the dimensionless variable $x_{2}$ is related to the others ones, namely

[TABLE]

Using this equation, one can straightforwardly find from (21) and (22) that the auxiliary variables become

[TABLE]

Hence, the dynamical system is closed and becomes a 3-dimensional one. From (25)–(27), one gets that the final form of the dynamical system for this model becomes

[TABLE]

This dynamical system contains eight critical points, but only seven satisfy the condition $1\leq x_{4}=\Omega_{\rm m}^{\rm eff}\leq 0$. We wil now concentrate in the dust case ($w=0$). These critical points contain different types of cosmological scenarios depending on the parameters $m$ and $n$. The values of these critical points among with their effective state parameter and their acceleration conditions are displayed in Table 2. In Table 3 is displayed the stability criteria for each point. Let us first describe the main properties of the critical points to then analyse them one by one. The first four critical points $P_{1},..,P_{4}$ are points which are governed by the geometrical dark energy components coming from $f(R,\mathbf{T^{2}})$ gravity. The origin of the phase space, the critical point $P_{1}$, represents a radiation era and depending on the parameters can describe a saddle point. Moreover, the critical points $P_{2}$ and $P_{7}$ also represent a radiation dominated era since their effective state parameter is equal to $1/3$. Depending on the parameters, these points can have different stability properties, but there is a family of parameters for $m$ and $n$ which ensures that all these three critical points which represent a radiation dominated era are either saddle points or unstable points. The critical points $P_{3}$ and $P_{4}$ are purely geometric terms and exhibit similar cosmological behaviours, which can represent acceleration or not depending on the parameters. These points can represent a matter dominated era. The critical point $P_{6}$ always represent an accelerating scenario with a De-Sitter expansion, behaving as a cosmological constant. Each point can be summarised as follows

• •

Point $P_{1}$: The origin of the phase space exists for any values of the parameters $m$ and $n$ and always represents a radiation era. Since $w_{\rm eff}=1/3$, there is no acceleration for this critical point. This point can be either a saddle point or unstable.

• •

Point $P_{2}$: This point corresponds to a universe governed only by the $\mathbf{T}^{2}$ term and represents a radiation era with no acceleration. Its behaviour is similar to the point $P_{1}$ and it also cannot be stable and always exists.

• •

Point $P_{3}$: This critical point represents a universe fully governed by the curvature term and depending on the parameter $n$, can describe different eras and also different stability properties. For example, it can describe an attractor accelerating cosmological solution for $m>1/2\ \land(\frac{1}{2}<n<1\lor 2n\geq\sqrt{3}+1\lor 2n+\sqrt{3}\leq 1)$ or can describe a saddle matter (dust) dominated era if $\left(12n+\sqrt{73}=7\land m<\frac{1}{2}\right)\lor\left(12n=\sqrt{73}+7\land m\neq\frac{1}{2}\right)$. In addition, for $\frac{1}{2}<n<1\lor n>2$ can describe a super-accelerating universe since $\dot{H}=H^{2}(x_{3}-2)>0$. These solutions are also sometimes dubbed as “crossing the phantom divide line”. Moreover, when one adds the extra condition that $m>\frac{1}{2}$, this super-accelerating solution is stable. However, this point cannot represent inflation since it is not possible to find a combination of $m,n$ such that one gets an early accelerating repeller.

• •

Point $P_{4}$: This critical point has similar properties as $P_{3}$, i.e., it can describe accelerating solutions being stable or can describe a matter dominated era being a saddle point, or even can describe super-accelerating solutions depending on the parameters $m$ and $n$. The difference of this critical point with respect to $P_{3}$ is that $P_{4}$ does not represent a universe only governed by a curvature term since a contribution from the $\mathbf{T}^{2}$ is present.

• •

Point $P_{5}$: This critical point can be either stable or a saddle point and can represent a non-dust matter dominated era with $w_{\rm eff}>0$ for $0<m<\frac{1}{2}\land\frac{4m}{2m+1}\leq n<1$. Then, it follows that can also represent a (saddle/stable) stiff matter era when $0<m\leq\frac{1}{6}\land n=\frac{1}{2}$. It cannot represent a standard dust matter dominated era since we already assumed that $n\neq 1$. Due to the condition $0\leq\Omega_{\rm m}^{\rm eff}\leq 1$, this point cannot represent accelerating solutions either.

• •

Point $P_{6}$: This critical point represents a De-Sitter accelerating solution since $w_{\rm eff}=-1$. For different set of family of values of $(m,n)$, it can represent a late-time accelerating attractor.

• •

Point $P_{7}$: This critical point corresponds to a radiation era that can be either stable or saddle depending on the parameters.

The left panel of Fig. 7 shows the regions where the critical points $P_{3},..,P_{7}$ are stable (or spiral stable) where one can see that there is some overlapping regions where different points are stable, but there it is not possible that all of them can be stable. The right panel of Fig. 7 shows the regions where the critical points $P_{3},P_{4}$ and $P_{6}$ are stable are representing accelerating solutions. Thus, this model can reproduce different cosmological eras. For example, it can describe dark energy without a cosmological constant, or a super-accelerating era, or even a stiff dominated era.

Since the phase space is non-compact for this model, one also needs to check if the critical points at infinity are physical or not. See Appendix A.2 for more details about the method used for this. One can see that there are six critical points at infinity, but four of them have a divergent $w_{\rm eff}$, hence, they are not physical. The other two critical points represent a radiation era that cannot be stable. As a numerical example, the Fig. 8 shows the evolution of the effective state parameter for the model $f(R,\mathbf{T^{2}})=\alpha R^{2}+\beta(\mathbf{T^{2}})^{2}$ which corresponds to $m=n=2$. This model exhibits the three main eras of the Universe, starting form a saddle point with $w_{\rm eff}=1/3$, then a matter era with $w_{\rm eff}=0$, finalizing with a late-time de-Sitter behavior.

IV Conclusions

In this paper we have studied the dynamical features of a recent model of modified gravity theory which is known as the energy-momentum squared gravity model $f(R,\mathbf{T^{2}})$. In this model $R$ represents the scalar curvature and $\mathbf{T^{2}}$ the square of the energy-momentum tensor defined as $\mathbf{T^{2}}=T^{\mu\nu}T_{\mu\nu}$. After obtaining the field equations we then concentrated in the study of its cosmology for the standard flat FLRW spacetime. Then, we have analyzed the physical implications of such a modified gravity theory by employing the linear stability method for two specific $f(R,\mathbf{T^{2}})$ functions. The first scenario studied is represented by a model which takes into account a direct product between the scalar curvature $R$ and the energy momentum squared $\mathbf{T^{2}}$ at different real powers $n$ and $m$, with the corresponding gravity function defined as $f(R,\mathbf{T^{2}})=f_{0}R^{n}(\mathbf{T^{2}})^{m}$. In this representation we have assumed that $f_{0},m$ and $n$ are constant parameters. By introducing the dimensionless variables we have been able to represent the dynamics of such a modified gravity theory in flat FLRW as an autonomous system, determining the associated critical points and the corresponding stability properties. In this specific case we have observed that the phase space has a reduced complexity, with a high sensitivity to the values of the $m$ and $n$ parameters. As seen in the analysis, this specific case can recover different cosmological epochs, depending on the physical properties of the critical points obtained.

The first critical point is a general critical point which can recover any cosmological era, staring from an accelerated expansion to matter and radiation behavior, stiff or super–stiff behavior, depending on the values of the $m,n$ constants. The second critical point is associated to an accelerated expansion, the effective matter density and total equation of state showing a high dependence on the values of the $m,n$ constants. In this case, we have determined in the figures provided specific constraints to the values of the $m,n$ constants due to the physical existence conditions and dynamical features. Finally, the third critical point represents a cosmological epoch characterized by a decelerated expansion, an era of a limited interest in the present cosmology.

The second cosmological scenario analyzed corresponds to a different mathematical model of the energy-momentum squared gravity model $f(R,\mathbf{T^{2}})$ which takes into account the following decomposition $f(R,\mathbf{T^{2}})=\alpha R^{n}+\beta\mathbf{(T^{2})}^{m}$. In this case the parameters $\alpha,\beta,m$ and $n$ are assumed to be constants. This scenario shows a higher complexity of the phase space features, with various epochs corresponding to radiation, matter–dominated, de–Sitter, and solutions having accelerating or super–accelerating expansions. As noted in the manuscript, in the second cosmological scenario the values of the $m,n$ parameters dictate the phase space features and the corresponding dynamical properties associated. Furthermore, in our analysis we have obtained different numerical or relational constraints to the values of the $m$ and $n$ due to the existence conditions and physical features of the critical points corresponding to the second cosmological scenario.

The energy-momentum squared gravity theory $f(R,\mathbf{T^{2}})$ represents a recent proposal which takes into account the embeddedness of the energy-momentum squared scalar $\mathbf{T^{2}}=T^{\mu\nu}T_{\mu\nu}$ and have been studied using different approaches Board and Barrow (2017); Moraes and Sahoo (2018); Akarsu et al. (2018a, b); Keskin (2018). Our study is based on the linear stability method and represents a complementary analysis of the energy-momentum squared gravity theory by investigating the phase space features and the stability properties of the critical points, associated to various cosmological epochs. The analysis presented here showed that the energy-momentum squared gravity theory represents an interesting modified gravity model which can explain the current evolution of the Universe and the emergence of the accelerated expansion as a geometrical physical effect, a viable solution to the dark energy problem.

Acknowledgements.

S.B. is supported by Mobilitas Pluss N*∘* MOBJD423 by the Estonian government. P.R. acknowledges Inter University Centre for Astronomy and Astrophysics (IUCAA), Pune, India, for awarding Visiting Associateship. M.Marciu would like to thank M.C. for support and suggestions.

Appendix A Critical points at infinity

When a phase space is not compact, one needs to study if there are critical points at infinity. To do this, one can use the method available in Bahamonde et al. (2018); Bahamonde and Wright (2015); Leon and Fadragas (2012) by introducing compactified Poincaré variables. We follow this approach for the two models studied in this paper. For a 3-dimensional dynamical system with dimensionless variables $x,y$ and $z$, one can use the following Poincaré variables to compactify the phase space:

[TABLE]

where $r^{2}=x^{2}+y^{2}+z^{2}$ and then define $\rho=r/\sqrt{1+r^{2}}$, so that $\rho^{2}=X^{2}+Y^{2}+Z^{2}$. In these coordinates, the dynamics at infinity is recovered when $\rho\rightarrow 1$. It is then simpler to further introduce spherical coordinates

[TABLE]

where $0\leq\theta\leq 2\pi$, $0\leq\rho\leq 1$ and $0\leq\psi\leq\pi$.

A.1 Case 1: $f(R,\mathbf{T}^{2})=f_{0}R^{n}(\mathbf{T^{2}})^{m}$

Following (53) and (54), for this model one has to replace $x=x_{2}$, $y=x_{4}$ and $z=x_{5}$. In this case after the transformation of the dynamical equations  (24)-(27) using the compactified Poincaré variables, we shall consider further the spherical coordinate system at infinity. Due to the high complexity of the dynamical equations before attempting to perform the limits at infinity one needs to choose specific values for the $m$ and $n$ parameters. In the case where $m=3$ and $n=2$ we obtain the following relations at infinity in the leading terms:

[TABLE]

showing that for this specific model the angular part is decoupled. At infinity the critical points are obtained by determining the angularity of the dynamical system in the case where the right hand side of the evolution relations in the limit $\rho\rightarrow 1$ reduces to zero. In this case we have obtained the following critical points in the Poincaré variables $(X_{2},X_{4},X_{5})$ for the general case:

[TABLE]

The first two critical points at infinity, $P_{\infty}^{\pm 1}$ represent a radiation dominated epoch with a zero effective matter density parameter. The second solution at infinity, $P_{\infty}^{\pm 2}$ does not present viable cosmological features at infinity for different values of $\psi$ due to the divergence of either the total equation of state or the effective matter density parameter. However, the last critical points $P_{\infty}^{\pm 3}$ in the case where $\theta=0$ reduces to the first critical points, $P_{\infty}^{\pm 1}$ where the geometrical dark energy component dominates and mimics a radiation era.

A.2 Case 2: $f(R,\mathbf{T^{2}})=\alpha R^{n}+\beta\mathbf{(T^{2})}^{m}$

Following (53) and (54), for this model one has that $x=x_{3}$, $y=x_{4}$ and $z=x_{5}$. By transforming the dynamical system (50)-(52) into Poincaré variables, one gets that at the limit $\rho\rightarrow 1$ (infinity), the dynamical system for the leading terms becomes

[TABLE]

where primes denote differentiation with respect to $N=\log a$. One can notice that the angular part decouples. If one sets the right hand side of these equations equal to zero, one finds three sets of critical points in Poincaré variables $(X_{3},X_{4},X_{5})$, namely,

[TABLE]

Since $0\leq x_{4}\leq 1$, one has that $0\leq(1-\rho^{2})^{-1/2}\rho\sin\theta\sin\psi\leq 1$, and then all these critical points are in the phase space. From (29) one can see that the effective state parameter in the Poincaré variables is

[TABLE]

This quantity is divergent for $P_{\infty,\pm 2}$ and $P_{\infty,\pm 3}$, hence, these critical points at infinity are non-physical unless $\theta=0$ and $\psi=\pi/2$ for $P_{\infty,\pm 2}$ and $P_{\infty,\pm 3}$, respectively. However, for this special choice, the Jacobian of the transformation also diverges. Then, one can conclude that only the the critical points $P_{\infty,\pm 1}$ are physical. For these critical points, the effective state parameter is equal to $1/3$ which represents a radiation era. By going back to the dynamical system for the Poincaré variables, one gets that the Eigenvalues evaluated at $P_{\infty,\pm 1}$ are

[TABLE]

Then, $P_{\infty,\pm 1}$ cannot be stable, it is unstable if $\left(n<1\lor n>\frac{8}{7}\right)\land\left(m<\frac{1}{2}\lor m>\frac{11}{16}\right)$ and saddle otherwise.