1 Introduction

Among many interesting foundations of the modern cosmology, one of the most controversial of them is the late time acceleration of the Universe [1, 2]. The traditional solution to this observation is reconsidering the cosmological constant term in the Einstein–Hilbert theory [3]. However, the cosmological constant suffers from phenomenological/theoretical issues, persuading modern cosmologists to consider modified gravitational theories [4]. There are three main streams on modifying gravity at large scales, one is considering some extra degrees of freedom, such as scalar field [5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23] or a vector field [24,25,26,27,28,29,30,31,32,33,34,35,36], and the other is to modify the gravitational interaction itself, such as massive gravity [37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61], Weyl-Cartan theories [62,63,64,65,66,67,68,69,70], etc. The third category is to modify the gravity action to contain higher order term in curvature tensor, such as f(R) theories [71,72,73,74].

There is also another line of thought about modifying the Einstein–Hilbert theory, which contains non-minimal interactions between matter field and curvature. As a very simple example, one can add a non-minimal coupling between the Ricci scalar and the matter Lagrangian [75]

$$\begin{aligned} S=\int d^4x\sqrt{-g}\big [\kappa ^2(R-2\Lambda )+f(R,\mathcal {L}_m)+\mathcal {L}_m\big ]. \end{aligned}$$
(1)

There are however other possible choices which contains f(R, T) or \(f(R,T,R_{\mu \nu }T^{\mu \nu })\) gravities where T is the trace of the energy–momentum tensor [76,77,78,79,80,81,82,83,84]. Cosmological consequences as well as other concepts of these theories is vastly investigated in the literature. It should be noted that in all of the above theories, the consequence of adding an extra term to the Einstein–Hilbert action is to weaken the role of the cosmological constant, or at an idealized level to nullify its effect.

Recently, a subclass of the non-minimal coupling between matter and geometry attracts some interests, which contains higher order terms in the energy momentum tensor. The action of this theory can be written as [85,86,87]

$$\begin{aligned} S=\int d^4x\sqrt{-g}\big [\kappa ^2(R-2\Lambda )+\alpha (T_{\mu \nu }T^{\mu \nu })^\eta +\mathcal {L}_m\big ], \end{aligned}$$
(2)

where \(\eta \) and \(\alpha \) are constants. It is worth mentioning that the special case \(\eta =2\) is dubbed energy–momentum squared gravity (EMSG). For small and negative values of \(\eta \), the model can satisfy the observational data, at least at background level, without any cosmological constant term. However, it is shown that for other values of the parameter \(\eta \) some non-zero cosmological constant, or another dynamical field is needed [85,86,87].

Another interesting possibility on the non-minimal coupling between matter and geometry would be to consider derivative terms containing matter fields [88,89,90]. If one noted that the matter Lagrangian \(\mathcal {L}_m\) is a scalar field and add some Galileon-like terms to the Einstein–Hilbert action, one would get the accelerated expansion of the Universe without the need for cosmological constant [89, 90]. Also, at the background level, such models satisfy the conservation equation for the energy–momentum tensor, which is promising since other non-minimal matter couplings could not fulfill such a property. The simplest action in this line could be written as

$$\begin{aligned} S=\int d\,^4x\sqrt{-g}\bigg [\kappa ^2 R-\alpha \,\nabla _\mu f\,\nabla ^\mu f\bigg ]+S_m, \end{aligned}$$
(3)

where \(f=f(\mathcal {L}_m)\) is an arbitrary function of the matter Lagrangian. It is shown that one can obtain an accelerated expansion of the Universe for power law function \(f\propto (\mathcal {L}_m)^\eta \), with small and negative \(\eta \) and positive values of \(\alpha \).

In this paper, we will explore a new possibility on the non-minimal coupling of matter and geometry. This contains a coupling between the Ricci tensor and the energy–momentum tensor in the form

$$\begin{aligned} \mathcal {L}\supseteq R_{\mu \nu }T^{\alpha \mu }T_\alpha ^{~\nu }. \end{aligned}$$
(4)

The above interaction term can be seen as a generalization of the EMSG theory [85,86,87]. This is true because one can write the interaction term of the EMSG as

$$\begin{aligned} \alpha g_{\mu \nu }T^{\alpha \mu }T_\alpha ^{~\nu }. \end{aligned}$$

Now considering the fact that in the maximally symmetric space-times we have \(R_{\mu \nu }=\alpha g_{\mu \nu }\), one can write the interaction term of the EMSG in this space-times as

$$\begin{aligned} R_{\mu \nu }T^{\alpha \mu }T_\alpha ^{~\nu }. \end{aligned}$$

The interaction term we are considering in this paper can then be considered as a generalization of the EMSG interaction terms to general space-times. There is also a difference between this term and the term \(RT^2\) where T is the trace of the energy–momentum tensor. The later term is a subset of f(R, T) gravity, which our interaction term does not belong to it.

As we have discussed earlier, the interaction term (4) does not satisfy the conservation of the energy–momentum tensor. In this paper, we will impose the conservation of the energy–momentum tensor by adding a Lagrange multiplier term to the action. The resulting theory would have the effects of the non-minimal matter/geometry couplings in a way that the energy–momentum itself is conserved.

The paper is organized as follows. In the next section we introduce the model and obtain the field equations. In Sect. 3 we will consider the background cosmology of the model and in Sect. 4 we investigate the dynamical system analysis of the model. In Sect. 5 we obtain the dynamical equation for the matter density contrast at first order in perturbation variables. We then use the observational data on the Hubble parameter H and also for \(f\sigma _8\) to estimate the best fit values of the model parameter. We conclude in Sect. 6.

2 The model

Let us consider the action functional of the form

$$\begin{aligned} S&=\int d^4x\sqrt{-g}\big [\kappa ^2(R-2\Lambda )\nonumber \\&\quad +\alpha R_{\mu \nu }T^{\alpha \mu }T_\alpha ^{~\nu }+A_\mu \nabla _\nu T^{\mu \nu }+\mathcal {L}_m\big ], \end{aligned}$$
(5)

where \(\Lambda \) is the cosmological constant, \(\alpha \) is an arbitrary constant with mass dimension \(M^{-6}\), and \(A_\mu \) is a Lagrange multiplier ensuring that the theory respects energy–momentum conservation. It should be noted that, there are strong constraints on the conservation of the matter sources. Here, we have alleviate this problem by introducing a Lagrange multiplier term to insure that the energy–momentum tensor remains conserved. The same procedure have been done before, for example in f(R, T) theories [84]. One can then consider the vector \(A_\mu \) as a modulator potential which determines the required force we need to keep the matter conserved and not convert directly to geometry.

Varying the action (5) with respect to the metric gives

$$\begin{aligned}&\kappa ^2(G_{\mu \nu }+\Lambda g_{\mu \nu })=\frac{1}{2}T_{\mu \nu }+\frac{1}{2}g_{\mu \nu }(\mathcal {L}_m g_{\alpha \beta }\nonumber \\&\quad -T_{\alpha \beta })\nabla ^\alpha A^\beta +\frac{1}{2}A^\alpha \nabla _\alpha T_{\mu \nu }+T_{\alpha (\mu }\nabla _{\nu )}A^\alpha \nonumber \\&\quad +\alpha g_{\mu \nu }\Bigg (T^{\alpha \beta }\nabla _\beta \nabla _\gamma T_\alpha ^{~\gamma }+\frac{1}{2}\nabla _\beta T_{\alpha \gamma }\nabla ^\gamma T^{\alpha \beta }\nonumber \\&\quad +R_{\alpha \gamma \beta \sigma }\Bigg (\mathcal {L}_m g^{\gamma \sigma }T^{\alpha \beta }-\frac{1}{2}T^{\alpha \beta }T^{\gamma \sigma }\Bigg )\Bigg )\nonumber \\&\quad +\alpha \Big (R^\alpha _{~(\mu }T_{\nu )}^{~\beta }T_{\alpha \beta }-2\mathcal {L}_m R^\alpha _{~(\mu }T_{\nu )\alpha }\nonumber \\&\quad -R_{\alpha \beta \gamma (\mu }T_{\nu )}^{~\alpha }T^{\beta \gamma }+R^{\alpha \beta }(T_{\alpha \mu }T_{\beta \nu }-T_{\alpha \beta }T_{\mu \nu })\nonumber \\&\quad +T^\alpha _{~(\mu }\Box T_{\nu )\alpha }-T^{\alpha \beta }\nabla _\beta \nabla _{(\mu }T_{\nu )\alpha }\nonumber \\&\quad +\nabla _\alpha T_{\mu \beta }\nabla _\alpha T_\nu ^{~\beta }-\nabla ^\beta T^\alpha _{~(\mu }\nabla _{\nu )}T_{\alpha \beta }\Big ). \end{aligned}$$
(6)

Also, taking the covariant divergence of the equation of motion (6) will determine the equation of motion of the Lagrange multiplier \(A_\mu \)

$$\begin{aligned}&\frac{1}{2}\mathcal {L}_m A^\alpha R_{\alpha \nu }-\frac{1}{2}A^\alpha R_\alpha ^{~\beta }T_{\nu \beta }+\frac{1}{2}\mathcal {L}_m\Box A_\nu -\frac{1}{2}T_{\nu \alpha }\Box A^\alpha \nonumber \\&\quad +\frac{1}{2}\nabla \alpha A^\beta (\nabla _\nu T_{\alpha \beta }-\nabla _\alpha T_{\nu \beta }-\nabla _\beta T_{\nu \alpha })+\nabla _{[\alpha }\mathcal {L}_m\nabla _{\nu ]}A^\alpha \nonumber \\&\quad +\alpha \Bigg [T^{\beta \gamma }T_\nu ^{~\alpha }\nabla _\alpha R_{\beta \gamma }+2T_{\beta (\alpha }R_{\nu )}^\beta \nabla ^\alpha \mathcal {L}_m\nonumber \\&\quad +\frac{1}{2}T_\nu ^{~\beta }\nabla ^\alpha R(\mathcal {L}_m g_{\alpha \beta }-T_{\alpha \beta })+\mathcal {L}_mT^{\alpha \beta }\nabla _\beta R_{\nu \alpha }+\mathcal {L}_mR^{\alpha \beta }\nabla _\beta T_{\nu \alpha }\nonumber \\&\quad -T^{\beta \gamma }T_\nu ^{~\alpha }\nabla _\gamma R_{\alpha \beta }+2R^{\alpha \beta }\nabla _\gamma T_{\beta [\alpha }T_{\nu ]}^{~\gamma }\nonumber \\&\quad -\frac{1}{2}T^{\alpha \beta }T^{\gamma \delta }(\nabla _\nu R_{\alpha \gamma \beta \delta }+2\nabla _\delta R_{\nu \alpha \beta \gamma })-\nabla _\nu (R_{\alpha \beta }T^{\alpha \beta }\mathcal {L}_m)\Bigg ]=0. \end{aligned}$$
(7)

Varying the action (5) with respect to the Lagrange multiplier \(A_\mu \) gives

$$\begin{aligned} \nabla _\mu T^{\mu \nu }=0. \end{aligned}$$
(8)

One can write the metric field equation as

$$\begin{aligned} \kappa ^2(G_{\mu \nu }+\Lambda _{eff}g_{\mu \nu })=\frac{1}{2}(T_{\mu \nu }+T^{eff}_{\mu \nu }), \end{aligned}$$
(9)

where we have defined the effective cosmological constant and energy–momentum tensor as

$$\begin{aligned} \kappa ^2\Lambda _{eff}&=\kappa ^2\Lambda +\frac{1}{2}(\mathcal {L}_m g_{\alpha \beta }-T_{\alpha \beta })\nabla ^\alpha A^\beta \nonumber \\&\quad +\alpha \Bigg (T^{\alpha \beta }\nabla _\beta \nabla _\gamma T_\alpha ^{~\gamma }+\frac{1}{2}\nabla _\beta T_{\alpha \gamma }\nabla ^\gamma T^{\alpha \beta }\nonumber \\&\quad +R_{\alpha \gamma \beta \sigma }\Bigg (\mathcal {L}_m g^{\gamma \sigma }T^{\alpha \beta }-\frac{1}{2}T^{\alpha \beta }T^{\gamma \sigma }\Bigg )\Bigg ), \end{aligned}$$
(10)

and

$$\begin{aligned} T_{\mu \nu }^{eff}&=A^\alpha \nabla _\alpha T_{\mu \nu }+2T_{\alpha (\mu }\nabla _{\nu )}A^\alpha \nonumber \\&\quad +2\alpha \Big (R^\alpha _{~(\mu }T_{\nu )}^{~\beta }T_{\alpha \beta }-2\mathcal {L}_m R^\alpha _{~(\mu }T_{\nu )\alpha }-R_{\alpha \beta \gamma (\mu }T_{\nu )}^{~\alpha }T^{\beta \gamma }\nonumber \\&\quad +R^{\alpha \beta }(T_{\alpha \mu }T_{\beta \nu }-T_{\alpha \beta }T_{\mu \nu })-T^{\alpha \beta }\nabla _\beta \nabla _{(\mu }T_{\nu )\alpha }\nonumber \\&\quad +T^\alpha _{~(\mu }\Box T_{\nu )\alpha }+\nabla _\alpha T_{\mu \beta }\nabla _\alpha T_\nu ^{~\beta }-\nabla ^\beta T^\alpha _{~(\mu }\nabla _{\nu )}T_{\alpha \beta }\Big ). \end{aligned}$$
(11)

As can be seen from the above relations, the effective energy–momentum tensor and the effective cosmological constant, have contributions from the non-minimal interactions between the Lagrange multiplier \(A_{\mu }\) and the matter field, and also from the non-minimal interactions of the matter field and geometry. Also, it should be noted that the effective energy–momentum tensor includes derivative couplings of the matter sources. This is a generic property of theories with non-minimal matter-geometry couplings [76,77,78,79,80,81,82,83] and also theories with derivative matter couplings [89, 90].

3 Cosmological consequences

Let us assume that the Universe is described by the flat FRW metric

$$\begin{aligned} ds^{2}=a^{2}(-dt^{2}+d\vec {x}^2), \end{aligned}$$
(12)

where \(a=a(t)\) is the scale factor and t is the conformal time. We define the Hubble parameter as \(H={\dot{a}}/a\), describing the rate of expansion of the Universe. Here dot denotes derivative with respect to t.

Let us also assume that the Universe is filled with a perfect fluid, with Lagrangian density \(\mathcal {L}_{m}=-\rho \) and the energy–momentum tensor

$$\begin{aligned} T_{\nu }^{\mu }=\mathrm {diag}\left( -\rho ,p,p,p\right) , \end{aligned}$$
(13)

where \(\rho \) is the energy density of the baryonic matter and p is its thermodynamics pressure. With the above assumptions, the energy–momentum conservation equation (8) takes the form

$$\begin{aligned} {\dot{\rho }}+3H(\rho +p)=0. \end{aligned}$$
(14)

For the Lagrange multiplier, we choose \(A_{\mu }=(A_0(t),\vec {0})\). It should be noted that the Lagrange multiplier is a vector field over FRW space-time. As a result the isotropy and homogeneity of the space-time implies that the vector field \(A_{\mu }\) has only a non-zero temporal component in Cartesian coordinates.

One can then obtain the Friedmann and Raychaudhuri equations as

$$\begin{aligned} 6\kappa ^2 H^2=a^2(\rho +2\kappa ^2\Lambda )+6\alpha H(\rho +p){{\dot{\rho }}}, \end{aligned}$$
(15)

and

$$\begin{aligned}&4\kappa ^2({\dot{H}}-H^2)=-a^2(\rho +p)+A_0\big ({\dot{p}}-H(\rho +p)\big )\nonumber \\&\quad -8Hp{\dot{p}}+2({\dot{H}}+2H^2)(p^2-\rho ^2)\nonumber \\&\quad -4({\dot{H}}-H^2)p\rho +\frac{d^2}{dt^2}(p^2+\rho ^2). \end{aligned}$$
(16)

The covariant divergence of the metric field equation (7) gives

$$\begin{aligned} \left( \frac{3}{2}HA_0-3\alpha p(2H^2+{\dot{H}})\right) \Bigg ({{\dot{\rho }}}+{\dot{p}}+2H(\rho +p)\Bigg )=0. \end{aligned}$$
(17)
Fig. 1
figure 1

The Hubble and deceleration parameters as a function of the redshift z for \(\beta \times 10^{5}=(0.1,0.3,0.6,1.1)\) indicated as dotted, long-dashed, dashed and dot-dashed respectively. The \(\Lambda \)CDM curve is depicted as a red curve

Full size image

It can be easily verified that the second parenthesis would not be zero for matter fields with \(p>-\rho /3\). As a result, one can obtain \(A_0\) from the first parenthesis as

$$\begin{aligned} A_0=2\alpha p\left( 2H+\frac{{\dot{H}}}{H}\right) . \end{aligned}$$
(18)

It is interesting to note that if the Universe is filled only by non-relativistic matter with \(p=0\), the Lagrange multiplier vanishes on top of FRW Universe.

In the following, we will assume that the Universe is filled by radiation with equation of state \(p_r=\rho _r/3\) and dust with equation of state \(p_m=0\). As a result, one has \(\rho =\rho _m+\rho _r\) and \(p=\rho _r/3\), where r/m stand for radiation/dust components, respectively.

Now, define the following set of dimensionless quantities

$$\begin{aligned} \tau&=H_0t,\quad H=H_0h,\quad {\bar{A}}_0=H_0A_0,\nonumber \\ {\bar{\rho }}_i&=\frac{\rho _i}{6\kappa ^2H_0^2},\quad \Omega _{\Lambda 0}=\frac{\Lambda }{3H_0^2},\quad \beta =\alpha \kappa ^2H_0^4, \end{aligned}$$
(19)

where \(i=r,m\). From (14), one can write conservation equations for dust and radiation separately. These equations can then be solved to obtain

$$\begin{aligned} {\bar{\rho }}_r=\frac{\Omega _{r0}}{a^4},\qquad {\bar{\rho }}_m=\frac{\Omega _{m0}}{a^3}, \end{aligned}$$
(20)

where \(\Omega _{r0}=0.53\times 10^{-4}\) and \(\Omega _{m0}=0.305\) are present time density parameters for radiation and dust, respectively [91].

Using dimensionless parameters (19) and also equations (18) and (20), one can obtain the Hubble parameter from Friedmann and Raychaudhuri equations as

$$\begin{aligned} h(z)^2=\frac{(1+z)\Omega _{m0}+(1+z)^2\Omega _{r0}+\frac{\Omega _{\Lambda 0}}{(1+z)^2}}{1+4\beta (1+z)^6f(z)}, \end{aligned}$$
(21)

where we have defined

$$\begin{aligned} f(z)=18\Omega _{m0}^2+45(1+z)\Omega _{m0}\Omega _{r0}+31(1+z)^2\Omega _{r0}^2, \end{aligned}$$
(22)

and we have used the redshift variable defined as

$$\begin{aligned} 1+z=\frac{1}{a}. \end{aligned}$$
(23)

In Fig. 1 we have plotted the Hubble and deceleration parameters, defined as

$$\begin{aligned} q=(1+z)\frac{d \ln h(z)}{dz}, \end{aligned}$$

for four different values of the model parameter \(\beta \times 10^{5}=(0.1,0.3,0.6,1.1)\). It should be noted that in Sect. 5 we will obtain the best estimation of the parameter \(\beta \) using observational data. The best fit value of the parameter is \(\beta =0.6\times 10^{-5}\) and up to \(2\sigma \) confidence level \(0.1<\beta <1.11\). As a result the abovee values are chosen to be in the \(2\sigma \) confidence interval. We have also plotted the \(\Lambda \)CDM curve as a red solid curve in these figures. This is obtained by taking \(\beta =0\) in the field equations. The observational data of the Hubble parameter together with their errors is also shown [92]. It can be seen from the figures that the non-minimal coupling between matter and geometry affects the cosmological behavior of the Universe at redshifts greater than \(z\sim 2\). For the values of the parameter \(\beta \) in \(2\sigma \) confidence interval, the Universe will decelerate less than standard \(\Lambda \)CDM theory. Larger values for \(\beta \) will make the Universe to accelerate at these redshifts, representing also itself in a smaller values of the Hubble parameter. For negative values of \(\beta \), the Universe will have more deceleration at early times, implying that the Universe was larger than the \(\Lambda \)CDM prediction at that times. We have not plotted this case in Fig. 1 since these values are not in the \(2\sigma \) confidence interval. In summary one could see that the late time observational data would be fulfilled with this non-minimal matter/geometry coupling model. So, more observational evidence would be needed to decide which model will fit the Universe more.

Fig. 2
figure 2

The behavior of the difference between effective cosmological constant with its \(\Lambda \)CDM value as a function of the redshift z for \(\beta \times 10^{5}=(0.1,0.3,0.6,1.1)\) indicated as dotted, long-dashed, dashed and dot-dashed respectively

Full size image

In Fig. 2 we have plotted the relative difference between the effective cosmological constant (10) with its \(\Lambda \)CDM value \(\Omega _\Lambda =0.694\) [91]. One can see from the figure that the effective cosmological constant becomes equal to the \(\Lambda \)CDM value at redshifts smaller than \(z\sim 1\), indicating that the theory becomes identical to the \(\Lambda \)CDM model. However, the value of the effective cosmological constant differs from the \(\Lambda \)CDM value for redshifts larger than unity. For values of the parameter \(\beta \) in the \(2\sigma \) confidence interval we have \(\Lambda _{eff}>\Lambda \), implying that the Universe has more acceleration at redshifts \(z>1\) with respect to the \(\Lambda \)CDM model as was discussed before. Also, negative values of the parameter \(\beta \) makes the effective cosmological constant to become smaller than its \(\Lambda \)CDM value and the Universe experiences more deceleration at early times.

Fig. 3
figure 3

The behavior of the Lagrange multiplier \({\bar{A}}_0\) as a function of the redshift z for \(\beta \times 10^{5}=(0.1,0.3,0.6,1.1)\) indicated as dotted, long-dashed, dashed and dot-dashed respectively

Full size image

In Fig. 3, we have plotted the evolution of the Lagrange multiplier \({\bar{A}}_0\) as a function of the redshift z using equation (18). One can see from the figure that for the values of parameter \(\beta \) in the \(2\sigma \) confidence range, the Lagrange multiplier is positive. It should also be noted that for small \(\beta \) values, the Lagrange multiplier tends to zero for \(z\rightarrow 0\). However, for larger values of \(\beta \), the Lagrange multiplier remains non-vanishing. This comes back to the fact that as we take the smaller values of the \(\beta \) parameter, the theory tends to the standard \(\Lambda \)CDM model which is conservative. As a result there is no need for an extra force to keep the theory conservative. By letting \(\beta \) to adopt larger values, the theory differs significantly from the \(\Lambda \)CDM theory and one should non-trivially keep it conservative through a Lagrange multiplier. At redshifts larger than unity, as we have discussed earlier, the theory deviates from the \(\Lambda \)CDM theory and one needs a non-zero Lagrange multiplier for every values of \(\beta \).

Fig. 4
figure 4

The phase space portrait of the dynamical system (26). We have plotted the \((\Omega _m,\Omega _r)\) planes for three different values of \(\Omega _\alpha =-1\) (left), \(\Omega _\alpha =0\) (center) and \(\Omega _\alpha =1\) (right). The corresponding fixed points are also shown

Full size image

4 Dynamical system analysis

In this section, we will consider the dynamical analysis of the system of Eqs. (15) and (16). Let us introduce the following set of dynamical variables

$$\begin{aligned} \Omega _m&=\frac{\rho _m a^2}{6\kappa ^2H^2},\qquad \Omega _r=\frac{\rho _r a^2}{6\kappa ^2H^2},\nonumber \\ \Omega _\Lambda&=\frac{\Lambda a^2}{3H^2},\qquad \quad \Omega _\alpha =\frac{\alpha \kappa ^2H^4}{a^4}. \end{aligned}$$
(24)

It should be noted that the the first three variables are the standard matter density parameters for dust, radiation and the cosmological constant respectively. The last one is related to the non-minimal coupling between matter and geometry which is related to the parameter \(\alpha \). Substituting the above variables in the Friedmann equation (15), one can obtain \(\Omega _\Lambda \) as

$$\begin{aligned} \Omega _\Lambda =1-\Omega _m-\Omega _r+\Omega _\alpha (72\,\Omega _m^2+127\,\Omega _r^2+180\,\Omega _m\Omega _r), \end{aligned}$$
(25)

where we have used the conservation equations for dust and radiation (20). It is then evident that the system has three dynamical variables. Using the Raychaudhuri equation (16), one can obtain the following set of autonomous 3D dynamical system

$$\begin{aligned} \frac{d\Omega _m}{dN}&=-A\Omega _m\Big (3+648\,\Omega _\alpha \Omega _m^2+4\,\Omega _r(341\,\Omega _\alpha \Omega _r-1)\nonumber \\ {}&\quad +3\,\Omega _m(660\,\Omega _\alpha \Omega _r-1)\Big ), \end{aligned}$$
(26a)
$$\begin{aligned} \frac{d\Omega _r}{dN}&=-A\Omega _r\Big (4+720\,\Omega _\alpha \Omega _m^2+4\,\Omega _r(372\,\Omega _\alpha \Omega _r-1)\nonumber \\ {}&\quad +3\,\Omega _m(660\,\Omega _\alpha \Omega _r-1)\Big ), \end{aligned}$$
(26b)
$$\begin{aligned} \frac{d\Omega _\alpha }{dN}&=2A\Omega _\alpha \Big (432\,\Omega _\alpha \Omega _m^2+4\,\Omega _r(248\,\Omega _\alpha \Omega _r-1) \nonumber \\ {}&\quad +3\,\Omega _m(420\,\Omega _\alpha \Omega _r-1)\Big ), \end{aligned}$$
(26c)

where \(N=\ln a\) and we have defined

$$\begin{aligned} A=1+4\Omega _\alpha (18\,\Omega _m^2+45\,\Omega _m\Omega _r+31\,\Omega _r^2). \end{aligned}$$
(27)

The effective equation of state parameter \(\omega _{eff}\) can also be calculated as

$$\begin{aligned} \omega _{eff}&\equiv -\frac{1}{3}-\frac{2}{3}\frac{{\dot{H}}}{H^2}\nonumber \\ {}&=-\frac{1}{3}A\Big (3+648\,\Omega _\alpha \Omega _m^2+4\,\Omega _r(341\,\Omega _\alpha \Omega _r-1)\nonumber \\ {}&\quad ~+3\,\Omega _m(600\,\Omega _\alpha \Omega _r-1)\Big ), \end{aligned}$$
(28)

The above dynamical system, has three different fixed points which we will discuss in the following.

4.1 Matter dominated fixed point

The first fixed point of the system (26) is

$$\begin{aligned} (\Omega _m,\Omega _r,\Omega _\alpha )=(1,0,0), \end{aligned}$$

with the effective equation of state parameter \(\omega _{eff}=0\). As a result, in this point we have a dust dominated Universe. The eigenvalues of this fixed point are \((-6,3,-1)\), indicating that the dust fixed point is a saddle point.

4.2 Radiation dominated fixed point

We also have a fixed point of the system (26) corresponding to the radiation dominated phase

$$\begin{aligned} (\Omega _m,\Omega _r,\Omega _\alpha )=(0,1,0), \end{aligned}$$

since the equation of state parameter is \(\omega _{eff}=1/3\). The eigenvalues of this fixed point are \((-8,4,1)\), indicating that the radiation fixed point is a saddle point.

4.3 de Sitter fixed line

The last fixed point of the system (26) is

$$\begin{aligned} (\Omega _m,\Omega _r,\Omega _\alpha )=(0,0,\Omega _\alpha ), \end{aligned}$$

which is true for every value of the variable \(\Omega _\alpha \). As a result, we have a fixed line with the equation of state parameter \(\omega _{eff}=-1\). This corresponds to the de Sitter expansion of the Universe. The eigenvalues of this fixed line is \((-4,-3,0)\), indicating that the de Sitter fixed line is stable.

In Fig. 4, we have plotted the \((\Omega _m,\Omega _r)\) phase space planes for three different values \(\Omega _\alpha =-1,0,1\). We have also shown the fixed points in the figures. It should be noted that the dust and radiation fixed points are located at \(\Omega _\alpha =0\) plane, which can be seen in the second plot. The de Sitter line is the line perpendicular to the planes and crosses the (0, 0) point. As one can see from the figures, the de Sitter fixed line is an attractor and all lines will end at this fixed point eventually. It is interesting to note that if one starts from the radiation dominated fixed point there is a flow which let the Universe to transform to dust dominated phase and then continue its evolution to the de Sitter stable fixed point. The theory can then in principle has a history such that the radiation and matter dominated phases occurs before the Universe falls into the late time de Sitter evolution. As a result the thermal history of the Universe can be recovered in this model.

5 Matter density perturbations

Let us consider the first order scalar perturbations of the model (6). The perturbed metric can be written in the Newtonian gauge as

$$\begin{aligned} ds^2=a^2(t)\Big [-(1+2\varphi )dt^2+(1-2\psi )d\vec {x}^2\Big ], \end{aligned}$$
(29)

where \(\varphi \) and \(\psi \) are the Bardeen potentials. The perturbed energy momentum tensor can also be written as

$$\begin{aligned}&\delta T^0_0=-\delta \rho \equiv -\rho \delta ,\nonumber \\ {}&\delta T^0_i=(1+w)\rho \partial _i v,\quad \nonumber \\ {}&\delta T^i_j=\delta ^i_jc_s^2\rho \delta . \end{aligned}$$
(30)

In the above expression, we have defined the matter density contrast \(\delta \) as \(\delta =\delta \rho /\rho \). Also v is the scalar mode of the velocity perturbation associated with the matter sector, \(c_s\) is the sound speed and w is the equation of state parameter of baryonic matter \(p=w\rho \). In the following, we will assume that the Universe at the time where we are performing the perturbations is in the dust dominated phase. As a result, the perturbed and unperturbed matter content of the Universe satisfy \(\delta p/\delta \rho =c_s^2=p/\rho \), where \(c_s^2=0\).

For the vector field \(A_\mu \) we define the first order scalar perturbation as

$$\begin{aligned} \delta A_\mu =a(\mathcal {A}_0,-\partial _i\mathcal {A}). \end{aligned}$$
(31)

Perturbing the conservation Eq. (8) to first order in perturbation variables gives

$$\begin{aligned} {{\dot{\delta }}}+\left( 3H+\frac{{{\dot{\rho }}}}{\rho }\right) (\delta -2\varphi )-3{{\dot{\psi }}}+\theta =0, \end{aligned}$$
(32)

and

$$\begin{aligned} {{\dot{\theta }}}+H\theta -k^2\varphi =0, \end{aligned}$$
(33)

where we have defined \(\theta =\nabla _i\nabla ^i v\), and we transform to the Fourier coordinates, where \(\vec {k}\) is the wave vector. Using the above equations and also the background field equation (14), one obtains a differential equation of the evolution of the matter density contrast as

$$\begin{aligned} \ddot{\delta }+H{{\dot{\delta }}}+k^2\varphi -3H{{\dot{\psi }}}-3\ddot{\psi }=0. \end{aligned}$$
(34)

Let us consider the sub-horizon limit of the theory where \(k\gg H\). In this limit the (00) component of the metric field equation and the (0) component of vector field equation at first order in perturbations reduces to

$$\begin{aligned}&2\alpha H\rho ^2{{\dot{\delta }}}+(a^2\rho +4\alpha k^2\rho ^2-24\alpha H^2\rho ^2)\delta \nonumber \\&\quad +k^2a\rho \mathcal {A}+4\kappa ^2 k^2\psi =0, \end{aligned}$$
(35)
$$\begin{aligned}&a\mathcal {A}+\alpha \rho (\varphi +3\psi )=0. \end{aligned}$$
(36)

Also, the \(i\ne j\) components of the metric field equation (6) reads

$$\begin{aligned} \kappa ^2(\varphi -\psi )=a\rho \mathcal {A}. \end{aligned}$$
(37)

From Eqs. (35) to (37). one can obtain the perturbation variables \(\varphi \), \(\psi \) and \(\mathcal {A}\) as

$$\begin{aligned} \psi =\frac{\kappa ^2+\alpha \rho ^2}{\kappa ^2-3\alpha \rho ^2}\varphi ,\quad \mathcal {A}=-\frac{4\alpha \kappa ^2\rho }{a(\kappa ^2-3\alpha \rho ^2)}\varphi . \end{aligned}$$
(38)

and

$$\begin{aligned} k^2\varphi =&-\frac{\rho (\kappa ^2-3\alpha \rho ^2)}{4\kappa ^4}\Big [2\alpha H\rho {{\dot{\delta }}}\nonumber \\&+(a^2+4\alpha k2\rho -24\alpha H^2\rho )\delta \Big ]. \end{aligned}$$
(39)

Using Eqs. (39) in (34), one obtains for the evolution of the matter density contrast

$$\begin{aligned}&\ddot{\delta }+\Bigg [1-\frac{\alpha }{2\kappa ^2}(\kappa ^2-3\alpha \rho ^2)\rho ^2\Bigg ]H{{\dot{\delta }}}\nonumber \\&\quad -\frac{\kappa ^2-3\alpha \rho ^2}{4\kappa ^4}\Big [a^2+4\alpha (\kappa ^2-6H^2)\rho \Big ]\rho \delta =0. \end{aligned}$$
(40)

It should be noted that in the case \(\alpha =0\), the above equation reduces to the standard equation in \(\Lambda \)CDM theory. One can see from the above equation that the presence of the non-minimal matter/geometry coupling affects the effective gravitational constant and also the friction term in the evolution equation of the density contrast.

Transforming to dimensionless variables (19), noting that the quantity \(\delta \) is dimensionless by its own and using the expression (20) for the dust energy density \(\rho _m\), one obtains

$$\begin{aligned}&\delta ^{\prime \prime }+\left[ \frac{h^\prime }{h}+18\beta \,\Omega _{m0}^2(1+z)^5\left( 1-108\beta \,\Omega _{m0}^2(1+z)^6\right) \right] \delta ^\prime \nonumber \\&\quad -\frac{3}{2}\frac{\Omega _{m0}}{(1+z)h^2}\Big (1-108\beta \,\Omega _{m0}^2(1+z)^6\Big )\nonumber \\&\quad \times \Big (1+24\beta \,\Omega _{m0}(1+z)^5(\gamma ^2-6h^2)\Big )\delta =0, \end{aligned}$$
(41)

where we have transformed to the redshift coordinates and prime denotes derivative with respect to the redshift z. It should be noted that the evolution equation of the matter density contrast depends on \(\gamma =k/H_0\).

In order to solve the above equation, we will use the same initial conditions as in \(\Lambda \)CDM theory in deep matter dominated era in which

$$\begin{aligned} \frac{d\delta }{d\ln a}|_{z_\star }=\delta |_{z_\star }, \end{aligned}$$
(42)

where \(z_\star \) is some point in the deep matter dominated era which we will assume to be \(z_\star =7.1\).

In order to compare the model with observational data, we will use the data set on the Hubble parameter in the redshift range \(0<z<2\) [92] and also the observational data on \(f\sigma _8\) [93]

$$\begin{aligned} f \sigma _8 \equiv \sigma _8(z)f(z), \end{aligned}$$
(43)

where \(\sigma _8(z)=\sigma _8^0\, \delta (z)/\delta (0)\), and \(\sigma _8^0\) is a model dependent constant. We have defined the growth rate of matter perturbations as

$$\begin{aligned} f=\frac{d \ln \delta }{d\ln a}=-(1+z)\frac{\delta '}{\delta }. \end{aligned}$$
(44)

We estimate the values of \(\sigma _8^0\) and also the model parameter \(\beta \) by maximizing the Likelihood function

$$\begin{aligned} {\mathcal {L}}={\mathcal {L}}_{0}e^{-\chi ^2/2}, \end{aligned}$$
(45)

where \({\mathcal {L}}_0\) is the normalization constant and the quantity \(\chi ^2\) in our case is given by

$$\begin{aligned} \chi ^2&=\chi _H^2+\chi _{f\sigma _8}^2\nonumber \\&=\sum _i\left( \frac{H_{i,obser}-H_{i,theory}}{\sigma _i}\right) ^2\nonumber \\&\quad +\sum _j\left( \frac{f\sigma _{8j,obser}-f\sigma _{8j,theory}}{\sigma _j}\right) ^2, \end{aligned}$$
(46)

where \(H_{i,theory}\) and \(f\sigma _{8j,theory}\) are the theoretical values for the observables \(H_{i,obser}\) and \(f\sigma _{8j,obser}\) and \(\sigma _i\) is the error of the ith data. It should be noted that the two sets of data on H and \(f\sigma _8\) that we used in this paper are independent, so the total Likelihood function is obtained by multiplying the individual Likelihoods of the two sets.

Table 1 The best fit values together with their \(1\sigma \) and \(2\sigma \) confidence intervals of the model parameter \(\beta \) and also for \(\sigma _8^0\)
Full size table

In Table 1 we have summarized the best fit values of the parameters \(\beta \) and \(\sigma _8^0\) and their \(1\sigma \) and \(2\sigma \) confidence intervals. It should be noted that for \(\sigma _8^0\) at the best fit point, we have \(\chi ^2/dof=0.44\).

In Fig. 5, we have plotted the evolution of the \(f\sigma _8\) as a function of redshift z for the best fit values shown in Table 1 and \(\gamma =1.4\).

Fig. 5
figure 5

The evolution of the \(f\sigma _8\) as a function of redshift z. The Solid curves correspond to the \(\Lambda \)CDM theory and the dashed line is associated with our model with \(\beta =0.6\times 10^{-5}\) and \(\gamma =1.4\). We have also plotted the observational data with their errors in the figure

Full size image

One can see from Fig. 5 that the qualitative behavior of the \(f\sigma _8\) function is similar to the \(\Lambda \)CDM model. For small redshift values, one can see that the value of the \(f\sigma _8\) is less than that of \(\Lambda \)CDM theory. This implies that this model predicts slower growth rate of matter fields. However, as an observational side, both theories are satisfactory and so, more observations would be needed to favor one of them.

6 Conclusions and final remarks

In this paper, we have considered a modified theory of gravity containing non-minimal coupling between matter and geometry. The new term can be considered as a generalization of the energy–momentum squared gravity by coupling it to the Ricci tensor. As, we have a non-minimal coupling between matter fields and gravity, the energy–momentum tensor is no longer conserved. We have fixed this problem by adding a constraint through Lagrange multiplier to ensure that the equation \(\nabla _\mu T^{\mu \nu }=0\) holds. The Lagrange multiplier in fact is a vector field and can be considered as a potential to keep the matter components to behave conservative. In the FRW background we have plotted the Lagrange multiplier in Fig. 3. In this figure, it can be seen that the Lagrange multiplier is always positive, which means that we should do positive work to keep the theory conservative. Also, the Lagrange multiplier tends to zero at the present time. Closer we are to the early times, stronger force we should imply to make the theory conservative. Moreover, the covariant divergence of the metric field equation becomes the dynamical equation for the Lagrange multiplier. In the FRW Universe, this equation becomes an algebraic equation for the Lagrange multiplier and one can obtain the Lagrange multiplier as a function of pressure and the Hubble parameter (18). This confirms that the source for the conservative force is in fact the non-minimal nature of the matter/geometry coupling itself. In the case of dust dominated Universe, the Lagrange multiplier vanishes. This is true only at the background level since at the first order perturbation level, the Lagrange multiplier is always non-zero unless the energy density vanishes.

In this paper, we have found the best fit values for the model paramter \(\alpha \) and also for \(\sigma _8^0\) by using two sets of independent data corresponding to the Hubble parameter and also the \(f\sigma _8\) function. In order to satisfy observational data, the best fit value of the model paramter \(\alpha \) and also its 2\(\sigma \) confidence interval shows that this parameter should be small and positive. In this range, the behavior of the Hubble parameter coincides with the \(\Lambda \)CDM theory for redshifts smaller than \(z\sim 2\). However for larger redshifts the Hubble parameter becomes smaller. This shows that the present theory predicts larger Hubble radius at that times compared to the \(\Lambda \)CDM model. The analysis of the deceleration parameter shows that the Universe at early times could be in the accelerating phase for large enough values of the parameter \(\alpha \). For the best fit value of \(\alpha \), however, the behavior of the Universe at early times is qualitatively the same as \(\Lambda \)CDM theory, despite the fact that we have less deceleration in the present theory.

Fig. 6
figure 6

The behavior of the effective equation of state parameter as a function of the redshift z for \(\beta \times 10^{5}=(0.1,0.3,0.6,1.1)\) indicated as dotted, long-dashed, dashed and dot-dashed respectively

Full size image

It is well-known that in the semi-classical regime, the vacuum expectation value of the energy–momentum tensor renormalizes the cosmological constant \(\Lambda \), as well as the gravitational constant G. In this paper, we have considered a non-minimal coupling between matter and geometry through the interaction of the Ricci tensor with the energy–momentum tensor squared. In the semi-classical regime this new coupling would contribute to the renormalization of the coefficients of the Ricci tensor and also to derivatives of the Riemann tensor. Remembering the definition of the effective cosmological constant (10), one expects that quantum corrections to the matter fields, renormalizes \(\Lambda _{eff}\) as well. It should also be noted that the effective energy–momentum tensor does not alter the exact moment of phase transitions at early times. This is because the theory is constrained to have a conservation of the energy–momentum tensor as we have implied by introducing a Lagrange multiplier to the action. However, the behavior of the matter abundances in radiation, dust and dark energy dominated phases differ from the standard \(\Lambda \)CDM model.

We have also analyzed the dynamical evolution of the model. The theory is a three dimensional autonomous system. We have shown that the theory has three fixed points similar to \(\Lambda \)CDM theory, corresponding to the dust, radiation and de Sitter expansion. However, in the present model, the fixed point corresponding to the de Sitter expansion is changed to a fixed line which is the \(\Omega _\alpha \) axes in the phase diagram of the theory. However, in the \(\Omega _\alpha =0\) plane the qualitative behavior of the theory is the same as in the \(\Lambda \)CDM theory.

We have also considered the first order perturbation analysis of the theory. Since we have a non-minimal matter-geometry coupling, the anisotrpic stress \(\eta =(\varphi -\psi )/\varphi \) is non-vanishing. From Eq. (41), one can see that the non-minimal coupling adds some new terms proportional to at least \((1+z)^5\) to the evolution equation of the density contrast. This shows that the matter/geometry coupling would be important for large values of z. For small redshifts, these terms are negligible and the qualitative behavior of the density contrast is similar to the \(\Lambda \)CDM theory. This can also be seen from the plot of \(f\sigma _8\) in Fig. 5 where one can see that the late time observational data could be satisfied in this theory. It should be noted that the value of \(f\sigma _8\) is smaller than the corresponding value of \(\Lambda \)CDM mode.

Let us note about the \(H_0\) and \(\sigma _8\) tensions in this model. In general, dynamical dark energy theories with \(\omega <-1\) could in principle loosen the \(H_0\) tension but they worsen the \(\sigma _8\) tension. Conversely, theories with \(\omega >-1\) does not change the status of \(H_0\) tension but they loosen the \(\sigma _8\) tension. In Fig. 6, we have plotted the \(\omega _{eff}\) as a function of the red-shift z for \(\beta \times 10^{5}=(0.1,0.3,0.6,1.1)\). One can see from the figure that our model belongs to the sub-class of dynamical dark-energy models with \(\omega >-1\). As a result we expect that the theory would make the \(H_0\) tension better but consequently will make the \(\sigma _8\) tension worse. However, more detailed analysis would be needed to investigate this point.

At last, we would like to note that the non-minimal coupling between matter and geometry could in principle explains late time cosmology, but more analysis is needed to choose the best theory among many late time cosmological models.