Late-time cosmological evolution of a general class of f(R, T) gravity with minimal curvature-matter coupling

TL;DR

This paper investigates the late-time cosmological behavior of a broad class of f(R,T) gravity models with minimal curvature-matter coupling, analyzing stability and evolution using dynamical systems methods.

Contribution

It introduces a dynamical systems approach to study the stability of late-time solutions in f(R,T) gravity with energy-momentum tensor violation.

Findings

Identification of stability conditions for late-time cosmological solutions.

Analysis of how functions g(R) and h(-T) influence cosmic evolution.

Insights into the role of energy-momentum violation in cosmological dynamics.

Abstract

In this work, we study the late-time cosmological solutions of f(R,T)=g(R)+h(-T) models assuming that the conservation of the energy-momentum tensor (EMT) is violated. We perform our analysis through constructing an autonomous dynamical system for the equations of motion. We study the stability properties of solutions via considering linear perturbations about the related equilibrium points. These parameters which are constructed out of the functions g(R) and h(-T) play the main role in finding the late time behavior of the solutions.

Tables

Table 1. Table 1 : The fixed points solutions of f ​ ( 𝖱 , 𝖳 ) = g ​ ( 𝖱 ) + h ​ ( − 𝖳 ) 𝑓 𝖱 𝖳 𝑔 𝖱 ℎ 𝖳 f({\sf R},{\sf T})=g({\sf R})+h(-{\sf T}) gravity.

Point	Coordinates $(x_1,x_2,x_3,x_4,x_5)$	$w^{\text{(eff)}}$	${\Omega }^{\mathrm{𝖽𝗎𝗌𝗍}}$	${\Omega }^{\text{(}\text{DE}\text{)}}$
${𝖯}_1$	$(-4,5-x_4,0,x_4,x_5)$	$\frac{1}{3}$	$-x_5$	$1+x_5$
${𝖯}_2$	$(0,-1-x_4,2,x_4,x_5)$	$-1$	$-x_5$	$1+x_5$
${𝖯}_3$	$(-4,\frac{4\left(n+1\right)\left(5n+9\right)}{{\left(2n+3\right)}^2},0,\frac{4n+9}{{\left(2n+3\right)}^2},x_5)$	$\frac{1}{3}$	$-x_5$	$1+x_5$
${𝖯}_4$	$(0,-\frac{4{\left(m+1\right)}^2{n}^2+\left(m\left(4m+21\right)+14\right)n+3\left(5m+4\right)}{{\left(m+1\right)}^2{\left(2n+3\right)}^2},2,$
	$\frac{-m\left(8mn+9m+3n+3\right)+2n+3}{{\left(m+1\right)}^2{\left(2n+3\right)}^2},x_5)$	$-1$	$-x_5$	$1+x_5$
${𝖯}_5$	$(-1,0,0,0,0)$	$\frac{1}{3}$	$2$	$-1$
${𝖯}_6$	$(\frac{6m\left(n+1\right)}{\left(m+1\right)\left(2n+3\right)},-\frac{m\left(4n+6\right)+n+3}{{\left(m+1\right)}^2\left(2n+3\right)},\frac{m\left(4n+6\right)+n+3}{\left(m+1\right)\left(2n+3\right)},$
	$\frac{-m\left(8mn+9m+3n+3\right)+2n+3}{{\left(m+1\right)}^2{\left(2n+3\right)}^2},$
	$-\frac{2\left(n+1\right)\left(m\left(8mn+9m+3n+3\right)-2n-3\right)}{{\left(m+1\right)}^2{\left(2n+3\right)}^2})$	$-\frac{m\left(2n+3\right)+1}{\left(m+1\right)\left(2n+3\right)}$	$0$	$1$
${𝖯}_7$	$(-\frac{2\left(n+3\right)}{2n+3},0,0,\frac{4n+9}{{\left(2n+3\right)}^2},\frac{2\left(n+1\right)\left(4n+9\right)}{{\left(2n+3\right)}^2})$	$\frac{1}{3}$	$0$	$1$
${𝖯}_8$	$(\frac{3m}{m+1},-\frac{4m+1}{2{\left(m+1\right)}^2},2-\frac{3}{2\left(m+1\right)},0,0)$	$-1+\frac{1}{m+1}$	$\frac{2-m\left(8m+3\right)}{2{\left(m+1\right)}^2}$	$\frac{m\left(10m+7\right)}{2{\left(m+1\right)}^2}$

Table 2. Table 2 : The fixed points solutions of f ​ ( 𝖱 , 𝖳 ) = 𝖱 α + ( − 𝖳 ) β 𝑓 𝖱 𝖳 superscript 𝖱 𝛼 superscript 𝖳 𝛽 f({\sf R},{\sf T})={\sf R}^{\alpha}+(-{\sf T})^{\beta} gravity.

Point	Coordinates $(x_1,x_2,x_4)$	$w^{\left(\mathrm{𝖾𝖿𝖿}\right)}$	${\Omega }^{\mathrm{𝖽𝗎𝗌𝗍}}$	${\Omega }^{\text{(}\text{DE}\text{)}}$
${𝖯}^{\left(𝖬\right)}$	$(\frac{3m}{m+1},-\frac{4m+1}{2{\left(m+1\right)}^2},0)$	$\frac{m\left(10m+7\right)}{2{\left(m+1\right)}^2}$	$\frac{2-m\left(8m+3\right)}{2{\left(m+1\right)}^2}$	$\frac{1}{m+1}-1$
${𝖯}^{\left(𝖲\right)}$	$(0,\frac{2}{m+1},\frac{2}{m+1}-1)$	$-1$	$\frac{2\left(m-1\right)\left(n+1\right)}{m+1}$	$\frac{-2mn-m+2n+3}{m+1}$
${𝖯}^{\left(\mathrm{𝖣𝖤}\right)}$	$(\frac{6m\left(n+1\right)}{\left(m+1\right)\left(2n+3\right)},-\frac{m\left(4n+6\right)+n+3}{{\left(m+1\right)}^2\left(2n+3\right)},$
	$\frac{-m\left(8mn+9m+3n+3\right)+2n+3}{{\left(m+1\right)}^2{\left(2n+3\right)}^2})$	$-\frac{m\left(2n+3\right)+1}{\left(m+1\right)\left(2n+3\right)}$	$0$	$1$

Equations

$$\displaystyle S=\int\sqrt{-g}d^{4}x\left[\frac{1}{2\kappa^{2}}f\Big{(}{\sf R},{\sf T}\Big{)}+{\sf L}^{\rm{(m)}}\right],$$

$${\mathsf{T}}_{\mu \nu }\equiv -\frac{2}{\sqrt{-g}}\frac{\delta \left[\sqrt{-g}{\mathsf{L}}^{(\mathrm{m})}\right]}{\delta g^{\mu \nu }}.$$

$${\mathsf{T}}_{\mu \nu }=-\rho \frac{d{\mathsf{L}}_m}{d\rho }{\mathsf{U}}_{\mu }{\mathsf{U}}_{\nu }+\left[{\mathsf{L}}_m-\rho \frac{d{\mathsf{L}}_m}{d\rho }\right]g_{\mu \nu },$$

$${\mathsf{T}}_{\mu \nu }=\rho {\mathsf{U}}_{\mu }{\mathsf{U}}_{\nu }.$$

$$\displaystyle F({\sf R},{\sf T}){\sf R}_{\mu\nu}-\frac{1}{2}f({\sf R},{\sf T})g_{\mu\nu}+\Big{(}g_{\mu\nu}\square-\triangledown_{\mu}\triangledown_{\nu}\Big{)}F({\sf R},{\sf T})=$$

$$\displaystyle\Big{(}\kappa^{2}-{\mathcal{F}}({\sf R},{\sf T})\Big{)}{\sf T}_{\mu\nu}-\mathcal{F}({\sf R},{\sf T})\mathbf{\Theta_{\mu\nu}},$$

$${\mathbf{\Theta }}_{\mu \nu }\equiv g^{\alpha \beta }\frac{\delta {\mathsf{T}}_{\alpha \beta }}{\delta g^{\mu \nu }},$$

$$\mathcal{F}(\mathsf{R},\mathsf{T})\equiv \frac{\partial f(\mathsf{R},\mathsf{T})}{\partial \mathsf{T}}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{mbox}and\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}F(\mathsf{R},\mathsf{T})\equiv \frac{\partial f(\mathsf{R},\mathsf{T})}{\partial \mathsf{R}}.$$

$$\displaystyle ds^{2}=-dt^{2}+a^{2}(t)\Big{(}dr^{2}+r^{2}d\Omega^{2}\Big{)},$$

$$\displaystyle 3H^{2}F({\sf R},{\sf T})+\frac{1}{2}\Big{(}f({\sf R},{\sf T})-F({\sf R},{\sf T}){\sf R}\Big{)}+3\dot{F}({\sf R},{\sf T})H=$$

$$\displaystyle\Big{(}\kappa^{2}+{\mathcal{F}}({\sf R},{\sf T})\Big{)}\rho,$$

$$\displaystyle 2F({\sf R},{\sf T})\dot{H}+\ddot{F}({\sf R},{\sf T})-\dot{F}({\sf R},{\sf T})H=-\Big{(}\kappa^{2}+{\mathcal{F}}({\sf R},{\sf T})\Big{)}\rho,$$

$$3H^{2}F(\mathsf{R},\mathsf{T})=8\pi G\left(\rho +{\rho }^{\text{(DE)}}\right),$$

$$\displaystyle-2\dot{H}F({\sf R},{\sf T})=8\pi G\Big{(}\rho+\rho^{\textrm{({\sf DE})}}+p^{\textrm{({\sf DE})}}\Big{)},$$

$$\displaystyle 8\pi G\rho^{\textrm{({\sf DE})}}\equiv{\mathcal{F}}({\sf R},{\sf T})\rho^{({\sf m})}-3\dot{F}({\sf R},{\sf T})H-\frac{1}{2}\Big{(}f({\sf R},{\sf T})-F({\sf R},{\sf T}){\sf R}\Big{)},$$

$$\displaystyle 8\pi Gp^{\textrm{({\sf DE})}}\equiv\ddot{F}({\sf R},{\sf T})+2\dot{F}({\sf R},{\sf T})H+\frac{1}{2}\Big{(}f({\sf R},{\sf T})-F({\sf R},{\sf T}){\sf R}\Big{)}.$$

$${\dot{\rho }}^{\text{(DE)}}+3H\left(1+w^{\text{(DE)}}\right){\rho }^{\text{(DE)}}=0.$$

$$w^{\text{(eff)}}={\Omega }^{\text{(DE)}}{w}^{\text{(DE)}},$$

$$f(\mathsf{R},\mathsf{T})=g(\mathsf{R})+h(-\mathsf{T}).$$

$$1+\frac{g}{6H^2{g}^{\prime }}+\frac{h}{6H^2{g}^{\prime }}-\frac{\mathsf{R}}{6H^2}+\frac{\dot{g^{\prime }}}{H{g}^{\prime }}=\frac{{\kappa }^2\rho }{3H^2{g}^{\prime }}-\frac{h^{\prime }\rho }{3H^2{g}^{\prime }},$$

$$2\frac{\dot{H}}{H^2}+\frac{\ddot{g^{\prime }}}{H^2{g}^{\prime }}-\frac{\dot{g^{\prime }}}{H{g}^{\prime }}=-\frac{{\kappa }^2\rho }{H^2{g}^{\prime }}+\frac{h^{\prime }\rho }{H^2{g}^{\prime }}.$$

$$x_1\equiv -\frac{\dot{g}(\mathsf{R})}{H{g}^{\prime }(\mathsf{R})},$$

$$x_2\equiv -\frac{g(\mathsf{R})}{6H^2{g}^{\prime }(\mathsf{R})},$$

$$x_3\equiv \frac{\mathsf{R}}{6H^2}=\frac{\dot{H}}{H^2}+2,$$

$$x_4\equiv -\frac{h(-\mathsf{T})}{6H^2{g}^{\prime }(\mathsf{R})},$$

$$x_5\equiv \frac{\mathsf{T}{h}^{\prime }(-\mathsf{T})}{3H^2{g}^{\prime }(\mathsf{R})},$$

$${\Omega }^{\mathsf{dust}}\equiv \frac{{\kappa }^2\rho }{3H^2{g}^{\prime }(\mathsf{R})},$$

$$w^{\text{(eff)}}=\frac{1}{3}\left(1-2x_3\right).$$

$${\Omega }^{(\mathsf{DE})}=x_1+x_2+x_3+x_4+x_5,$$

$${\Omega }^{\mathsf{dust}}+{\Omega }^{(\mathsf{eff})}=1.$$

$$m\equiv \frac{\mathsf{R}{g}^{\prime \prime }(\mathsf{R})}{g^{\prime }(\mathsf{R})},$$

$$r\equiv -\frac{\mathsf{R}{g}^{\prime }(\mathsf{R})}{g(\mathsf{R})}=\frac{x_3}{x_2},$$

$$n\equiv -\frac{\mathsf{T}{h}^{\prime \prime }(-\mathsf{T})}{h(-\mathsf{T})},$$

$$s\equiv -\frac{\mathsf{T}{h}^{\prime }(-\mathsf{T})}{h(-\mathsf{T})}=\frac{x_5}{2x_4},$$

$$\frac{1}{2}\mathcal{F}{\nabla }_{\mu }\mathsf{T}+{\mathsf{T}}_{\mu \nu }{\nabla }^{\mu }\mathcal{F}=0,$$

$$\dot{\mathcal{F}}(\mathsf{R},\mathsf{T})=\frac{3}{2}H\mathcal{F}(\mathsf{R},\mathsf{T}).$$

Full text

1 Introduction

One of the recent major challenges in modern cosmology is the late time behavior of the Universe. Using distant type Ia supernovae as standard candles, it has been discovered that the Universe is currently experiencing an accelerated expansion [1, 2, 3, 4]. Such a discovery led to the widespread acceptance of the idea that the Universe is dominated by a mysterious substance (with negative pressure) named “dark energy” (DM) that drives this accelerating expansion. It is a general belief that the concordance model [5] which includes simply cosmological constant within the Einstein’s field equation is also well consistent with the most observational data. Therefore, with this assumption for which all analyses of the recent data have been performed, the equation of state parameter (EoS) of DE is found to be $w^{\rm{({\sf DE})}}=−1.006\pm 0.045$ [4]. However, there are some substantial difficulties within the cosmological constant description of late time behavior of the Universe. The well-known cosmological constant problem states that the concordance model suffers from a fine-tuning issue related to its energy scale, if it is attributed to vacuum energy [6, 7]. Indeed, from a theoretical viewpoint, the necessity of explaining DE, as well as the second dominant component of the Universe, i.e., dark matter (DM) [8], raises the fundamental question of whether the general theory of relativity (GR) could provide a befitting quantitative description of the Universe on all scales. In general, there are two basic approaches that could provide a setting for theoretical explanation of the accelerated expansion of the Universe. The first procedure consists of modifying the matter content of the Universe by introducing a DE sector, starting either with a canonical scalar field, a phantom field, or the combination of both fields in a unified model and proceeding afterwards to more complicated scenarios; see [9] and references therein for more details and reviews. The second approach is to modify the gravitational sector itself (see e.g. [10, 11, 12, 13]), motivated by this fundamental assumption that at large astrophysical and cosmological scales, standard GR may not describe correctly the dynamical evolution of the Universe. To deal with this issue, several efforts have been made, among which gravity theories extending GR have attracted a great deal of interest over the past decades. In the context of modified gravity theories (MGT), a geometric description for DM can be provided and the accelerated expansion can be achieved at late times, thus, the problem of cosmological constant may be resolved (for a detailed dynamical system analysis of several cosmological models in the context of MGT see [14].) Various methods have been proposed so far, in order to modify the gravitational action (cf. [15] for a consummate review), giving rise to different classes of MGT. From a historical perspective, in going beyond the GR theory, the first step has been to generalize the geometric part of the Einstein-Hilbert action. One of the simplest forms with the aim of modifying the dynamical behavior of gravity and consequently reaching cosmological scenarios different to GR, is to consider functions of the Ricci scalar ${\sf R}$, dubbed $f({\sf R})$ theories [11]. Recently, this picture has been extended to incorporate the trace of energy momentum tensor (EMT) within the gravitational sector of the action. Such a class of MGT known as “$f({\sf R},{\sf T})$ gravity theories” which allow for matter-geometry coupling are significant as they can provide, from a fundamental theoretical perspective, a complete theoretical description for the late time acceleration of the Universe, without resorting to the existence of DE [16]. These models can also offer some alternative explanations for the nature of DM [17, 18] and the problem of flattening of the galaxy rotation curves [19]. Since the advent of this idea, different features of $f({\sf R},{\sf T})$ gravity have been investigated and in particular many efforts have been devoted to studying the cosmological aspects of this theory [16, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 31, 32, 33, 34] and very recently, quantum cosmology of $f({\sf R},{\sf T})$ gravity has been proposed in [35] where the authors have investigated the evolution of the wave function of early Universe by virtue of the Wheeler-de Witt equation. In [20, 22], the background evolution of the Universe has been studied under the assumption that the EMT is conserved. Taking into account the conservation of EMT as a basic constraint, the authors have shown that $f({\sf R},{\sf T})$ models with a minimal coupling between the geometrical and matter sectors yield the late time accelerated expansion of the Universe. However, most of these models lead to the present value for the EoS parameter as $w^{\rm{({\sf DE})}}=-1/2$, which is not observationally acceptable. Therefore, it is reasonable to seek alternative ways in order to cure this problem.

Very recently, an interesting result has been reported in [36] where the authors have found that in the context of unimodular gravity, violation of EMT conservation can lead to the emergence of an effective cosmological constant. Regardless of how the conservation of EMT can be violated,111For example, in the context of classical GR, the violation of EMT conservation could provide a setting to deal with the entropy problem of standard cosmology, see e.g. [37] and references therein. Such a violation may also occur by modifying quantum mechanical models and also some proposals based on the causal set approach to quantum gravity [36]. it is a novel idea to reconsider previous unsuccessful models under this assumption. It should be emphasized that, though the violation of EMT conservation may not always lead to a physically reasonable accelerated expansion scenario at late time, it can be a novel task to investigate the effects of violation of EMT conservation on the late time behavior of a gravitation theory. Motivated by this idea, we study minimally coupled models of the form $f({\sf R},{\sf T})={\sf R}+\Lambda({\sf T})$ as presented in [31], ignoring the usual relation for EMT conservation. These types of models can account for, at least, modification of the Einstein-Hilbert action by putting aside the assumption that the matter Lagrangian plays a subordinate and passive role only. In the present paper we extend the previous work to include models with arbitrary functions of the Ricci curvature scalar. We shall see that the violation of EMT conservation could give rise to acceptable accelerated expansion at late times, contrary to the situations where the conservation of EMT is respected. Note that in $f({\sf R})$ gravity, applying the Bianchi identity to the geometrical parts of the field equation leads to a null result and, thus, the condition on EMT conservation cannot be relaxed. From this point of view, $f({\sf R},{\sf T})$ gravity has more chance for designing viable cosmological models.

Up to now, the dynamical system approach has been widely used to describe cosmological aspects of different gravity theories [14]. Using this powerful method, one can find cosmological solutions and investigate the nature of the fixed points and their stability properties. However, this approach has not been completely utilized in the framework of $\sf{f(R,T)}$ gravity. Nevertheless, dynamical system techniques have been exploited to investigate cosmological behavior of models of the form $\sf{f(R,T)=R+C\sqrt{T}}$ [20, 22], models that deal with the Einstein static Universe [25] and late time solutions for $\sf{f(R,T)=R+h(T)}$ gravity models [30, 31]. The present work could be among the first studies which deal with the dynamical system approach in the context of $\sf{f(R,T)=g(R)+h(T)}$ models with the aim to describe the late time cosmological acceleration. We would like to add that there are different works in the literature where the authors have considered cosmological acceleration in ${\sf f(R,T)}$ gravity, however, they have employed algebraic manipulations or have obtained cosmological solutions by applying indirect procedures such as reconstructing ${\sf f(R,T)}$ models for a presumed cosmological behavior [38, 39, 40, 41, 42, 43, 44, 45]. A significant trait of our work is a parametrization of the variables which is used in the dynamical field equations. We shall show that in $f(\sf{R},\sf{T})=g(\sf{R})+h(\sf{T})$ gravity models, one can recast the functions constructed out of curvature and the trace of EMT, i.e., $g(\sf{R})$ and $h(\sf{T})$ (at least for well-defined ones) into the four dimensionless parameters $(m,r)$ and $(n,s)$, respectively. More exactly, one can rewrite the expressions given for $g(\sf{R})$ and $h(\sf{T})$ in terms of equivalent dimensionless $m(r)$ and $n(s)$ functions. Therefore, when these parameters appear in the dynamical equations, the critical points and their stability properties could determine consistent ranges of these parameters and thus the validity of the underlying model. Most studies in the literature either deal with particular models from the beginning of the problem or present a general form of the field equations and finally employ them for some specific types of $g(\sf{R})$ and $h(\sf{T})$ functions. For instance, authors of [46, 47, 48, 49] have discussed models with $h(\sf{T})=\alpha\sf{T}$. The case of a power law form, $h(\sf{T})=\sf{T}^{\alpha}$, has been considered in [50, 51] and investigation of some non-trivial forms has been carried out in [45, 52, 53]. Briefly speaking, we use a formalism which works at least for a wide well-defined kinds of Lagrangians. This formalism was introduced for the first time in [54].

This paper is arranged as follows: In Sect. 2, the field equations for $f({\sf R},{\sf T})$ gravity together with essential variables will be introduced. In this section we also discuss the violation of EMT conservation. In Sect. 3 we construct an autonomous system of equations and consider the corresponding fixed point solutions. Moreover, we study the cosmological behavior of the minimal $f({\sf R},{\sf T})$ models by relaxing the conservation of EMT. In Sect. 4, we explore models with power law form for the Lagrangian. Finally in Sect. 5, we summarize our results.

2 Field equations of $f({\sf R},{\sf T})$ gravity

In the present section, we review the field equations of $f({\sf R},{\sf T})$ gravity and introduce necessary variables which shall be used in the rest of paper. Since we are interested in studying the late time solutions of a general class of $f({\sf R},{\sf T})=g({\sf R})+h({\sf T})$ models, we assume only the pressure-less fluid as the matter content. We call this particular model the “minimal $f({\sf R},{\sf T})$ model” for, only the minimal coupling between geometrical and matter sectors is taken into account. $f({\sf R},{\sf T})$ gravity theories have been initially introduced by the following action [16]

[TABLE]

where ${\sf R}$, ${\sf T}\equiv g^{\mu\nu}{\sf T}_{\mu\nu}$, ${\sf L}^{\rm{(m)}}$ are defined as the Ricci curvature scalar, the trace of energy momentum tensor and the Lagrangian of pressure-less matter (here, pressure-less fluid), respectively. $g$ is the determinant of the metric, $\kappa^{2}\equiv 8\pi G$ and we set $c=1$. The EMT of matter, ${\sf T}_{\mu\nu}$ is defined as the Euler-Lagrange expression of the matter Lagrangian, i.e.,

[TABLE]

A point here which begs a few elucidation is the inclusion of trace of EMT within the curvature part of the action. This might, at first sight, be seen a little ambiguous as one usually expects that the matter Lagrangian or a function of it must be coupled to an arbitrary function of the Ricci scalar. In this respect, non-trivial gravitational Lagrangians have been proposed where the gravitational sector directly interacts with matter sector through a non-minimal coupling process, the so-called $f({\sf R},{\sf L}_{m})$ gravity [55]. Various choices for matter Lagrangian would lead to different gravitational theories and thus have various cosmological scenarios as the outcome. In this regard, one possible way to take into account the footprint of matter within the gravitational Lagrangian is to include the trace of EMT within the gravitational sector of the action.222 The dependence of the gravitational action on the trace of EMT may be due to the presence of quantum effects (conformal anomaly), or of some exotic imperfect quantum fluids; see e.g, [17] for more details. This can be regarded as a possible way of defining the gravitational Lagrangian (1) where one assumes that the contributions due to matter that appear within the argument of the function $f$ solely correspond to the EMT constructed out of the material part of the Lagrangian, as defined via (2). We then might intuitively imagine that the matter existing within the spacetime fabric (as descried by the Lagrangian ${\sf L}_{m}$) creates curvature. The matter content then interacts minimally with spacetime curvature through the metric and curvature mutually has an effect on the matter distribution within the spacetime. However, the EMT trace is taken a priori as the only agent for direct interaction between matter and curvature. The a priori contribution owing to matter in an uncommon interaction with spacetime curvature may also have some connections with known subjects such as a geometrical description of the curvature induced due to matter distribution, describing physical forces geometrically, and a geometrical origin for the matter content of the Universe; see, e.g., [56] and the references therein. A specific class of $f({\sf R},{\sf L}_{m})$ gravity theories can be constructed assuming that the matter Lagrangian is a function of energy density of matter $\rho$ only, so that ${\sf L}_{m}={\sf L}_{m}(\rho)$. The EMT then reads [57]

[TABLE]

where ${\sf U}_{\mu}$ being the four-vector velocity of the fluid which satisfies the condition ${\sf U}_{\mu}{\sf U}^{\mu}=-1$. Considering, for instance, the choice of ${\sf L}_{m}=-\rho$ for matter Lagrangian [58, 59], we get the EMT for a pressure-less fluid as

[TABLE]

The EMT trace for such a fluid (which is our case study) is found as ${\sf T}=-\rho$, and this could justify the way of the appearance of the EMT trace within the gravitational sector of the action.

Varying action (1) with respect to the metric field leads to the following field equation [16]

[TABLE]

where

[TABLE]

and for the sake of convenience, we have defined the following functions for derivatives with respect to the trace ${\sf T}$ and the Ricci curvature scalar ${\sf R}$, as

[TABLE]

We assume that the spacetime geometry is described by the spatially flat Friedmann–Lemaître–Robertson–Walker (FLRW) metric, which in spherical coordinates is given by

[TABLE]

where $a$ denotes the scale factor of the Universe and $d\Omega^{2}$ being the line element on a unit two-sphere. The field equation (2) then gives rise to

[TABLE]

as the modified Friedmann equation,

[TABLE]

as the modified Raychaudhuri equation and $H$ indicates the Hubble parameter.

In the case of $f({\sf R},{\sf T})$ gravity one may define a DE component which is responsible for the accelerated expansion at the late times. To this aim, the Friedmann equation, (2), and the Raychaudhuri equation (10) can be redefined as [22]333This representation has been already used in the literature for $f({\sf R})$ gravity; for more details we consult the reader to [60, 61, 62, 63, 64, 65, 66, 67, 68].

[TABLE]

and

[TABLE]

where in equations (11) and (12) the first term is related to the pressure-less fluid and $\rho^{\textrm{({\sf DE})}}$ and $p^{\textrm{({\sf DE})}}$ denote the density and pressure of DE component, respectively, which are defined as

[TABLE]

and

[TABLE]

To see the consistency of the above definitions, one can substitute the Lagrangian for GR, i.e., $f({\sf R},{\sf T})={\sf R}$ into equations (13) and (14). This leaves us with the usual Friedman equations instead of those given in (11) and (12), and also $\rho^{\textrm{({\sf DE})}}=0,p^{\textrm{({\sf DE})}}=0$ in the case of GR. In fact, these definitions introduce the DE component as a real fluid444As we shall see, interactions of these two components needs to define an effective behavior and therefore effective pressure and density..

In this case, the EoS parameter for DE is defined as $w^{\textrm{({\sf DE})}}\equiv p^{\textrm{({\sf DE})}}/\rho^{\textrm{({\sf DE})}}$, as usual. One can check that definitions (13) and (14) guarantee the continuity equation for DE component, i.e.,

[TABLE]

In the context of modified gravity theories (specially in $f({\sf R})$ gravity), an effective EoS parameter would be defined very similar to the GR case. In GR, the EoS parameter for a perfect fluid is defined as $w=p/\rho$. Using this definition in modified gravity, one can consider an effective behavior of fluids and correspondingly define an effective EoS as $w^{\textrm{({\sf eff})}}\equiv-1-2\dot{H}/3H^{2}$. Performing some straightforward calculations we conclude that

[TABLE]

where we have defined a density parameter for DE component as $\Omega^{\textrm{({\sf DE})}}\equiv 8\pi G\rho^{\textrm{({\sf DE})}}/3H^{2}F$. The effective EoS includes information about gravitational interactions for different components (pressure-less and DE in our case). More precisely, one can define an effective density and pressure using $3H^{2}=8\pi G\rho^{\textrm{({\sf eff})}}$ and $-2\dot{H}=8\pi G(\rho^{\textrm{({\sf eff})}}+p^{\textrm{({\sf eff})}})$, thereby, it is possible to describe the cosmological consequences of the theory with the help of the effective behavior of the pressure-less and DE components.555Redefinition of the Friedman equations of $f({\sf R},{\sf T})$ gravity in terms of an conserved effective fluid could lead to interesting results. This approach have let us to discover some novel features of $f({\sf R},{\sf T})$ gravity [69].

Next, we proceed to consider a Lagrangian with minimal coupling between the Ricci scalar and the trace of EMT, i.e.,666We do not include a coupling constant, since, such a coefficient could be absorbed into the dynamical system variables.

[TABLE]

Because of our metric signature, we have ${\sf T}=-\rho$ for the pressure-less matter, therefore, to avoid ambiguity due to the negative sign, we hereafter study $h(-{\sf T})$ functions. For this class of $f({\sf R},{\sf T})$ models we obtain $\mathcal{F}({\sf R},{\sf T})=-h^{\prime}(-{\sf T})$ and $F({\sf R},{\sf T})=g^{\prime}({\sf R})$, hence the field equations (2) and (10) yield

[TABLE]

and

[TABLE]

Let us now define a few variables and parameters which are useful for recasting the field equations (18) and (19) into a closed dynamical system. These variables are defined as

[TABLE]

where we have used the expression ${\sf R}=6(\dot{H}+2H^{2})$ for the Ricci scalar within the definition (22). Using definition (22) along with the mentioned definition of EoS, we obtain

[TABLE]

Also, using the definitions (13), (20)-(25) and the density parameter of DE which is defined in (16), we get

[TABLE]

Furthermore, we use the following definitions

[TABLE]

where all primes denote differentiation with respect to the argument. Expressions (29)-(32) show that eliminating ${\sf R}$ from (29) and (30) leads to the function $m=m(r)$ and also eliminating ${\sf T}$ from (31) and (32) gives the function $n=n(s)$. As we mentioned before, parameters $r$ and $m$ have originally been used in [54], in order to describe cosmological solutions of $f({\sf R})$ gravity models. Here, we extend their method to the case of minimal $f({\sf R},{\sf T})$ model. Models with $m=m(r)$ and $n=n(s)$, instead of $g({\sf R})$ and $h(-{\sf T})$, can be very suitable for dynamical system analysis.

In [20, 22], the authors have comprehensively considered cosmological solutions of $f({\sf R},{\sf T})$ gravity when the conservation of EMT is respected. To see the consequences of EMT conservation, one can apply the Bianchi identity to the field equation (2). Then imposing the usual form for the EMT conservation of matter current i.e., $\nabla_{\alpha}{\sf T}^{\alpha\beta}=0$ which leads to $\dot{\rho}+3H\rho=0$, and applying the Bianchi identity to the field equation (2), we arrive at the following simple constraint:

[TABLE]

whereby substituting for the metric components gives

[TABLE]

As is seen, equation (34) places a constraint on the functionality of $f({\sf R},{\sf T})$, which in the case of the minimal coupling directly leads to a constraint on the trace dependent sector (as we shall see below), since, as is well known, in pure $g({\sf R})$ gravity the Bianchi identity is automatically satisfied, therefore leading to the conservation of EMT. Substituting the Lagrangian (17) in constraint (34) gives

[TABLE]

whereby after a straightforward algebra we get

[TABLE]

where $C_{1}$ and $C_{2}$ are constants of integration. Therefore, constraint (34) determines the form of $f({\sf R},{\sf T})$ function and in the minimally coupled case this constraint gives the relation (36). We note that this solution with the use of definition (31) corresponds to $n=-1/2$. The behavior of these types of models has been fully studied in [20, 22] and further issues discussed in [31] using a dynamical system approach. It is shown that despite the appearance of accelerated expansion at late times, the corresponding dark energy cannot be suitably matched with the observational data. In these models a DE with $w^{\rm{({\sf DE})}}_{0}=-1/2$ is responsible for the accelerated expansion of the Universe. Hence, it may be a good idea to explore the cosmological behavior of the model under violation of EMT conservation. To this aim, we rewrite the field equation (2) for specific models given by the choice (17) for $f({\sf R},{\sf T})$ function, as follows:

[TABLE]

Then, we apply the Bianchi identity to the field equation (37) whence we obtain the following covariant conservation equation

[TABLE]

where the argument of $h^{\prime}(-{\sf T})$ has been dropped for abbreviation. Note that equation (38) is more general than equation (33), which holds when we do not assume the conservation of matter EMT, i.e., $\nabla_{\alpha}{\sf T}^{\alpha\beta}=0$. Substituting for the components of EMT of a dust into equation (38), we obtain

[TABLE]

where we have used $\rho=-{\sf T}$. Comparing equation (35) and (39) shows that, in the case of relaxing the conservation of EMT, a more complicated constraint would determine the form of $h(-{\sf T})$ function. Once the function $h(-{\sf T})$ is determined, the dependency of $-{\sf T}$ and consequently $\rho$ on the scale factor can be calculated. On the other hand, in two situations equation (39) could simply lead to the usual conservation equation, i.e., $\dot{\rho}+3H\rho=0$; when $h(-{\sf T})=0$ that is, in GR case and $g({\sf R})$ gravity. The next situation could occur when the two expressions in parentheses are equal, giving then the solution (36). Otherwise, there is a modified version of the conservation equation. Therefore, equation (39) completely depicts the behavior of the trace/matter density with respect to the scale factor. After determining ${\sf T}(a)$ via using equation (39), one can solve for the equations of motion (2) and (10) to obtain the scale factor and then discuss the cosmological consequences. More precisely, we can rewrite equation (39) as follows:

[TABLE]

where ${\sf T}_{0}$ and $a_{0}$ are the present values. Principally, the above relation can give the trace as a function of the scale factor, however, the left hand side integral may not easily be solved in general. In other words, even if the integration is carried out, finding the density as an explicit function of scale factor may not be possible, at least analytically. Such a situation can occur when the EoS is not taken simply as that of a perfect fluid. However, from (40) for $h(-{\sf T})=\chi^{2}(-{\sf T})^{\alpha}$ and a dust fluid we obtain

[TABLE]

where $C$ is a constant of integration. As is seen, solution (41) may not be further simplified for arbitrary values of $\alpha$. For applications in the dynamical system approach which we shall utilize in later sections, we rewrite equation (39) in terms of the dimensionless variables, as follows:

[TABLE]

Equation (42) reduces to the standard form for $n=-1/2$, which corresponds to the conserved model (36). It is noteworthy to briefly point out an important thermodynamical trait of conservation equation (39). As the author of [17] has beautifully detailed this issue, one can write down the conservation equation in $f({\sf R},{\sf T})$ gravity in order to explain the matter creation process (we do not enter the details of calculations). That is, the following equations can be obtained:

[TABLE]

where $\Gamma$ is the particle creation rate. This equation denotes no particle creation when $\mathcal{F}=0$, e.g., in the case of $f({\sf R},{\sf T})=g({\sf R})+\Lambda$ gravity, where $\Lambda$ is a cosmological constant. We can interpret this modified conservation equation as indicating an irreversible matter creation process caused by an irreversible energy flow from the gravitational field to the matter constituents. Briefly, in the context of $f({\sf R},{\sf T})$ gravity, it is possible that the spacetime transmutes to matter. Using the concepts of quantum cosmology, the particles which are created in this way may be in the form of some scalar particles (bosons) which specify the DM content of the Universe. One can also rewrite the above conservation equation as follows:

[TABLE]

The above form of conservation equation suggests a new term that is called the creation pressure, $p_{\textrm{c}}$ [70]. From thermodynamic point of view, this means that the coupling of geometry and matter generates (or can be translated as) a new pressure effect (which may even be considered as an effective bulk viscosity [17]) that participates in a (gravitational) particle creation procedure. Therefore, unlike some other models (such as GR and $g({\sf R})$ gravity) within which the problem of particle creation must be considered based on an admissible physical processes, in the framework of models which allow for coupling between matter and geometry (such as $f({\sf R},{\sf T})$ gravity) the irreversible thermodynamic processes are completely determined by the gravitational action. In the next section, we shall use equation (42) in order to construct the dynamical system representation of equations (18) and (19) and study its cosmological consequences.

3 Late-time behavior in $f({\sf R},{\sf T})=g({\sf R})+h({\sf T})$ gravity

In this section we study cosmological consequences of models of the type $f({\sf R},{\sf T})=g({\sf R})+h({\sf T})$ allowing for the EMT conservation to be violated, that is we suppose $\nabla_{\alpha}{\sf T}^{\alpha\beta}\neq 0$. Equations (18) and (19) in terms of dimensionless variables (20)-(32) are obtained:

[TABLE]

respectively. Therefore, relations (42)-(44) would help us to obtain an autonomous system of equations of motion for dimensionless variables $x_{1}-x_{5}$. Evolutionary equations are then obtained as follows:

[TABLE]

As can be seen, both $n$ and $m$ parameters appear within the above system of equations, that is, fixed point solutions and their properties depending on the parameters ($n$ and $m$) which exhibit the nature of the underlying model. In Table 1 we have summarized equilibrium points and their cosmological features for the above system.

Before discussing the solutions and their stability properties, there are some issues that should be pointed out. Note that, since the radiation fluid is absent in the present study, we exclude fixed points with $w^{({\sf eff})}=1/3$ from our considerations. As Table 1 indicates, there are some fixed points for which the coordinates and cosmological features are not completely determined. More precisely, fixed point ${\sf P}_{2}$ is the solution of a set of equations (45)-(49) for arbitrary values of $x_{4}$ and $x_{5}$ and ${\sf P}_{4}$ is a solution for arbitrary values of $x_{5}$. Naturally, the question may arise of whether it is possible to determine exactly the physical properties of the fixed points. The answer is positive! There are some mathematical and physical criteria that must be met. Since at the critical points we have $dx_{i}/dN=0,\leavevmode\nobreak\ i=1,\cdots,5$, in general, this condition must hold for every function, namely $\mathcal{W}(x_{1},x_{2},x_{3},x_{4},x_{5})$. For instance, $r=r(x_{2},x_{3})$ and $s=s(x_{4},x_{5})$ must be stationary at the fixed point. This means that

[TABLE]

and

[TABLE]

whereby using definitions (21)–(24) and (29)–(32), we have

[TABLE]

and

[TABLE]

where we have defined

[TABLE]

As a result, solutions must accept the conditions $m=-r-1$ (provided that $m(r)\neq 0$), $r=0$ or $x_{1}=0$ and $n=s-1$ or $s=0$. The condition $r=x_{3}/x_{2}=0$ leads to $x_{3}=0$ and condition $s=x_{5}/2x_{4}=0$ gives $x_{5}=0$. This fact restricts the value of $x_{5}$ to null, and with this choice the location of point ${\sf P}_{4}$ in the phase space will be completely determined. Moreover, ${\sf P}_{2}$ and ${\sf P}_{4}$ must have $x_{5}=0$ in order to completely represent a dominant de Sitter phase. Therefore, one of the novel features developed by violation of EMT conservation is the appearance of a new de Sitter solution (${\sf P}_{4}$) as compared to the case of pure $f({\sf R})$ gravity, for which the properties depend upon the matter sector of the Lagrangian.777$f({\sf R})$ gravity theories accept a de Sitter solution in the late times with coordinates $(0,-1,2)$ as authors of [54] have shown. Another interesting result is that a DE fixed point (${\sf P}_{6}$) solution is achieved within this framework. This point indicates a dominant DE solution with an effective equation of state which depends on both the geometrical and the matter part of the Lagrangian via the parameters $m$ and $n$. Let us now proceed with exploring the stability properties of the fixed points ${\sf P}_{2}$, ${\sf P}_{4}$, ${\sf P}_{6}$ and ${\sf P}_{8}$.

• •

The equilibrium points ${\sf P}_{2}$ and ${\sf P}_{4}$

These two points have the same eigenvalues given as

[TABLE]

Therefore, these points are representative of stable equilibrium points if $m$ and $n$ parameters satisfy the following intervals888Indeed, the last two eigenvalues given in (55) will become negative for $\frac{16}{25}\leqslant m<1$ signaling that these two points are stable nodes. While, for $0<m<\frac{16}{25}$, the last two eigenvalues become complex with pure negative real part for which the corresponding point is called stable focus in the literature. However, these type of fixed points are asymptotically stable.

[TABLE]

The stability region for these solutions has been presented in gray color in the right panel of Figure 1. I is seen that, for every curve lying within the gray region, ${\sf P}_{2}$ or ${\sf P}_{4}$ will represent a late time solution.

• •

The equilibrium Point ${\sf P}_{6}$

The eigenvalues for this equilibrium point are obtained as

[TABLE]

where

[TABLE]

where $m^{\prime}\equiv dm(r)/dr$ and $n^{\prime}\equiv dn(s)/ds$ must be calculated at the desired fixed point. There are two class of solutions in which the fixed point ${\sf P}_{6}$ is a stable point. The first class corresponds to those for which ${\sf P}_{6}$ is a stable node and those that ${\sf P}_{6}$ is a spiral source. We have plotted the stability region for this point in the left panel of Figure 1.

Point ${\sf P}_{6}$ is the only critical point which lies both on the $m=-r-1$ and $n=s-1$ curves, that is for this fixed point we have $x_{3}/x_{2}=r=(-1-m)$ and $x_{5}/2x_{4}=s=n+1$. This means that, ${\sf P}_{6}$ is a DE solution only for the roots of equations $n(s_{u})=s_{u}-1,\leavevmode\nobreak\ u=1,...,N$ and $m(r_{v})=r_{v}-1,\leavevmode\nobreak\ v=1,...,\mathcal{N}$. In other words, if for a model specified in terms of the functions $m(r)$ and $n(s)$, there exit a class of $s_{u}$ and $r_{v}$ roots, we then can conclude that the model would accept ${\sf P}_{6}$ as a DE solution at late times. Besides, conditions on $n^{\prime}$ or $m^{\prime}$ as shown in the left panel of Figure 1 must be satisfied for these roots. For example, for a model admitting an $s_{u}$ and $r_{v}$ solution (in blue colored region in Figure 1) for which either $n^{\prime}(s_{u})<1,m^{\prime}(r_{v})>-1$ or $n^{\prime}(s_{u})>1,m^{\prime}(r_{v})<-1$ holds, a DE solution consistent with Plank 2015 data exists.

• •

The equilibrium Point ${\sf P}_{8}$

By considering the matter density and EoS parameters for this fixed point, we find out that it gives a matter dominated solution in the limit $m\to 0$. The eigenvalues for this point are obtained as

[TABLE]

This fixed point is unstable under the following conditions

[TABLE]

We note that, in order to avoid a short matter dominated era or even the absence of it, we must chose the limit $m\to 0^{+}$. For $m\to 0^{-}$, the second eigenvalue becomes a large positive number and hence, in this case ${\sf P}_{8}$ indicates a very transient matter dominated era. Numerical manipulations show that, in order to have a proper matter era, that is for $m\to 0^{+}$, the condition $m^{\prime}>-1$ must be satisfied, otherwise, matter dominated solutions would be disappeared. Actually, numerical studies reveal that, for models with a ${\sf P}_{6}$ solution lying within the red area in the left panel of Figure 1, there is no matter dominated solution. There is only a direct transition from either of the critical points with $w^{({\sf eff})}=1/3$ to ${\sf P}_{2},{\sf P}_{4}$ or ${\sf P}_{6}$.

Mathematically accepted cosmological solutions are those that include a true connection between the matter dominated fixed point and late time solutions which correspond to de Sitter points ${\sf P}_{2}$ or ${\sf P}_{4}$ or the DE fixed point ${\sf P}_{6}$. Thus, the possible cosmological solutions are classified as follows:

• •

Solutions which correspond to trajectories from ${\sf P}_{8}$ with $m^{\prime}(r\to-1^{-})>-1$ to ${\sf P}_{2}$ or ${\sf P}_{4}$ with $n>0$ or $n<-1/2$.

For a certain model which is specified by its $m(r)$ and $n(s)$ functions, physically justified transitions in phase space occur for those conditions which make ${\sf P}_{8}$ an unstable fixed point and ${\sf P}_{2}$ or ${\sf P}_{4}$ a stable one, simultaneously. From the right panel of Figure 1, there are two stable regions for ${\sf P}_{2}$ and ${\sf P}_{4}$; The regions are defined by intervals $0<m(r_{v})<1,\leavevmode\nobreak\ n(s_{u})<-1/2$ and $0<m(r_{v})<1,\leavevmode\nobreak\ n(s_{u})>0$ for the roots $r_{v}$ and $s_{u}$. Since, the models must be designed in such a way that the condition $m^{\prime}(r\to-1^{-})>-1$ holds, all de Sitter solutions can be connected to matter ones. In the region with $n(s)<0$, cosmological solutions with a de Sitter phase at late times can be classified as two distinct categories; a) initial values can be set such that ${\sf P}_{6}$ solution would not be satisfied. In this case, there is only a direct transition from matter to de Sitter solution. b) However, initial values can also be chosen so that an unstable ${\sf P}_{6}$ solution be realized. In such cases an intermediate transient ${\sf P}_{6}$ phase for $-0.5<n<0$999See the white narrow band in the right panel of Figure 1. would appear. To illustrate these two types of solutions we have drawn in Figure 2, the effective EoS parameter for a class of models with $g({\sf R})={\sf R}\log[\alpha{\sf R}]$ which corresponds to $m(r)=r+1/r$, and $h(-{\sf T})=c_{2}{\rm exp}{[-c_{1}(-{\sf T})-\beta(-{\sf T})^{4}/12]}$ giving $n(s)=s-\beta/s$, for $\beta=0.75$ and two different sets of initial values for the variables $x_{1}-x_{5}$101010The effective equation of state in figure 2 is independent of the constants $\alpha$, $c_{1}$ and $c_{2}$.. In this example, we always have $m^{\prime}(r)>-1$ and $n^{\prime}(s)>1$ leading to an unstable ${\sf P}_{6}$ fixed point which is recognized by the root $s_{v}=\beta=0.75$ for which $n(0.75)=-0.25$ and correspondingly gives $w^{\rm{({\sf DE})}}_{0}=-0.4$. Thus, the $n(s)$ curve evolves into $n<-1/2$ or $n>0$ regions in which the de Sitter fixed point is stable (see blue curve in the right panel of Figure 1).

For models with $n>0$ there is a proper cosmological solution corresponding to a direct transition from matter to de Sitter phases.

• •

Solutions which correspond to trajectories from ${\sf P}_{8}$ with $m^{\prime}(r\to-1^{-})>-1$ to ${\sf P}_{6}$.

Figure 1 shows that a necessary stability condition for ${\sf P}_{6}$ is $-3/2<n<0$. However, as it is already emphasized, there is no proper matter dominated solution in the region where $-1<n<0$. Models with $n<-1$ has been already accepted as de Sitter solution. This is because the curve $n(s)$ for all values approaching to a ${\sf P}_{6}$ root in the region $-3/2<n(s_{v})<-1$, lies within the stable de Sitter range. Briefly speaking, there is no proper cosmological solution with a DE phase at late times which corresponds to a ${\sf P}_{6}$ solution.

There are some classes of minimal models for which the above theme slightly differs; models that each of gravitational or matter sectors or both of them within the Lagrangian are in the form of power law. In these cases, $m^{\prime}$ or $n^{\prime}$ or both vanish, which means that the pairs ($m,r$) or ($n,s$) or both of them are constant. Hence, since $x_{3}=rx_{2}$ and $x_{5}=2sx_{4}$, the number of independent variables would reduce from five to four or three. This fact leads to some changes in the fixed point solutions and their stability of the dimensionally reduced dynamical system. To illustrate this issue, we shall investigate the class of models for which both parts of their Lagrangian are power law. This is the subject of next section.

4 Pure minimal case $f({\sf R},{\sf T})={\sf R}^{\alpha}+(-{\sf T})^{\beta}$

In this section, we briefly study the late time solutions of pure minimally coupled models of the type $f({\sf R},{\sf T})={\sf R}^{\alpha}+(-{\sf T})^{\beta}$. For these models we obtain $r=x_{3}/x_{2}=-\alpha$ and $s=x_{5}/2x_{4}=\beta$. That is, for these cases, there are only three independent variables $x_{1}$, $x_{2}$ and $x_{4}$. Hence, this dimensional reduction demands slightly different method. Here, because $m^{\prime}=0$ and $n^{\prime}=0$, the stability properties of fixed points are described by the $m$ and $n$ parameters only. For these types of models there are six fixed point solutions of which three are of physical interest, as shown in Table 2.

It is interesting to note that a de Sitter and a DE solution have appeared again. The first fixed point can denote the matter dominated era only for $m=(\alpha-1)\to 0^{+}$ which leads to $\alpha\to 1^{+}$. Note that the eigenvalues of this point are given by the first three parts of expression (• ‣ 3), that is, the matter dominated point is well behaved only for $m\to 0^{+}$ and is unstable only when $n<0$ as the other two eigenvalues are negative for $m\to 0^{+}$. The stability conditions for de Sitter point ${\sf P}^{{\sf(S)}}$ are the same as (56) which means that only for $n<-1/2$, the matter dominated point ${\sf P}^{{(\sf M)}}$ can be connected to de Sitter point. Since the matter density of a de Sitter fixed point depends on the model parameters, this fact implies that to have a thoroughly de Sitter dominant phase (i.e., $\Omega^{\sf dust}=0$ and $\Omega^{\rm{({\sf DE})}}$=1) we must have $n\to-1$ for arbitrary non-zero $m$. Therefore, the model with $m\to 0^{+}$ and $n=-1$ includes a true connection from a long enough matter dominated era to a stable de Sitter era at late times. The value $n\to-1$ corresponds to $\beta=0$, from which in this case we have a cosmological constant instead of the function $h(-{\sf T})$.

The eigenvalues for ${\sf P}^{{}^{\rm{({\sf DE})}}}$ are given by the first three parts given in expression (• ‣ 3) which show that this point is stable within the same region as shown in Figure 1, regardless of the the conditions on $m^{\prime}$ and $n^{\prime}$. Therefore, all models with $m>0$ and $-3/2<n<0$ could provide a true connection from matter dominated to DE fixed points. However, the blue colored zone in Figure 1 restricts models to those that are consistent with the observational data. This means that, only models with $n\to-1$ are acceptable. To show the behavior of such a model with two late time solutions we have plotted in Figure 3, the phase space diagrams for two different values $n=-0.7,-1.8$ with $m=+0.001$ in $(x_{2},x_{4})$ plane. The red point denotes the de Sitter fixed point, black one shows the matter dominated point and the green point indicates the DE one. The purple trajectories show physically justified connections. As can be seen, the DE point is unstable for $n<-3/2$ which makes ${\sf P}^{{\sf(S)}}$ as the only late time solution. For $-3/2<n<-1/2$ both ${\sf P}^{{\sf(S)}}$ and ${\sf P}^{{}^{\rm{({\sf DE})}}}$ are stable, nevertheless, for $-1<n<-1/2$ the trajectories will be trapped by ${\sf P}^{{}^{\rm{({\sf DE})}}}$ (see Figure 3) and for $n<-1$ will be attracted by ${\sf P}^{{}^{\sf(S)}}$. For $n=-1$ the two points coincide. A particular case has occurred; the coordinates of the DE fixed point and its eigenvalues (see the results given in (• ‣ 3)) are singular for $n=-3/2$. In [31], we have discussed in detail the cosmological behavior of this case by algebraic treatments and have shown that in this case the late time solution corresponds to a de Sitter era. As a final remark, we point out that, from 2015 Planck data [4] which has determined the present value of the EoS of DE as $-1.051<w_{0}^{\rm{({\sf DE})}}<-0.961$, we may conclude that power law models can be admissible only for $m\to 0^{+}$ and $-1.025<n<-0.980$ which corresponds to $-0.024<\alpha<0.02$. Therefore, these types of modified gravity theories, by including two different late time solutions, deserve more investigations as alternatives of the $\Lambda$CDM model.

5 Concluding remarks

In this work, in the context of a particular class of $f({\sf R},{\sf T})$ gravity theories, we have discussed the issue of accelerated expansion of the Universe at the late time cosmological regimes under assumption that the EMT conservation is violated. Though the assumption on relaxing the EMT conservation may not necessarily provide a late time accelerated expansion scenario, we tried to seek for possible relations between such an assumption and the cosmic late time acceleration within the framework of $f({{\sf R},{\sf T}}$) gravity. Specifically, our attention was concentrated on those models for which the Lagrangian can be chosen with minimal coupling between the geometrical and matter sectors, which we call minimally coupled models and which can be written as $f({\sf R},{\sf T})=g({\sf R})+h({\sf T})$. We utilized the dynamical system approach in order to reformulate the field equations. We used a way of parameterizations of models via introducing the parameters $m={\sf R}g^{\prime\prime}({\sf R})/g^{\prime}({\sf R})$ and $r=-{\sf R}g^{\prime}({\sf R})/g({\sf R})$ for the geometrical part of the Lagrangian and $n=-{\sf T}h^{\prime\prime}(-{\sf T})/h(-{\sf T})$ and $s=-{\sf T}h^{\prime}(-{\sf T})/h(-{\sf T})$ for the matter part. Thus, an $f({\sf R},{\sf T})$ model can be determined by $m(r)$ and $n(s)$ functions. These parameters appear within the coordinates of equilibrium points of the system and also within their corresponding eigenvalues. Therefore, all cosmological properties of the system, as derived through exploring the eigenvalues, will depend upon these parameters. In fact, these parameters play a fundamental role for determining the cosmological viability of $f({\sf R},{\sf T})$ models under relaxing the condition on EMT conservation. The application of this approach to explain the late time cosmological acceleration has been introduced in [54] in the framework of $f(\sf{R})$ gravity models. It has also been used in $f({\sf R},{\sf T})={\sf R}+\alpha{\sqrt{\sf T}}$ as the only case which respects the conservation of EMT [20] and within the framework of $f({\sf R},{\sf T})={\sf R}+\alpha\Lambda(\sf{T})$ gravity [31], to extract cosmological solutions. Therefore, the current study may be regarded as an attempt to complete the previous investigations on the cosmological late time behavior of the minimally coupled $f({\sf R},{\sf T})$ gravity models. The study of the dynamical system representation of the field equations shows that there are three fixed points that can represent the late time behavior of the model; two de Sitter points together with a critical point for which $\Omega^{\rm{({\sf DE})}}=1$ and $\Omega^{\sf{(m)}}=0$ and its effective EoS depends on the mentioned parameters. We call the latter fixed point the DE solution. Having investigated and analyzed the properties of fixed points, we found that three types of transition between matter to either of de Sitter or DE fixed points can occur. Therefore, depending on the underlying model as well as the initial values, there are three classes of cosmological solutions, so that, for all cases, in order to have a prolonged matter dominated era the condition $m(r\lessapprox-1)=0^{+}$ must hold. One can construct models in such a way that a “secure” transition from the matter to a de Sitter solution occurs. These types of models have to satisfy the condition $n>0$. Some models include a DE fixed point at late times if the conditions $-1/2<n<0$, $n^{\prime}>1$ and $m^{\prime}(r\lessapprox-1)<-1$ are fulfilled. However, in these models the effective EoS does not match the observational data. Finally, in some models, before reaching the de Sitter phase there is a transient phase of the DE solution, as well. All models of this class must respect the condition $n<-1/2$. Note that, in the last two cases, one could set the initial values in order to obtain a direct transition to a de Sitter phase. However, care must be taken; since for the de Sitter fixed points the matter and DE densities depend on the fixed point coordinates, the constants of the model must be chosen such that they lead to a dominant de Sitter phase at late times.

In summary, we can conclude that violation of EMT conservation in the minimal $f({\sf R},{\sf T})$ gravity models could provide conditions under which the underlying model includes a de Sitter expansion at late times, at least for some subset of space parameters of the model. This fact indicates that allowing for the violation of EMT conservation, the situation of minimal $f({\sf R},{\sf T})$ gravity theories can be encouraging enough for cosmological models and merits further investigation, in comparison with the same models when EMT conservation is taken for granted. As the authors of [20, 22, 33] have already elucidated, respecting the EMT conservation leads to a DE description at late times but the effective EoS for such models would be found as $w^{({\sf eff})}=-1/2$, which could not be observationally accepted.

We also have considered models with power law sections within the Lagrangian as $f({\sf R},{\sf T})={\sf R}^{\alpha}+(-{\sf T})^{\beta}$. It is found that these types of models can describe the recently accelerated expansion of the Universe for $\alpha\to 1^{+}$ and $-0.024<\beta<0.02$.

As the final remark, we would like to highlight that the problem of de Sitter solution in the late time and more importantly the transition from matter dominated era to de Sitter era has been investigated previously in some other work. However, most was based on a reconstruction method in order to obtain models which include the desired behaviors [38, 39, 40, 42]. Nevertheless, the method presented in the herein work may provide a more comprehensive plan to investigate all late time solutions as well as their stability.