Generalized Curvature-Matter Couplings in Modified Gravity

Abstract

In this work, we review a plethora of modified theories of gravity with generalized curvature-matter couplings. The explicit nonminimal couplings, for instance, between an arbitrary function of the scalar curvature R and the Lagrangian density of matter, induces a non-vanishing covariant derivative of the energy-momentum tensor, implying non-geodesic motion and, consequently, leads to the appearance of an extra force. Applied to the cosmological context, these curvature-matter couplings lead to interesting phenomenology, where one can obtain a unified description of the cosmological epochs. We also consider the possibility that the behavior of the galactic flat rotation curves can be explained in the framework of the curvature-matter coupling models, where the extra terms in the gravitational field equations modify the equations of motion of test particles and induce a supplementary gravitational interaction. In addition to this, these models are extremely useful for describing dark energy-dark matter interactions and for explaining the late-time cosmic acceleration.

1. Introduction

Recent observations of supernovae, together with the Wilkinson Microwave Anisotropy Probe (WMAP) and Sloan Digital Sky Survey (SDSS) data, lead to the remarkable conclusion that our universe is not just expanding, but has begun to accelerate [1,2,3]. The resolution of this fundamental question is extremely important for theoretical cosmology, looking beyond the standard theory of gravity. The standard model of cosmology has favored the dark energy models as fundamental candidates responsible for the cosmic expansion. However, it is clear that these questions involve not only gravity, but also particle physics. String theory provides a synthesis of these two branches of physics and is widely believed to be moving towards a viable quantum gravity theory. One of the key predictions of string theory is the existence of extra spatial dimensions. In the braneworld scenario, motivated by recent developments in string theory, the observed three-dimensional universe is embedded in a higher-dimensional spacetime [4]. The new degrees of freedom belong to the gravitational sector and can be responsible for the late-time cosmic acceleration.

Indeed, detailed theoretical and phenomenological analysis of the relation between the effective field yielding dark energy, non-canonical Lagrangians and non-linear gauge kinetic functions have been extensively explored in the literature. Note that generalizations of the action functional can be approached in several ways. For instance, prescriptions consist in replacing the linear scalar curvature term in the Einstein–Hilbert action by a function of the scalar curvature, $f\left(R\right)$, or by more general scalar invariants of the theory. This class of theories is often termed higher-order gravity theories [5,6,7]. In this context, infra-red modifications of General Relativity have been extensively explored, including four-dimensional modifications to the Einstein–Hilbert action, and the consistency of various candidate models has been extensively analyzed. All modified gravity theories induce observational signatures at the post-Newtonian level, which are translated by the parameterized post-Newtonian (PPN) metric coefficients arising from these extensions of General Relativity (GR). Thus, generalizations of the Einstein–Hilbert Lagrangian, including quadratic Lagrangians that involve second order curvature invariants have also been extensively explored [5]. While the latter modified theories of gravity offer an alternative explanation to the standard cosmological model [8], it offers a paradigm for nature fundamentally distinct from dark energy models of cosmic acceleration, even those that perfectly mimic the same expansion history. It is a fundamental question to understand how one may observationally and theoretically differentiate these modified theories of gravity from dark energy models. Thus, one should test these models against large-scale structure and lensing, astrophysical and laboratory measurements, as well as laboratory and space-based equivalence principle experiments. These tests from the Solar System, the large-scale structure and lensing essentially restrict the range of allowed modified gravity models and, thus, offer a window into understanding the perplexing nature of the cosmic acceleration, and of gravity itself.

In this review, we consider an interesting possibility, which includes non-minimal couplings between the scalar curvature and the matter Lagrangian density, introduced in [9], and its extensions and generalizations. The specific cases of a curvature coupling to dark energy [10,11,12], the Maxwell field [13] and an explicit curvature-Yang–Mills coupling [14] was also explored in the context of inflation and of the late-time cosmic acceleration. Indeed, in the context of $f\left(R\right)$ modified theories of gravity, it was shown that an explicit coupling between an arbitrary function of the scalar curvature R and the Lagrangian density of matter induces a non-vanishing covariant derivative of the energy-momentum, implying non-geodesic motion and, consequently, leads to the appearance of an extra force [9]. The latter extra force is orthogonal to the four-velocity, and the corresponding acceleration law was obtained in the weak field limit. Connections with Modified Newtonian Dynamics (MOND) and with the Pioneer anomaly were further discussed. These curvature-matter coupling theories include the more evolved generalizations of $f(R,L_m)$ [15] and $f(R,T)$ gravities [16]. Amongst other features, these models allow for an explicit breaking of the equivalence principle (EP), which is highly constrained by Solar System experimental tests [17,18], by imposing a matter-dependent deviation from geodesic motion.

Note that the weak equivalence principle is considered one of the pillars of General Relativity, and in fact, even a sizable part of the modified-gravity community considers this principle as truly fundamental [19]. This fact has to be stressed, because it demonstrates somehow a limitation for this class of theories. However, it has been recently reported, from data of the Abell Cluster A586, that the interaction of dark matter and dark energy does imply the violation of the equivalence principle [20]. Notice that the violation of the equivalence principle is also found as a low-energy feature of some compactified version of higher-dimensional theories. Indeed, as emphasized by Thibault Damour, it is important to note that the EP is not one of the “universal” principles of physics [21]. It is a heuristic hypothesis, which was introduced by Einstein in 1907, and used by him to construct his theory of General Relativity. In modern language, the (Einsteinian) Equivalence Principle (EP) consists in assuming that the only long-range field with gravitational-strength couplings to matter is a massless spin-2 field. Modern unification theories, and notably String Theory, suggest the existence of new fields (in particular, scalar fields: “dilaton” and “moduli”) with gravitational-strength couplings. In most cases, the couplings of these new fields violate the EP. If the field is long-ranged, these EP violations lead to many observable consequences, such as the variation of fundamental “constants”, the non-universality of free fall and the relative drift of atomic clocks, amongst others. The best experimental probe of a possible violation of the EP is to compare the free-fall acceleration of different materials. Further tests of this principle remain important and relevant for new physics and do indeed strongly restrict the parameters of the considered theory [22,23]. However, it is important to note that, in this context, the EP does not in principle rule out the specific theory.

More specifically, $f(R,L_m)$ gravity further generalizes $f\left(R\right)$ gravity by assuming that the gravitational Lagrangian is given by an arbitrary function of the Ricci scalar R and of the matter Lagrangian $L_m$ [15]. This may also be considered a maximal extension of the Einstein–Hilbert action. Here, we use the term “maximal extension” in a strict mathematical sense, and we define it as “an element in an ordered set that is followed by no other”. For instance, if we consider the set $\left\{\mathrm{Grav}\mathrm{Theor}\right\}$ of gravity theories as given by $\{\mathrm{Grav}\mathrm{Theor}\}=\left\{\mathrm{standard}\mathrm{general}\mathrm{relativity},f\left(R\right)\mathrm{gravity},\mathrm{linear}\mathrm{curvature}-\mathrm{matter}\mathrm{coupling},...,f\left(R,L_m\right)\right\}$, the $f\left(R,L_m\right)$ gravity theory represent the maximal extension of the set of all gravitational theories constructed in a Riemann space and with the action depending on the Ricci scalar and the matter Lagrangian only. In [15], the gravitational field equations of the $f\left(R,L_m\right)$ gravity theory in the metric formalism were obtained, as well as the equations of motion for test particles, which follow from the covariant divergence of the energy-momentum tensor. The equations of motion for test particles can also be derived from a variational principle in the particular case in which the Lagrangian density of the matter is an arbitrary function of the energy-density of the matter only. In [15], the Newtonian limit of the equation of motion was also considered, and a procedure for obtaining the energy-momentum tensor of the matter was presented. In particular, the gravitational field equations and the equations of motion for a particular model in which the action of the gravitational field has an exponential dependence on the standard general relativistic Hilbert–Einstein Lagrange density were also derived. The $f(R,L_m)$ gravitational theory was further generalized by considering the novel inclusion of a scalar field and a kinetic term constructed from the gradients of the scalar field, respectively [24]. Specific models with a nonminimal coupling between the scalar field and the matter Lagrangian were further explored. We emphasize that these models are extremely useful for describing an interaction between dark energy and dark matter and for explaining the late-time cosmic acceleration.

In the context of $f(R,T)$ modified theories of gravity, the gravitational Lagrangian is given by an arbitrary function of the Ricci scalar R and of the trace of the energy-momentum tensor T [16]. The gravitational field equations in the metric formalism were obtained, and it was shown that these field equations explicitly depend on the nature of the matter source. The field equations of several particular models, corresponding to some explicit forms of the function $f(R,T)$, were also presented. Furthermore, in [16], the Newtonian limit of the equation of motion was also analyzed, and constraints on the magnitude of the extra-acceleration were obtained by analyzing the perihelion precession of the planet Mercury in the framework of the present model. An interesting specific case, namely the $f(R,T^{\varphi })$ model, was also analyzed in detail, where $T^{\varphi }$ is the trace of the stress-energy tensor of a self-interacting scalar field. The cosmological implications of the model were briefly considered. We refer the reader to [16] for more details.

The non-minimally curvature-matter-coupled $f(R,T)$ gravitational theory was further generalized by considering the inclusion of a contraction of the Ricci tensor with the matter energy-momentum tensor [25,26,27]. We emphasize that examples of such couplings can be found in the Einstein–Born–Infeld theories [28] when one expands the square root in the Lagrangian. An interesting feature of this theory is that in considering a traceless energy-momentum tensor, i.e., $T=0$, the field equations of $f(R,T)$ gravity reduces to those of $f\left(R\right)$ gravity theories, while the presence of the $R_{\mu \nu }{T}^{\mu \nu }$ coupling term still entails a non-minimal coupling to matter. In [25], the Newtonian limit of the $f(R,T,R_{\mu \nu }{T}^{\mu \nu })$ gravitational theory was considered, and an explicit expression for the extra-acceleration, which depends on the matter density, was obtained in the small velocity limit for dust particles. The so-called Dolgov–Kawasaki instability [29] was also analysed in detail, and the stability conditions of the model with respect to local perturbations was obtained. A particular class of gravitational field equations can be obtained by imposing the conservation of the energy-momentum tensor [25]. In this context, the corresponding field equations for the conservative case was also derived by using a Lagrange multiplier method, from a gravitational action that explicitly contains an independent parameter multiplying the divergence of the energy-momentum tensor [25]. The cosmological implications of the model were extensively investigated for both the conservative and non-conservative cases, and several classes of analytical solutions were obtained. In [26], the Friedmann-Lemaître-Robertson-Walker (FRLW) cosmological dynamics for several versions of the $f(R,T,R_{\mu \nu }{T}^{\mu \nu })$ gravity theory was also considered. The reconstruction of the above action was explicitly analyzed, including the numerical reconstruction for the occurrence of the ΛCDM (Lambda-Cold Dark Matter) universe. De Sitter universe solutions in the presence of non-constant fluids were also presented, and the problem of matter instability was further discussed.

All of the above gravitational modifications are based on the curvature description of gravity. However, an interesting and rich class of modified gravity arises by modifying the action of the equivalent torsional formulation of General Relativity. The latter approach has been denoted the “Teleparallel Equivalent of General Relativity” (TEGR) [30,31,32,33,34,35] and consists of replacing the torsion-less Levi–Civita connection by the curvature-less Weitzenböck one and using the vierbein instead of the metric as the fundamental field. In this formulation, the Lagrangian of the theory is constructed by contractions of the torsion tensor. Thus, in analogy to $f\left(R\right)$ gravity, if one wishes to modify gravity in this formulation, the simplest approach would be in generalizing the torsion scalar $\mathcal{T}$ to an arbitrary function $f\left(\mathcal{T}\right)$ [36,37,38].

Now, in the context of TEGR models, one may also construct an extension of $f\left(\mathcal{T}\right)$ gravity with the inclusion of a non-minimal torsion-matter coupling in the action [39]. The resulting theory is a novel gravitational modification, since it is different from both $f\left(\mathcal{T}\right)$ gravity, as well as from the non-minimal curvature-matter-coupled theory. The cosmological application of this new theory proves to be very interesting. In particular, an effective dark energy sector was obtained, where the equation-of-state parameter can be quintessence or phantom-like, or cross the phantom-divide, while for a large range of the model parameters, the Universe results in a de Sitter, dark-energy-dominated, accelerating phase. Additionally, early-time inflationary solutions were also obtained, and thus, one can provide a unified description of the cosmological history. We refer the reader to [39] for more details.

In addition to the latter, non-minimal torsion-matter coupling, an alternative extension of $f\left(\mathcal{T}\right)$ gravity was explored by considering a general coupling of the torsion scalar $\mathcal{T}$ with the trace of the matter energy-momentum tensor T [40]. The resulting $f(\mathcal{T},T)$ theory is a new modified gravity, since it is different from all of the existing torsion- or curvature-based constructions. Applied to a cosmological framework, it also leads to interesting phenomenology. In particular, one can obtain a unified description of the initial inflationary phase, the subsequent non-accelerating, matter-dominated expansion and, then, the transition to late-time accelerating phase. In the far future, the Universe results either in a de Sitter exponential expansion or in eternal power-law accelerated expansions. A similar analysis was investigated in [41], where using a perturbational approach, the stability of the solutions and, in particular, of the de Sitter phase was explored. Furthermore, the constraints imposed by the energy conditions were also considered. We refer the reader to [40,41] for more details.

The possibility that the behavior of the rotational velocities of test particles gravitating around galaxies can be explained in the framework of modified gravity models with nonminimal curvature-matter coupling was considered in [42]. Generally, the dynamics of test particles around galaxies, as well as the corresponding mass deficit, is explained by postulating the existence of dark matter. The extra terms in the gravitational field equations with curvature-matter coupling modify the equations of motion of test particles and induce a supplementary gravitational interaction. Starting from the variational principle describing the particle motion in the presence of the non-minimal coupling, the expression of the tangential velocity of a test particle, moving in the vacuum on a stable circular orbit in a spherically symmetric geometry, was derived. The tangential velocity depends on the metric tensor components, as well as on the coupling function between matter and geometry. The Doppler velocity shifts were also obtained in terms of the coupling function. If the tangential velocity profile is known, the coupling term between matter and curvature can be obtained explicitly in an analytical form. The functional form of this function was obtained in two cases, for a constant tangential velocity and for an empirical velocity profile obtained from astronomical observations, respectively. Therefore, these results open the possibility of directly testing the modified gravity models with a non-minimal curvature-matter coupling by using direct astronomical and astrophysical observations at the galactic or extra-galactic scale. These issues will be reviewed in detail below.

It is rather important to discuss the theoretical motivations for these curvature-matter coupling theories more carefully. For instance, one may wonder if there is a fundamental theory or model that reduces to one of these theories in some limit. Indeed, for the specific case of the linear nonminimal curvature-matter coupling, one can show that this theory can be expressed as a scalar-tensor theory with two scalar fields [43]. These are reminiscent of the extensions of scalar-tensor gravity, which include similar couplings, such as the theories considered by Damour and Polyakov [44], where it was shown that string-loop modifications of the low-energy matter couplings of the dilaton may provide a mechanism for fixing the vacuum expectation value of a massless dilaton. In addition to this, the results presented in [44] provide a new motivation for trying to improve by several orders of magnitude the various experimental tests of Einstein’s equivalence principle, such as the universality of free fall and the constancy of the fundamental constants, amongst others. In addition to this, note that the modern revival of the Kaluza–Klein theory often leads to the introduction of several scalar fields, i.e., the “compacton” [45]. Superstring theory also leads to several scalar fields coupled to the macroscopic distribution of energy; indeed, the “dilaton” is already present in ten dimensions, and the “compacton” comes from a dimensional reduction [45]. Quantum-motivated, higher order generalizations of Einstein’s General Relativity, under some conditions, are also equivalent to adding several scalar fields to the Einstein–Hilbert action [46]. In addition to this, we note that the nonminimal curvature-matter coupling also arises from one-loop vacuum-polarization effects in the formulation of quantum electrodynamics in curved spacetimes [47]. All of the above considerations make it natural to consider modified fields of gravity containing several scalar fields, which, in some cases, are equivalent, in the present context, to the presence of a nonminimal curvature-matter couplings. These considerations further motivate the analysis of the coupling between matter and curvature.

In this work, based on published literature, we review these generalized curvature-matter coupling-modified theories of gravity. The paper is outlined in the following manner: In Section 2, we briefly introduce the linear curvature-matter coupling and present some of its interesting features. In Section 3, we generalize the latter linear curvature-matter coupling by considering the maximal extension of the Einstein–Hilbert action, which results in the $f\left(R,L_m\right)$ gravitational theory, and in Section 4, we extend the theory with the inclusion of general scalar field and kinetic term dependencies. In Section 5, we consider another extension of general relativity, namely $f(R,T)$ modified theories of gravity, where the gravitational Lagrangian is given by an arbitrary function of the Ricci scalar R and of the trace of the energy-momentum tensor T. In Section 6, we include an explicit invariant Ricci-energy-momentum tensor coupling and explore some of its astrophysical and cosmological phenomenology. This latter theory is motivated considering a traceless energy-momentum tensor, $T=0$; the gravitational field equations for the $f(R,T)$ theory reduce to that of the $f\left(R\right)$ gravity, and all non-minimal couplings of gravity to the matter field vanish, while the inclusion of the $R_{\mu \nu }{T}^{\mu \nu }$ term still allows a nonminimal coupling. The possibility of explaining dark matter as a consequence of the curvature-matter coupling is considered in Section 7. Finally, in Section 8, we briefly discuss our results and conclude our review paper.

2. Gravity Theories with Linear Curvature-Matter Coupling

In the present Section, we review some of the basic features and properties of the simplest class of models involving a nonminimal coupling between matter and geometry, which is linear in the matter Lagrangian. For the sake of comparison of the different approaches for the description of the gravitational interaction, we also briefly introduce the $f\left(R\right)$ modified gravity theory.

2.1. $f\left(R\right)$ Gravity

One of the simplest extensions of standard general relativity, based on the Hilbert–Einstein action, is the $f\left(R\right)$-modified theory of gravity, whose action takes the form [48,49]:

$$S=\int \left[\frac{1}{16\pi G}f\left(R\right)+L_m\right]\sqrt{-g}{d}^{4}x$$

(1)

where $f\left(R\right)$ is an arbitrary analytical function of the Ricci scalar R and $L_m$ is the Lagrangian density corresponding to matter.

Varying the action with respect to the metric $g_{\mu \nu }$ yields the field equations of $f\left(R\right)$ gravity:

$$F{R}_{\mu \nu }-\frac{1}{2}f{g}_{\mu \nu }-\left({\nabla }_{\mu }{\nabla }_{\nu }-g_{\mu \nu }\square \right)F=8\pi G{T}_{\mu \nu }$$

(2)

where we have denoted $F=df/dR$, and $T_{\mu \nu }$ is the standard, minimally coupled, matter energy-momentum tensor. Note that the covariant derivative of the field equations and of the matter energy-momentum tensor vanishes for all $f\left(R\right)$ by means of generalized Bianchi identities [9,50]. By contracting the field equations, Equation (2), we obtain the useful relation:

$$3\square F+FR-2f=8\pi GT$$

(3)

where T is the trace of the energy-momentum tensor and from which one verifies that the Ricci scalar is now a fully dynamical degree of freedom.

By introducing the Legendre transformation $\left\{R,f\right\}\to \left\{\varphi ,V\right\}$ defined as:

$$\varphi \equiv F\left(R\right),V\left(\varphi \right)\equiv R\left(\varphi \right)F-f\left(R\left(\varphi \right)\right)$$

(4)

the field equations of $f\left(R\right)$ gravity can be reformulated as [51,52,53,54]:

$$R_{\mu \nu }-\frac{1}{2}{g}_{\mu \nu }R=8\pi \frac{G}{\varphi }{T}_{\mu \nu }+{\theta }_{\mu \nu }$$

(5)

where:

$${\theta }_{\mu \nu }=-\frac{1}{2}V\left(\varphi \right)g_{\mu \nu }+\frac{1}{\varphi }\left({\nabla }_{\mu }{\nabla }_{\nu }-g_{\mu \nu }\square \right)\varphi$$

(6)

Using these variables, Equation (3) takes the following form:

$$3\square \varphi +2V\left(\varphi \right)-\varphi \frac{dV}{d\varphi }=8\pi GT$$

(7)

In this scalar-tensor representation, the field equations of $f\left(R\right)$ gravity can be derived from a Brans–Dicke-type gravitational action, with parameter $\omega =0$, given by:

$$S=\frac{1}{16\pi G}\int \left[\varphi R-V\left(\varphi \right)+L_m\right]\sqrt{-g}{d}^{4}x$$

(8)

The only requirement for the $f\left(R\right)$ model equations to be expressed in the form of a Brans–Dicke theory is that $F\left(R\right)$ be invertible, that is, $R\left(F\right)$ exists [55]. This condition is necessary for the construction of $V\left(\varphi \right)$.

The modification of the standard Einstein–Hilbert action leads to the appearance in the field equations of an effective gravitational constant $G_{\mathrm{eff}}=G/\varphi$, which is a function of the curvature. Secondly, a new source term for the gravitational field, given by the tensor ${\theta }_{\mu \nu }$, is also induced. The tensor ${\theta }_{\mu \nu }$ is determined by the trace of the energy-momentum tensor via Equation (7), which thus acts as an independent physical parameter determining the metric of the space-time.

2.2. Linear Non-Minimal Curvature-Matter Coupling

The action of the $f\left(R\right)$ modified gravity theory can be generalized by introducing in the action the linear nonminimal coupling between matter and geometry. The action of this modified gravity theory is given by [9]:

$$S=\int \left\{\frac{1}{2}{f}_1\left(R\right)+\left[1+\lambda f_2\left(R\right)\right]L_m\right\}\sqrt{-g}{d}^{4}x$$

(9)

where the factors $f_i\left(R\right)$ (with $i=1,2$) are arbitrary functions of the Ricci scalar R and $L_m$ is the matter Lagrangian density, and the coupling constant λ determines the strength of the interaction between $f_2\left(R\right)$ and the matter Lagrangian.

Now, varying the action with respect to the metric $g_{\mu \nu }$ provides the following field equations:

$$\begin{array}{c}\hfill F_1\left(R\right)R_{\mu \nu }-\frac{1}{2}{f}_1\left(R\right)g_{\mu \nu }-{\nabla }_{\mu }{\nabla }_{\nu }{F}_1\left(R\right)+g_{\mu \nu }\square F_1\left(R\right)=-2\lambda F_2\left(R\right)L_m{R}_{\mu \nu }\\ \hfill +2\lambda ({\nabla }_{\mu }{\nabla }_{\nu }-g_{\mu \nu }\square )L_m{F}_2\left(R\right)+[1+\lambda f_2\left(R\right)]T_{\mu \nu }^{\left(m\right)}\end{array}$$

(10)

where $F_i\left(R\right)=f_i^{\prime }\left(R\right)$, and the prime represents the derivative with respect to the scalar curvature. The matter energy-momentum tensor is defined as:

$$T_{\mu \nu }^{\left(m\right)}=-\frac{2}{\sqrt{-g}}\frac{\delta \left(\sqrt{-g}{L}_m\right)}{\delta \left(g^{\mu \nu }\right)}$$

(11)

Throughout this paper, we use the metric formalism. However, the field equations and the equations of motion were derived, for massive test particles in modified theories of gravity for a linear curvature-matter coupling, using the Palatini formalism in [56].

A general property of these nonminimal curvature-matter coupling theories is the non-conservation of the energy-momentum tensor. This can be easily verified by taking into account the covariant derivative of the field Equation (10), the Bianchi identities, ${\nabla }^{\mu }{G}_{\mu \nu }=0$, and the following identity $(\square {\nabla }_{\nu }-{\nabla }_{\nu }\square )F_i=R_{\mu \nu }{\nabla }^{\mu }{F}_i$, which then imply the following relationship:

$${\nabla }^{\mu }{T}_{\mu \nu }^{\left(m\right)}=\frac{\lambda F_2}{1+\lambda f_2}\left[g_{\mu \nu }{L}_m-T_{\mu \nu }^{\left(m\right)}\right]{\nabla }^{\mu }R$$

(12)

Thus, the coupling between the matter and the higher derivative curvature terms may be interpreted as an exchange of energy and momentum between both. Analogous couplings arise after a conformal transformation in the context of scalar-tensor theories of gravity and also in string theory. In the absence of the coupling, one verifies the conservation of the energy-momentum tensor [50], which can also be verified from the diffeomorphism invariance of the matter part of the action. It is also interesting to note that, from Equation (12), the conservation of the energy-momentum tensor is verified if $f_2\left(R\right)$ is a constant or the matter Lagrangian is not an explicit function of the metric.

In order to test the motion in our model, we consider for the energy-momentum tensor of matter a perfect fluid $T_{\mu \nu }^{\left(m\right)}=\left(\rho +p\right)U_{\mu }{U}_{\nu }+p{g}_{\mu \nu }$, where ρ is the overall energy density and p, the pressure, respectively. The four-velocity, $U_{\mu }$, satisfies the conditions $U_{\mu }{U}^{\mu }=-1$ and ${\nabla }_{\nu }{U}^{\mu }{U}_{\mu }=0$. We also introduce the projection operator $h_{\mu \lambda }=g_{\mu \lambda }+U_{\mu }{U}_{\lambda }$ from which one obtains $h_{\mu \lambda }{U}^{\mu }=0$.

By contracting Equation (12) with the projection operator $h_{\mu \lambda }$, one deduces the following expression:

$$\begin{array}{c}\hfill \left(\rho +p\right)g_{\mu \lambda }{U}^{\nu }{\nabla }_{\nu }{U}^{\mu }-\left({\nabla }_{\nu }p\right)({\delta }_{\lambda }^{\nu }-U^{\nu }{U}_{\lambda })-\frac{\lambda F_2}{1+\lambda f_2}\left(L_m-p\right)\left({\nabla }_{\nu }R\right)({\delta }_{\lambda }^{\nu }-U^{\nu }{U}_{\lambda })=0\end{array}$$

(13)

Finally, contraction with $g^{\alpha \lambda }$ gives rise to the equation of motion for a fluid element:

$$\frac{D{U}^{\alpha }}{ds}\equiv \frac{d{U}^{\alpha }}{ds}+{\Gamma }_{\mu \nu }^{\alpha }{U}^{\mu }{U}^{\nu }=f^{\alpha }$$

(14)

where the extra force is given by:

$$\begin{array}{ccc}\hfill f^{\alpha }& =& \frac{1}{\rho +p}\left[\frac{\lambda F_2}{1+\lambda f_2}\left(L_m-p\right){\nabla }_{\nu }R+{\nabla }_{\nu }p\right]h^{\alpha \nu }\hfill \end{array}$$

(15)

As one can immediately verify, the extra force $f^{\alpha }$ is orthogonal to the four-velocity of the particle, $f^{\alpha }{U}_{\alpha }=0$ which can be seen directly, from the properties of the projection operator. This is consistent with the usual interpretation of the force, according to which only the component of the four-force that is orthogonal to the particle’s four-velocity can influence its trajectory.

It has also been shown that an $f\left(G\right)$-modified Gauss–Bonnet gravity with a non-minimal coupling to matter also induces an extra force, which is normal to their four-velocities and, as a result, moves along nongeodesic world-lines [57].

A particularly intriguing feature is that the extra force depends on the form of the Lagrangian density. Note that considering the Lagrangian density $L_m=p$, where p is the pressure, the extra-force vanishes [58]. It has been argued that this is not the unique choice for the matter Lagrangian density and that more natural forms for $L_m$, such as $L_m=-\rho$, do not imply the vanishing of the extra-force. Indeed, in the presence of nonminimal coupling, they give rise to two distinct theories with different predictions [59], and this issue has been further investigated in different contexts [60,61], including phantom energy [62]. In this context, a matter Lagrangian density as an arbitrary function of the energy-density of the matter only was explored [63] (this possibility has also been explored in five dimensions [64] and in the braneworld context [65]). It was also argued that the corresponding energy-momentum tensor of the matter in modified gravity models with non-minimal coupling is more general than the usual general-relativistic energy-momentum tensor for perfect fluids [66], and it contains a supplementary, equation of state-dependent term, which could be related to the elastic stresses in the body or to other forms of internal energy. Therefore, the extra force induced by the coupling between matter and geometry never vanishes as a consequence of the thermodynamic properties of the system or for a specific choice of the matter Lagrangian, and it is non-zero in the case of a fluid of dust particles. In the following subsection, we discuss in detail the problem of the matter Lagrangian and of the energy-momentum tensor in modified gravity theories with linear curvature-matter coupling [66].

2.3. The Matter Lagrangian and the Energy-Momentum Tensor in Modified Gravity with Non-Minimal Linear Coupling between Matter and Geometry

It is an interesting and novel feature of the $f\left(R\right)$ gravity with non-minimal curvature-matter coupling that the matter Lagrangian, and the energy-momentum tensor obtained from it, are not model-independent quantities, but they are completely and uniquely determined by the coupling function between matter and geometry. This important result can be obtained by deriving first the equations of motion of test particles (or test fluid) in the modified gravity model from a variational principle and then considering the Newtonian limit of the particle action for a fluid obeying a barotropic equation of state [66]. The energy-momentum tensor of the matter obtained in this way is more general than the usual general-relativistic energy-momentum tensor for perfect fluids, and it contains a supplementary term that may be related to the elastic stresses in the test fluid or to other sources of internal energy. Since we assume that the matter obeys a barotropic equation of state, the matter Lagrangian can be expressed either in terms of the density or in terms of the pressure, and in both representations, the physical description of the system is equivalent. Therefore, the presence (or absence) of the extra-force is independent of the specific form of the matter Lagrangian, and it never vanishes, except in the particular case of (un)physical matter systems with zero sound speed. In particular, in the case of dust particles, the extra-force is always non-zero.

We define the energy-momentum tensor of the matter as [67]:

$$T_{\mu \nu }=-\frac{2}{\sqrt{-g}}\left[\frac{\partial \left(\sqrt{-g}{L}_m\right)}{\partial g^{\mu \nu }}-\frac{\partial }{\partial x^{\lambda }}\frac{\partial \left(\sqrt{-g}{L}_m\right)}{\partial \left(\partial g^{\mu \nu }/\partial x^{\lambda }\right)}\right]$$

(16)

By assuming that the Lagrangian density $L_m$ of the matter depends only on the metric tensor components $g_{\mu \nu }$, and not on its derivatives, we obtain:

$$T_{\mu \nu }=L_m{g}_{\mu \nu }-2\frac{\partial L_m}{\partial g^{\mu \nu }}$$

(17)

By taking into account the explicit form of the field equations Equation (10), one obtains for the covariant divergence of the energy-momentum tensor the equation:

$${\nabla }^{\mu }{T}_{\mu \nu }=2\left\{{\nabla }^{\mu }ln\left[1+\lambda f_2\left(R\right)\right]\right\}\frac{\partial L_m}{\partial g^{\mu \nu }}$$

(18)

As a specific example of generalized gravity models with linear curvature-matter coupling, we consider the case in which matter, assumed to be a perfect thermodynamic fluid, obeys a barotropic equation of state, with the thermodynamic pressure p being a function of the rest mass density of the matter (for short: matter density) ρ only, so that $p=p\left(\rho \right)$. In this case, the matter Lagrangian density, which, in the general case, could be a function of both density and pressure, $L_m=L_m\left(\rho ,p\right)$, or of only one of the thermodynamic parameters, becomes an arbitrary function of the density of the matter ρ only, so that $L_m=L_m\left(\rho \right)$. Then, the matter energy-momentum tensor is obtained as [66]:

$$T^{\mu \nu }=\rho \frac{d{L}_m}{d\rho }{U}^{\mu }{U}^{\nu }+\left(L_m-\rho \frac{d{L}_m}{d\rho }\right)g^{\mu \nu }$$

(19)

where the four-velocity $U^{\mu }=d{x}^{\mu }/ds$ satisfies the condition $g^{\mu \nu }{U}_{\mu }{U}_{\nu }=-1$. To obtain Equation (19), we have imposed the condition of the conservation of the matter current:

$${\nabla }_{\nu }\left(\rho U^{\nu }\right)=0$$

(20)

and we have used the relation:

$$\delta \rho =\frac{1}{2}\rho \left(g_{\mu \nu }-U_{\mu }{U}_{\nu }\right)\delta g^{\mu \nu }$$

(21)

whose proof is given in the Appendix of [66]. With the use of the identity $U^{\nu }{\nabla }_{\nu }{U}^{\mu }=d^2{x}^{\mu }/d{s}^2+{\Gamma }_{\nu \alpha }^{\mu }{U}^{\nu }{U}^{\alpha }$, from Equation (12) and Equation (19), we obtain the equation of motion of a massive test particle, or of a test fluid in the modified gravity model with linear coupling between matter and geometry, as:

$$\frac{d{U}^{\mu }}{ds}+{\Gamma }_{\nu \alpha }^{\mu }{U}^{\nu }{U}^{\alpha }=f^{\mu }$$

(22)

where:

$$f^{\mu }=-{\nabla }_{\nu }ln\left\{\left[1+\lambda f_2\left(R\right)\right]\frac{d{L}_m\left(\rho \right)}{d\rho }\right\}\left(U^{\mu }{U}^{\nu }-g^{\mu \nu }\right)$$

(23)

The extra-force $f^{\mu }$, generated due to the presence of the coupling between matter and geometry, is perpendicular to the four-velocity, $f^{\mu }{U}_{\mu }=0$.

The equation of motion Equation (22) can be obtained from the variational principle [66]:

$$\delta S_p=\delta \int L_{p}ds=\delta \int \sqrt{Q}\sqrt{g_{\mu \nu }{U}^{\mu }{U}^{\nu }}ds=0$$

(24)

where $S_p$ and $L_p=\sqrt{Q}\sqrt{g_{\mu \nu }{U}^{\mu }{U}^{\nu }}$ are the action and the Lagrangian density for test particles (test fluid), respectively, and:

$$\sqrt{Q}=\left[1+\lambda f_2\left(R\right)\right]\frac{d{L}_m\left(\rho \right)}{d\rho }$$

(25)

To prove this result, we start with the Lagrange equations corresponding to the action Equation (24):

$$\frac{d}{ds}\left(\frac{\partial L_p}{\partial U^{\alpha }}\right)-\frac{\partial L_p}{\partial x^{\alpha }}=0$$

(26)

Since:

$$\frac{\partial L_p}{\partial U^{\alpha }}=\sqrt{Q}{u}_{\alpha }$$

(27)

and:

$$\frac{\partial L_p}{\partial x^{\alpha }}=\frac{1}{2}\sqrt{Q}{g}_{\mu \nu ,\alpha }{U}^{\mu }{U}^{\nu }+\frac{1}{2}\frac{Q_{,\alpha }}{Q}$$

(28)

where a comma indicates the derivative with respect to $x^{\lambda }$, a straightforward calculation gives the equations of motion of the particle as:

$$\frac{d^2{x}^{\mu }}{d{s}^2}+{\Gamma }^{\mu }{}_{\nu \alpha }{U}^{\nu }{U}^{\alpha }+\left(U^{\mu }{U}^{\nu }-g^{\mu \nu }\right){\nabla }_{\nu }ln\sqrt{Q}=0$$

(29)

By simple identification with the equation of motion of the modified gravity model with linear curvature-matter coupling, given by Equation (22), we obtain the explicit form of $\sqrt{Q}$, as given by Equation (25). When $\sqrt{Q}\to 1$, we reobtain the standard general relativistic equation for geodesic motion.

2.3.1. The Newtonian Limit

The variational principle Equation (24) can be used to study the Newtonian limit of the model. In the limit of the weak gravitational fields:

$$ds\approx \sqrt{1+2\varphi -{\overrightarrow{v}}^2}dt\approx \left(1+\varphi -\frac{{\overrightarrow{v}}^2}{2}\right)dt$$

(30)

where ϕ is the Newtonian potential and $\overrightarrow{v}$ is the usual tridimensional velocity of the fluid. By representing the function $\sqrt{Q}$ as:

$$\sqrt{Q}=\frac{d{L}_m\left(\rho \right)}{d\rho }+\lambda f_2\left(R\right)\frac{d{L}_m\left(\rho \right)}{d\rho }$$

(31)

in the first order of approximation, the equations of motion of the fluid can be obtained from the variational principle:

$$\delta \int \left[\frac{d{L}_m\left(\rho \right)}{d\rho }+\lambda f_2\left(R\right)\frac{d{L}_m\left(\rho \right)}{d\rho }+\varphi -\frac{{\overrightarrow{v}}^2}{2}\right]dt=0$$

(32)

and are given by:

$$\overrightarrow{a}=-\nabla \varphi -\nabla \frac{d{L}_m\left(\rho \right)}{d\rho }-\nabla U_E={\overrightarrow{a}}_N+{\overrightarrow{a}}_H+{\overrightarrow{a}}_E$$

(33)

where $\overrightarrow{a}$ is the total acceleration of the system:

$${\overrightarrow{a}}_N=-\nabla \varphi$$

(34)

is the Newtonian gravitational acceleration, and:

$${\overrightarrow{a}}_E=-\nabla U_E=-\lambda \nabla \left[f_2\left(R\right)\frac{d{L}_m\left(\rho \right)}{d\rho }\right]$$

(35)

is a supplementary acceleration induced due to the curvature-matter coupling. As for the term:

$${\overrightarrow{a}}_H=-\nabla \left[\frac{d{L}_m\left(\rho \right)}{d\rho }\right]$$

(36)

it has to be identified with the hydrodynamic acceleration term in the perfect fluid Euler equation:

$${\overrightarrow{a}}_H=-\nabla \frac{d{L}_m\left(\rho \right)}{d\rho }=-\nabla {\int }_{{\rho }_0}^{\rho }\frac{dp}{d\rho }\frac{d\rho }{\rho }$$

(37)

where ${\rho }_0$, an integration constant, plays the role of a limiting density.

2.3.2. The Matter Lagrangian in Modified Gravity Theories with Curvature-Matter Coupling

With the use of Equation (37), the matter Lagrangian in modified gravity theories with curvature-matter coupling can be obtained by a simple integration as [66]:

$$L_m\left(\rho \right)=\rho \left[1+\Pi \left(\rho \right)\right]-{\int }_{p_0}^{p}dp$$

(38)

where:

$$\Pi \left(\rho \right)={\int }_{p_0}^p\frac{dp}{\rho }$$

(39)

and we have normalized an arbitrary integration constant to one; $p_0$ is an integration constant, or a limiting pressure. The corresponding energy-momentum tensor of matter is given by:

$$T^{\mu \nu }=\left\{\rho \left[1+\Phi \left(\rho \right)\right]+p\left(\rho \right)\right\}{U}^{\mu }{U}^{\nu }-p\left(\rho \right)g^{\mu \nu }$$

(40)

respectively, where:

$$\Phi \left(\rho \right)={\int }_{{\rho }_0}^{\rho }\frac{p}{{\rho }^2}d\rho =\Pi \left(\rho \right)-\frac{p\left(\rho \right)}{\rho }$$

(41)

and with all of the constant terms included in the definition of p. By introducing the energy density of the body according to the definition:

$$\epsilon =\rho \left[1+\Phi \left(\rho \right)\right]$$

(42)

the energy-momentum tensor of a test fluid can be written in the modified gravity models with curvature-matter coupling in a form similar to the standard general relativistic case:

$$T^{\mu \nu }=\left[\epsilon \left(\rho \right)+p\left(\rho \right)\right]U^{\mu }{U}^{\nu }+p\left(\rho \right)g^{\mu \nu }$$

(43)

From a physical point of view, $\Phi \left(\rho \right)$ can be interpreted as the elastic (deformation) potential energy of the body, and therefore, Equation (40) corresponds to the energy-momentum tensor of a compressible elastic isotropic system [68]. The matter Lagrangian can also be written in the simpler form:

$$L_m\left(\rho \right)=\rho \Phi \left(\rho \right)$$

(44)

If the pressure does not have a thermodynamic or radiative component, one can take $p_0=0$. If the pressure is a constant background quantity, independent of the density, so that $p=p_0$, then $L_m\left(\rho \right)=\rho$, and the matter energy-momentum tensor takes the form corresponding to dust:

$$T^{\mu \nu }=\rho U^{\mu }{U}^{\nu }$$

(45)

Since matter is supposed to obey a barotropic equation of state, these results are independent of the concrete representation of the matter Lagrangian in terms of the thermodynamic quantities [66]. The same results are obtained by assuming $L_m=L_m\left(p\right)$; due to the equation of state, ρ and p are freely interchangeable thermodynamic quantities, and the Lagrangians expressed in terms of ρ and p only are completely equivalent. More general situations, in which the density and pressure are functions of the particle number and temperature, respectively, and the equation of state is given in a parametric form, can be analyzed in a similar way.

The forms of the matter Lagrangian and the energy-momentum tensor are strongly dependent on the equation of state of the test fluid. For example, if the barotropic equation of state is linear, $p=\left(\gamma -1\right)\rho$, $\gamma =$ constant, $1\le \gamma \le 2$, then:

$$L_m\left(\rho \right)=\rho \left\{1+\left(\gamma -1\right)\left[ln\left(\frac{\rho }{{\rho }_0}\right)-1\right]\right\}$$

(46)

and $\Phi \left(\rho \right)=\left(\gamma -1\right)ln\left(\rho /{\rho }_0\right)$, respectively. For the case of a polytropic equation of state $p=K{\rho }^{1+1/n}$, $K,n=$ constant, we have:

$$L_m\left(\rho \right)=\rho +K\left(\frac{n^2}{n+1}-1\right){\rho }^{1+1/n}$$

(47)

and $\Phi \left(\rho \right)=Kn{\rho }^{1+1/n}=np\left(\rho \right)$, respectively, where we have taken for simplicity ${\rho }_0=p_0=0$. For a test fluid satisfying the ideal gas equation of state:

$$p=\frac{k_B}{\mu }\rho T$$

(48)

where $k_B$ is Boltzmann’s constant, T is the temperature and μ is the mean molecular weight, we obtain:

$$L_m\left(\rho \right)=\rho \left\{1+\frac{k_{B}T}{\mu }\left[ln\left(\frac{\rho }{{\rho }_0}\right)-1\right]\right\}+p_0$$

(49)

In the case of a physical system satisfying the ideal gas equation of state, the extra acceleration induced by the presence of the non-minimal coupling between matter and geometry is given by:

$${\overrightarrow{a}}_E\approx -\lambda \frac{k_{B}T}{\mu }\nabla \left[f_2\left(R\right)ln\frac{\rho }{{\rho }_0}\right]$$

(50)

and it is proportional to the temperature of the fluid. It is also interesting to note that the limiting density and pressure ${\rho }_0$ and $p_0$ generate in the energy-momentum tensor some extra constant terms, which may be interpreted as dark energy.

In conclusion, the extra-force induced by the coupling between matter and geometry does not vanish for any specific choices of the matter Lagrangian. In the case of the dust, with $p=0$, the extra force is given by:

$$f^{\mu }=-{\nabla }_{\nu }ln\left[1+\lambda f_2\left(R\right)\right]\left(U^{\mu }{U}^{\nu }-g^{\mu \nu }\right)$$

(51)

and it is independent of the thermodynamic properties of the system, being completely determined by the geometry, kinematics and coupling. In the limit of small velocities and weak gravitational fields, the extra-acceleration of a dust fluid is given by:

$${\overrightarrow{a}}_E=-\lambda \nabla \left[f_2\left(R\right)\right]$$

(52)

The thermodynamic condition for the vanishing of the extra-force is:

$$\frac{\partial L_m}{\partial g^{\mu \nu }}=\frac{1}{2}\left(\frac{\partial L_m}{\partial \rho }\right)\rho \left(g_{\mu \nu }-U_{\mu }{U}_{\nu }\right)=0$$

(53)

only. If the matter Lagrangian is written as a function of the pressure, then:

$$\frac{\partial L_m}{\partial \rho }=\left(\frac{\partial L_m}{\partial p}\right)\left(\frac{\partial p}{\partial \rho }\right)$$

(54)

and for all physical systems satisfying an equation of state (or, equivalently, for all systems with a non-zero sound velocity), the extra-force is non-zero. Therefore, the curvature-matter coupling is introduced in the generalized gravity models with a curvature-matter coupling in a consistent way. The coupling determines all of the physical properties of the system, including the extra-force, the matter Lagrangian and the energy-momentum tensor, respectively.

Hence, we have shown that in $f\left(R\right)$ gravity with a non-minimal coupling, the matter Lagrangian and the corresponding energy-momentum tensor are not model and thermodynamic parameters (independent quantities), but they are completely and uniquely determined by the nature of the coupling between matter and geometry, which, in the present model, is given by the function $f_2\left(R\right)$. We have obtained this result by deriving first the equations of motion in the modified gravity model with curvature-matter coupling from a variational principle and, then, by taking the Newtonian limit of the particle action for a fluid obeying a barotropic equation of state. The energy-momentum tensor of the matter obtained in this way is more general than the usual general-relativistic energy-momentum tensor for perfect fluids, and it contains a supplementary term that may be related to the elastic stresses in the body, or to other sources of internal energy. The matter Lagrangian can be expressed either in terms of the density or in terms of the pressure, and in both representations, the physical description of the system is equivalent. Therefore, the presence (or absence) of the extra-force is independent of the specific form of the matter Lagrangian, and it never vanishes, except in the case of (un)physical systems with zero sound speed. In particular, in the case of dust particles, the extra-force is always non-zero.

2.4. Equivalence of the Modified Gravity Theory with Linear Matter-Geometry Coupling with an Anomalous Scaler-Tensor Theory

As we have shown in Section 2.1, pure $f\left(R\right)$ gravity, with action given by Equation (1), is equivalent to a scalar-tensor theory, with action given by Equation (8). The equivalence between the modified gravity models with linear geometry-matter coupling was established by Faraoni [69]. In the following, we show that the action Equation (9) is equivalent with a scalar-tensor Brans–Dicke-type theory, with a single scalar field, a vanishing Brans–Dicke parameter ω and an unusual coupling of the potential $U\left(\psi \right)$ of the theory to matter.

By introducing a new field ϕ, the action Equation (9) becomes:

$$S\int d^{4}x\sqrt{-g}\left\{\frac{f_1\left(\varphi \right)}{2}+\frac{1}{2}\frac{d{f}_1}{d\varphi }\left(R-\varphi \right)+\left[1+\lambda f_2\left(\varphi \right)\right]{\mathcal{L}}_m\right\}$$

(55)

Next, we further introduce the field $\psi \left(\varphi \right)\equiv f_1^{\prime }\left(\varphi \right)$ (with a prime denoting the differentiation with respect to ϕ), and we obtain for the action the expression:

$$S=\int d^{4}x\sqrt{-g}\left[\frac{\psi R}{2}-V\left(\psi \right)+U\left(\psi \right){\mathcal{L}}_m\right]$$

(56)

where:

$$\begin{array}{ccc}\hfill V\left(\psi \right)& =& \frac{\varphi \left(\psi \right)f_1^{\prime }\left[\varphi \left(\psi \right)\right]-f_1\left[\varphi \left(\psi \right)\right]}{2}\hfill \end{array}$$

(57)

$$\begin{array}{ccc}\hfill U\left(\psi \right)& =& 1+\lambda f_2\left[\varphi \left(\psi \right)\right]\hfill \end{array}$$

(58)

$\varphi \left(\psi \right)$ must be obtained by inverting $\psi \left(\varphi \right)\equiv f_1^{\prime }\left(\varphi \right)$. The actions Equation (9) and Equation (56) are equivalent when $f_1^{\prime \prime }\left(R\right)\ne 0$[69]. We can see this by setting $\varphi =R$, and then Equation (56) reduces trivially to Equation (9). On the other hand, the variation of Equation (55) with respect to ϕ gives:

$$\left(R-\varphi \right)f_1^{\prime \prime }\left(\varphi \right)+2\lambda f_2^{\prime }\left(\varphi \right){\mathcal{L}}_m=0$$

(59)

In a vacuum, we have ${\mathcal{L}}_m=0$, and Equation (59) gives $\varphi =R$ whenever $f_1^{\prime \prime }\ne 0$ [51,52,53,54]. However, in the presence of matter, there seem to be other possibilities, which, however, can be excluded as follows.

When ${\mathcal{L}}_m\ne 0$, the action Equation (9) and Equation (55) are equivalent if:

$$\left(R-\varphi \right)f_1^{\prime \prime }\left(\varphi \right)+2\lambda f_2^{\prime \prime }\left(\varphi \right){\mathcal{L}}_m\ne 0$$

(60)

When Equation (59) is satisfied, we have a pathological case. It corresponds to:

$$\lambda f_2\left(\varphi \right){\mathcal{L}}_m=\frac{f_1^{\prime }\left(\varphi \right)}{2}\left(\varphi -R\right)-\frac{f_1\left(\varphi \right)}{2}$$

(61)

However, if Equation (61) holds, then the action Equation (55) reduces to the trivial case of pure matter without the gravity sector. Then, it follows that in the modified gravity theories with linear geometry-matter coupling, the actions Equation (9) and Equation (56) are equivalent when $f_1^{\prime \prime }\left(R\right)\ne 0$, similar to the case of pure $f\left(R\right)$ gravity [51,52,53,54].

2.5. Further Theoretical Developments in Modified Gravity with Linear Curvature-Matter Coupling

This linear nonminimal curvature-matter coupling has been extensively explored in a plethora of contexts in the literature. For instance, the equations of motion of test bodies for a theory with nonminimal coupling by means of a multipole method [70] was also studied, and it was shown that the propagation equations for pole-dipole particles allow for a systematic comparison with the equations of motion of general relativity and other gravity theories.

The consequences that a non-minimal coupling between curvature and matter can have on the dynamics of perfect fluids has also been investigated [71]. It was argued that the presence of a static, axially-symmetric, pressureless fluid does not imply an asymptotically Minkowski space-time, such as in General Relativity. This feature can be attributed to a pressure mimicking mechanism related to the non-minimal coupling. The case of a spherically symmetric black hole surrounded by fluid matter was analyzed, and it was shown that under equilibrium conditions, the total fluid mass is about twice that of the black hole [71].

The consequences of the curvature-matter coupling on stellar equilibrium and constraints on the nonminimal coupling was considered, where particular attention was paid to the validity of the Newtonian regime and on the boundary and exterior matching conditions [72]. This explicit “anomalous” coupling of the Ricci curvature to matter has also raised the question of curvature instabilities, and in [69,73,74,75,76], constraints imposed by the energy condition and the conditions in order to avoid the notorious Dolgov–Kawasaki instability [29] were obtained. It has also been claimed that the curvature-matter action leads to a theory of gravity that includes higher order derivatives of the matter fields without introducing more dynamics in the gravity sector and, therefore, cannot be a viable theory for gravitation [77,78]. However, we emphasize that the results of [77,78] only apply to the specific case of $f_2\left(R\right)=R$ and not to action Equation (9) in general. For more generic functional forms, the theory propagates extra degrees of freedom and the conclusions of [77,78] do not apply. The relation between these theories and ordinary scalar-tensor gravity was also analyzed, as well as its implications for the equivalence principle [43].

In fact, the theoretical consistency of these nonminimal curvature-matter couplings was studied using a scalar field Lagrangian to model the matter content [79]. The conditions that the coupling does not introduce ghosts, classical instabilities or superluminal propagation of perturbations were derived. These consistency conditions were then employed to rule out or severely restrict the forms of the non-minimal coupling functions [79]. Lagrange–Noether methods have been used to derive the conservation laws for models in which matter interacts non-minimally with the gravitational field [80]. Furthermore, a covariant derivation of the equations of motion for test bodies for a wide class of gravitational theories with nonminimal coupling was presented in [81], encompassing a general interaction via the complete set of nine parity-even curvature invariants. The equations of motion for spinning test bodies in such theories were explicitly derived by means of Synge’s expansion technique. The authors’ findings generalize previous results in the literature and allowed for a direct comparison to the general relativistic equations of motion of pole-dipole test bodies. In [82], the authors derived multipolar equations of motion for gravitational theories with general nonminimal coupling in spacetimes admitting torsion. Their general findings allow for a systematic testing of whole classes of theories by means of extended test bodies. One peculiar feature of certain subclasses of nonminimal theories turns out to be their sensitivity to post-Riemannian spacetime structures, even in experiments without microstructured test matter. Weak field constraints have also been studied in detail [83,84], and wormhole solutions were further explored, where it is the higher curvature coupling terms that support these exotic geometries [85,86,87].

In a cosmological context, the perturbation equation of matter on subhorizon scales was deduced, and specific bounds on the theory from weak lensing observations and the primordial nucleosynthesis were obtained in order to constrain the parameters of the model [88]. It was also shown that a non-minimal coupling between the scalar curvature and the matter Lagrangian density may account for the accelerated expansion of the Universe [89,90,91] and provide, through mimicking, for a viable unification of dark energy and dark matter [89]. It was shown that a generalized non-minimal coupling between curvature and matter is compatible with Starobinsky inflation and leads to a successful process of preheating [92], and the problem of a cosmological constant was further explored [93].

The effects of the non-minimal curvature-matter coupling on the evolution of cosmological perturbations around a homogeneous and isotropic Universe and, hence, the formation of large-scale structure have also been analyzed [94]. This framework places constraints on the terms, which arise due to the coupling with matter and, in particular, on the modification in the growth of matter density perturbations. Approximate analytical solutions were obtained for the evolution of matter overdensities during the matter dominated era, and it was shown that these favor the presence of a coupling function that is compatible with the late-time cosmic acceleration.

The observations related to the growth of matter has also shown that there is a small, but finite window, where one can distinguish the non-minimally-coupled $f\left(R\right)$ models with the concordance ΛCDM [95]. The possibility that the behavior of the rotational velocities of test particles gravitating around galaxies can be explained in the framework of modified gravity models with nonminimal curvature-matter coupling has also been extensively explored [42,96,97].

In fact, the literature is extremely vast, and rather than enumerate all of the features of these models, we refer the reader to [98,99] for a review on the topic of the linear nonminimal curvature-matter coupling.

3. $\mathbf{f}\left(\mathbf{R},{\mathbf{L}}_{\mathbf{m}}\right)$ Gravity

In this section, we generalize the $f\left(R\right)$-type gravity models by assuming that the gravitational Lagrangian is given by an arbitrary function of the Ricci scalar R and of the matter Lagrangian $L_m$ [15]. This consists in a maximal extension of the Hilbert–Einstein action, and the action takes the following form:

$$S=\int f\left(R,L_m\right)\sqrt{-g}{d}^{4}x$$

(62)

where $f\left(R,L_m\right)$ is an arbitrary function of the Ricci scalar R, and of the Lagrangian density corresponding to matter, $L_m$. The energy-momentum tensor of the matter is defined by Equation (11). Thus, by assuming that the Lagrangian density $L_m$ of the matter depends only on the metric tensor components $g_{\mu \nu }$, and not on its derivatives, we obtain $T_{\mu \nu }=g_{\mu \nu }{L}_m-2\partial L_m/\partial g^{\mu \nu }$, which will be useful below.

Now, varying the action with respect to the metric yields the following field equation:

$$\begin{array}{ccc}& & f_R\left(R,L_m\right)R_{\mu \nu }+\left(g_{\mu \nu }\square -{\nabla }_{\mu }{\nabla }_{\nu }\right)f_R\left(R,L_m\right)\hfill \\ & & -\frac{1}{2}\left[f\left(R,L_m\right)-f_{L_m}\left(R,L_m\right)L_m\right]g_{\mu \nu }=\frac{1}{2}{f}_{L_m}\left(R,L_m\right)T_{\mu \nu }\hfill \end{array}$$

(63)

If $f\left(R,L_m\right)=R/2+L_m$ (the Hilbert–Einstein Lagrangian), we recover the standard Einstein field equation of general relativity, $R_{\mu \nu }-(1/2)g_{\mu \nu }R=T_{\mu \nu }$. For $f\left(R,L_m\right)=f_1\left(R\right)+f_2\left(R\right)G\left(L_m\right)$, where $f_1$ and $f_2$ are arbitrary functions of the Ricci scalar and G a function of the matter Lagrangian density, respectively, we reobtain the field equations of the modified gravity with arbitrary curvature-matter coupling, considered in [63].

The contraction of Equation (63) provides the following relation between the Ricci scalar R, the matter Lagrange density $L_m$ and the trace $T=T_{\mu }^{\mu }$ of the energy-momentum tensor:

$$\begin{array}{c}\hfill f_R\left(R,L_m\right)R+3\square f_R\left(R,L_m\right)-2\left[f\left(R,L_m\right)-f_{L_m}\left(R,L_m\right)L_m\right]=\frac{1}{2}{f}_{L_m}\left(R,L_m\right)T\end{array}$$

(64)

By eliminating the term ${\nabla }_{\mu }{\nabla }^{\mu }{f}_R\left(R,L_m\right)$ between Equation (63) and Equation (64), we obtain another form of the gravitational field equations as:

$$\begin{array}{ccc}& & f_R\left(R,L_m\right)\left(R_{\mu \nu }-\frac{1}{3}R{g}_{\mu \nu }\right)+\frac{1}{6}\left[f\left(R,L_m\right)-f_{L_m}\left(R,L_m\right)L_m\right]g_{\mu \nu }=\hfill \\ & & \frac{1}{2}{f}_{L_m}\left(R,L_m\right)\left(T_{\mu \nu }-\frac{1}{3}T{g}_{\mu \nu }\right)+{\nabla }_{\mu }{\nabla }_{\nu }{f}_R\left(R,L_m\right)\hfill \end{array}$$

(65)

By taking the covariant divergence of Equation (63), with the use of the mathematical identity [50]:

$$\begin{array}{c}\hfill {\nabla }^{\mu }\left[f_R\left(R,L_m\right)R_{\mu \nu }-\frac{1}{2}f\left(R,L_m\right)g_{\mu \nu }+\left(g_{\mu \nu }\square -{\nabla }_{\mu }{\nabla }_{\nu }\right)f_R\left(R,L_m\right)\right]\equiv 0\end{array}$$

(66)

we obtain for the divergence of the energy-momentum tensor $T_{\mu \nu }$, the following equation:

$$\begin{array}{ccc}\hfill {\nabla }^{\mu }{T}_{\mu \nu }& =& {\nabla }^{\mu }ln\left[f_{L_m}\left(R,L_m\right)\right]\left\{L_m{g}_{\mu \nu }-T_{\mu \nu }\right\}\hfill \\ & =& 2{\nabla }^{\mu }ln\left[f_{L_m}\left(R,L_m\right)\right]\frac{\partial L_m}{\partial g^{\mu \nu }}\hfill \end{array}$$

(67)

The requirement of the conservation of the energy-momentum tensor of matter, ${\nabla }^{\mu }{T}_{\mu \nu }=0$, gives an effective functional relation between the matter Lagrangian density and the function $f_{L_m}\left(R,L_m\right)$:

$${\nabla }^{\mu }ln\left[f_{L_m}\left(R,L_m\right)\right]\frac{\partial L_m}{\partial g^{\mu \nu }}=0$$

(68)

Thus, once the matter Lagrangian density is known, by an appropriate choice of the function $f\left(R,L_m\right)$, one can construct, at least in principle, conservative models with arbitrary curvature-matter dependence.

Now, assuming that the matter Lagrangian is a function of the rest mass density ρ of the matter only, from Equation (67), we obtain explicitly the equation of motion of the test particles in the $f\left(R,L_m\right)$ gravity model as:

$$\frac{D^2{x}^{\mu }}{d{s}^2}=U^{\nu }{\nabla }_{\nu }{U}^{\mu }=\frac{d^2{x}^{\mu }}{d{s}^2}+{\Gamma }_{\nu \lambda }^{\mu }{U}^{\nu }{U}^{\lambda }=f^{\mu }$$

(69)

where the world-line parameter s is taken as the proper time, $U^{\mu }=d{x}^{\mu }/ds$ is the four-velocity of the particle, ${\Gamma }_{\sigma \beta }^{\nu }$ are the Christoffel symbols associated with the metric and the extra-force $f^{\mu }$ is defined as:

$$f^{\mu }=-{\nabla }_{\nu }ln\left[f_{L_m}\left(R,L_m\right)\frac{d{L}_m\left(\rho \right)}{d\rho }\right]\left(U^{\mu }{U}^{\nu }-g^{\mu \nu }\right)$$

(70)

Note that as in the linear curvature-matter coupling, the extra-force $f^{\mu }$, generated by the curvature-matter coupling, is perpendicular to the four-velocity, $f^{\mu }{U}_{\mu }=0$. Due to the presence of the extra-force $f^{\mu }$, the motion of the test particles in modified theories of gravity with an arbitrary coupling between matter and curvature is non-geodesic. From the relation $U_{\mu }{\nabla }_{\nu }{U}^{\mu }\equiv 0$, it follows that the force $f^{\mu }$ is always perpendicular to the velocity, so that $U_{\mu }{f}^{\mu }=0$.

Building on this analysis, the geodesic deviation equation, describing the relative accelerations of nearby particles, and the Raychaudhury equation, giving the evolution of the kinematical quantities associated with deformations (expansion, shear and rotation) were extensively considered in this framework of modified theories of gravity with an arbitrary curvature-matter coupling, by taking into account the effects of the extra force [100,101]. As a physical application of the geodesic deviation equation, the modifications of the tidal forces due to the supplementary curvature-matter coupling were obtained in the weak field approximation. The tidal motion of test particles is directly influenced not only by the gradient of the extra force, which is basically determined by the gradient of the Ricci scalar, but also by an explicit coupling between the velocity and the Riemann curvature tensor. As a specific example, the expression of the Roche limit (the orbital distance at which a satellite will begin to be tidally torn apart by the body it is orbiting) was also obtained for this class of models. These aspects will be presented below.

Furthermore, the energy conditions and cosmological applications were explored [102]. In [103], the Wheeler–DeWitt equation of $f(R,L_m)$ gravity was analysed in a flat FRW (Friedmann-Robertson-Walker) universe, which is the first step of the study of quantum $f(R,L_m)$ cosmology. In the minisuperspace spanned by the FRW scale factor and the Ricci scalar, the equivalence of the reduced action was examined, and the canonical quantization of the $f(R,L_m)$ model was undertaken and the corresponding Wheeler–DeWitt equation derived. The introduction of the invariant contractions of the Ricci and Riemann tensors were further considered, and applications to black hole and wormhole physics were analyzed [104].

3.1. Solar System Tests of $f\left(R,L_m\right)$ Gravity

One of the basic predictions of the modified gravity theories with a curvature-matter coupling is the existence of an extra force, which makes the motion of the test particles non-geodesic. The existence of this force can also be tested at the level of the Solar System, by estimating its effects on the orbital parameters of the motion of the planets around the Sun. The impact on the planetary motion of the extra force can be obtained in a simple way by using the properties of the Runge–Lenz vector, defined as:

$$\overrightarrow{A}=\overrightarrow{v}\times \overrightarrow{L}-\alpha {\overrightarrow{e}}_r$$

(71)

where $\overrightarrow{v}$ is the velocity of the planet of mass m relative to the Sun, with mass $M_{\odot }$, $\overrightarrow{r}=r{\overrightarrow{e}}_r$ the two-body position vector, $\overrightarrow{p}=\mu \overrightarrow{v}$ the relative momentum, $\mu =m{M}_{\odot }/\left(m+M_{\odot }\right)$ the reduced mass:

$$\overrightarrow{L}=\overrightarrow{r}\times \overrightarrow{p}=\mu r^2\dot{\theta }\overrightarrow{k}$$

(72)

is the angular momentum, and $\alpha =Gm{M}_{\odot }$, respectively (see [63] and the references therein). For an elliptical orbit of eccentricity e, major semi-axis a and period T, the equation of the orbit is given by $\left(L^2/\mu \alpha \right)r^{-1}=1+ecos\theta$. The Runge–Lenz vector can be expressed as:

$$\overrightarrow{A}=\left(\frac{{\overrightarrow{L}}^2}{\mu r}-\alpha \right){\overrightarrow{e}}_r-\dot{r}L{\overrightarrow{e}}_{\theta }$$

(73)

and its derivative with respect to the polar angle θ is given by $d\overrightarrow{A}/d\theta =r^2\left[dV\left(r\right)/dr-\alpha /r^2\right]{\overrightarrow{e}}_{\theta }$, where $V\left(r\right)$ is the potential of the central force [63].

The gravitational potential term acting on a planet consists of the Post–Newtonian potential:

$$V_{PN}\left(r\right)=-\frac{\alpha }{r}-3\frac{{\alpha }^2}{m{r}^2}$$

(74)

plus the gravitational term induced by the general coupling between matter and geometry. Thus, we have:

$$\frac{d\overrightarrow{A}}{d\theta }=r^2\left[6\frac{{\alpha }^2}{m{r}^3}+m{\overrightarrow{a}}_E\left(r\right)\right]{\overrightarrow{e}}_{\theta }$$

(75)

where we have also assumed that $\mu \approx m$. Then, we obtain the change of $\Delta \varphi$ of the perihelion with a change of θ of $2\pi$ as:

$$\Delta \varphi =\left(\frac{1}{\alpha e}\right){\int }_0^{2\pi }\left|\dot{\overrightarrow{L}}\times d\frac{\overrightarrow{A}}{d\theta }\right|d\theta$$

(76)

which can be explicitly calculated as:

$$\Delta \varphi =24{\pi }^3{\left(\frac{a}{T}\right)}^2\frac{1}{1-e^2}+\frac{L}{8{\pi }^{3}me}\frac{{\left(1-e^2\right)}^{3/2}}{{\left(a/T\right)}^3}{\int }_0^{2\pi }\frac{a_E\left[L^2{\left(1+ecos\theta \right)}^{-1}/m\alpha \right]}{{\left(1+ecos\theta \right)}^2}cos\theta d\theta$$

(77)

where we have used the relation $\alpha /L=2\pi \left(a/T\right)/\sqrt{1-e^2}$. The first term of Equation (77) gives the expression of the standard general relativistic precession of the perihelion of the planets, while the second term gives the contribution to the perihelion precession due to the presence of the extra force generated by the coupling between matter and curvature.

As an example of the application of Equation (77), we consider the simple case for which the extra force may be considered as a constant, $a_E\approx$ constant, an approximation that could be valid for small regions of space-time. This case may also correspond to a MOND-type acceleration $a_E\approx \sqrt{a_0{a}_N}=\sqrt{G{M}_{\odot }{a}_0}/r$, where $a_0$ is a constant acceleration, which was proposed phenomenologically as a dynamical model for dark matter [105,106]. In the Newtonian limit, the extra-acceleration generated by the curvature-matter coupling can be expressed in a similar form [9]. With the use of Equation (77), one finds for the perihelion precession of a planet in the Solar System the expression:

$$\Delta \varphi =\frac{6\pi G{M}_{\odot }}{a\left(1-e^2\right)}+\frac{2\pi a^2\sqrt{1-e^2}}{G{M}_{\odot }}{a}_E$$

(78)

where we have also used Kepler’s third law, $T^2=4{\pi }^2{a}^3/G{M}_{\odot }$. For the planet Mercury $a=57.91\times {10}^{11}$ cm, and $e=0.205615$, respectively, while $M_{\odot }=1.989\times {10}^{33}$ g [19]. With the use of these numerical values, the first term in Equation (78) gives the standard general relativistic value for the precession angle, ${\left(\Delta \varphi \right)}_{GR}=42.962$ arcsec per century. On the other hand, the observed value of the precession of the perihelion of Mercury is ${\left(\Delta \varphi \right)}_{obs}=43.11\pm 0.21$ arcsec per century [19]. Therefore, the difference ${\left(\Delta \varphi \right)}_E={\left(\Delta \varphi \right)}_{obs}-{\left(\Delta \varphi \right)}_{GR}=0.17$ arcsec per century can be attributed to other physical effects. Hence, the observational constraint on the perihelion precession of the planet Mercury requires that the value of the constant acceleration $a_E$ induced by the curvature-matter coupling at the scale of the Solar System must satisfy the condition $a_E\le 1.28\times {10}^{-9}$ cm/s${}^2$. This value of $a_E$, obtained from the high precision Solar System observations, is somewhat smaller than the value of the extra acceleration $a_0\approx {10}^{-8}$ cm/s${}^2$, necessary to explain the “dark matter” properties, as well as the Pioneer anomaly [9]. However, this value does not rule out the possibility of the presence of some extra curvature-matter coupling-type gravitational effects, acting at both the Solar System and galactic levels, since the assumption of a constant extra acceleration may not be correct on larger astronomical scales.

3.2. The Geodesic Deviation Equation and the Raychaudhury Equation in $f\left(R,L_m\right)$ Gravity

In the present subsection, we present the equation of the geodesic deviation and the Raychaudhury equation in $f\left(R,L_m\right)$ gravity, which explicitly contain the effects of the curvature-matter coupling and of the extra force. Some of the physical implications of the geodesic deviation equation, namely, the problem of the tidal forces in $f\left(R,L_m\right)$ gravity is also considered, and the generalization of the Roche limit will be considered in the next subsection. In our presentation, we closely follow [100,101].

Consider a one-parameter congruence of curves $x^{\mu }\left(s;\lambda \right)$, so that for each $\lambda ={\lambda }_0=$ constant, $x^{\mu }\left(s,{\lambda }_0\right)$ satisfies Equation (69). We suppose the parametrization to be smooth, and hence, we can introduce the tangent vector fields along the trajectories of the particles as $U^{\mu }=\partial x^{\mu }\left(s;\lambda \right)/\partial s$ and $n^{\mu }=\partial x^{\mu }\left(s;\lambda \right)/\partial \lambda$, respectively. We also introduce the four-vector:

$${\eta }^{\mu }=\left[\frac{\partial x^{\mu }\left(s;\lambda \right)}{\partial \lambda }\right]\delta \lambda \equiv n^{\mu }\delta \lambda$$

(79)

joining points on infinitely close geodesics, corresponding to parameter values λ and $\lambda +\delta \lambda$, which have the same value of s. By taking into account Equation (69), we obtain the geodesic deviation equation (Jacobi equation), giving the second order derivative with respect to the parameter s of the deviation vector ${\eta }^{\mu }$ as [67,107]:

$$\frac{D^2{\eta }^{\mu }}{d{s}^2}=R_{\nu \alpha \beta }^{\mu }{\eta }^{\alpha }{U}^{\beta }{U}^{\nu }+{\eta }^{\alpha }{\nabla }_{\alpha }{f}^{\mu }$$

(80)

In the case $f^{\mu }\equiv 0$, we reobtain the standard Jacobi equation, corresponding to the geodesic motion of test particles in standard general relativity. The interest in the deviation vector ${\eta }^{\mu }$ derives from the fact that if $x_0^{\mu }\left(s\right)=x^{\mu }\left(s;{\lambda }_0\right)$ is a solution of Equation (80), then to first order $x_1^{\mu }\left(s\right)=x_0^{\mu }\left(s\right)+{\eta }^{\mu }$ is a solution of the geodesic equation, as well, since $x^{\mu }\left(s;{\lambda }_1\right)\approx x^{\mu }\left(s;{\lambda }_0\right)+n^{\mu }\left(s;{\lambda }_0\right)\delta \lambda \approx x^{\mu }\left(s;{\lambda }_0+\delta \lambda \right)$.

By taking into account the explicit form of the extra force given by Equation (70), in $f\left(R,L_m\right)$-modified theories of gravity with a curvature-matter coupling, the geodesic deviation equation can be written as:

$$\begin{array}{c}\hfill \frac{D^2{\eta }^{\mu }}{d{s}^2}=R_{\nu \alpha \beta }^{\mu }{\eta }^{\alpha }{U}^{\beta }{U}^{\nu }+{\eta }^{\alpha }{\nabla }_{\alpha }\left\{{\nabla }_{\nu }ln\left[f_{L_m}\left(R,L_m\right)\frac{d{L}_m\left(\rho \right)}{d\rho }\right]\left(L_m{g}^{\mu \nu }-T^{\mu \nu }\right)\right\}\end{array}$$

(81)

Explicitly, the geodesic deviation equation becomes:

$$\begin{array}{ccc}\hfill \frac{D^2{\eta }^{\mu }}{d{s}^2}& =& R_{\nu \alpha \beta }^{\mu }{\eta }^{\alpha }{U}^{\beta }{U}^{\nu }+{\eta }^{\alpha }\left\{{\nabla }_{\alpha }{\nabla }_{\nu }ln\left[f_{L_m}\left(R,L_m\right)\frac{d{L}_m\left(\rho \right)}{d\rho }\right]\right\}\left(L_m{g}^{\mu \nu }-T^{\mu \nu }\right)\hfill \\ & & +{\eta }^{\alpha }{\nabla }_{\nu }ln\left[f_{L_m}\left(R,L_m\right)\frac{d{L}_m\left(\rho \right)}{d\rho }\right]\left(g^{\mu \nu }{\nabla }_{\alpha }{L}_m-{\nabla }_{\alpha }{T}^{\mu \nu }\right)\hfill \end{array}$$

(82)

In the presence of an extra force, the Raychaudhury equation is obtained as [107]:

$$\dot{\theta }+\frac{1}{3}{\theta }^2+\left({\sigma }^2-{\omega }^2\right)={\nabla }_{\mu }{f}^{\mu }+R_{\mu \nu }{U}^{\mu }{U}^{\nu }$$

(83)

where $\theta ={\nabla }_{\nu }{U}^{\nu }$ is the expansion of the congruence of particles, ${\sigma }^2={\sigma }_{\mu \nu }{\sigma }^{\mu \nu }$ and ${\omega }^2={\omega }_{\mu \nu }{\omega }^{\mu \nu }$, respectively.

With the explicit use of the field equation Equation (63) and the expression of the extra force given by Equation (70), in $f\left(R,L_m\right)$-modified theories of gravity with an arbitrary coupling between curvature and matter, the Raychaudhury equation assumes the following generalised form:

$$\begin{array}{ccc}\hfill \dot{\theta }& =& -\frac{1}{3}{\theta }^2-\left({\sigma }^2-{\omega }^2\right)+\Lambda \left(R,L_m\right)+{\nabla }_{\mu }\left\{{\nabla }_{\nu }ln\left[f_{L_m}\left(R,L_m\right)\frac{d{L}_m\left(\rho \right)}{d\rho }\right]\left(L_m{g}^{\mu \nu }-T^{\mu \nu }\right)\right\}\hfill \\ & & +\frac{1}{f_R\left(R,L_m\right)}{U}^{\mu }{U}^{\nu }{\nabla }_{\mu }{\nabla }_{\nu }{f}_R\left(R,L_m\right)+\Phi \left(R,L_m\right)\left(T_{\mu \nu }{U}^{\mu }{U}^{\nu }-\frac{1}{3}T\right)\hfill \end{array}$$

(84)

where we have denoted:

$$\Lambda \left(R,L_m\right)=\frac{2f_R\left(R,L_m\right)R-f\left(R,L_m\right)+f_{L_m}\left(R,L_m\right)L_m}{6f_R\left(R,L_m\right)}$$

(85)

and:

$$\Phi \left(R,L_m\right)=\frac{f_{L_m}\left(R,L_m\right)}{f_R\left(R,L_m\right)}$$

(86)

respectively.

3.3. Tidal Forces and the Roche Limit

Tides are common astrophysical phenomena, and they are due to the presence of a gradient of the gravitational force field induced by a mass above an extended body or a system of particles. In the Solar System, tidal perturbations act on compact celestial bodies, such as planets, moons and comets. On larger scales than the Solar System, as, for example, in a galactic or cosmological context, one can observe tidal deformations or disruptions of a stellar cluster by a galaxy or in galaxy encounters [108,109,110]. In the relativistic theories of gravitation, as well as in Newtonian gravity, a local system of coordinates can be chosen, which is inertial, except for the presence of the tidal forces. In strong gravitational fields, relativistic tidal effects can lead to important physical phenomena, such as the emission of tidal gravitational waves [108,109,110].

In the following, we denote by a prime the reference frame in which all of the Christoffel symbols vanish. In such a system, one can always take the deviation vector component ${\eta }^{\prime 0}=0$, which means that the particle accelerations are compared at equal times. ${\eta }^{\prime i}$ is then the displacement of the particle from the origin [111]. Moreover, in the static/stationary case, in which the metric, the Ricci scalar and the thermodynamic parameters of the matter do not depend on time, we have $f^{\prime 0}=0$. With the use of the equation of motion, this condition implies $U^{\prime 0}=$ constant $=1$. Therefore, with these assumptions, in a static or stationary space-time, the equation of the geodesic deviation (the Jacobi equation) takes the form [100]:

$$\frac{d^2{\eta }^{\prime i}}{d{t}^{\prime 2}}=R_{0l0}^{\prime i}{\eta }^{\prime l}+R_{jlm}^{\prime i}{\eta }^{\prime l}{U}^{\prime j}{U}^{{}^{\prime }m}+{\eta }^{\prime l}\frac{\partial f^{\prime i}}{\partial x^{\prime l}}$$

(87)

Equation (87) can be reformulated as:

$$F^{\prime i}=\frac{d^2{\eta }^{\prime i}}{d{t}^{\prime 2}}=K_j^i{\eta }^{\prime j}$$

(88)

where $F^{\prime i}$ is the tidal force, and we have also introduced the generalized tidal matrix $K_l^i$ [108,109,110], which is defined as:

$$K_j^i=R_{0j0}^{\prime i}+R_{kjm}^{\prime i}{U}^{\prime k}{U}^{{}^{\prime }m}+\frac{\partial f^{\prime i}}{\partial x^{\prime j}}$$

(89)

The tidal force has the property:

$$\frac{\partial F^{\prime i}}{\partial {\eta }^{\prime j}}=K_j^i$$

(90)

and its divergence is given by:

$$\frac{\partial F^{\prime i}}{\partial {\eta }^{\prime i}}=K$$

(91)

where the trace K of the tidal matrix is:

$$K=K_i^i=R_{00}^{\prime }+R_{km}^{\prime }{U}^{\prime k}{U}^{{}^{\prime }m}+\frac{\partial f^{{}^{\prime }j}}{\partial x^{\prime j}}$$

(92)

With the use of the gravitational field equations Equation (63), we can express K as:

$$\begin{array}{ccc}\hfill K& =& \Lambda \left(R^{\prime },L_m^{\prime }\right){\eta }_{00}+\frac{1}{f_R\left(R^{\prime },L_m\right)}\frac{{\partial }^2}{\partial t^{\prime 2}}{f}_R\left(R^{\prime },L_m^{\prime }\right)+\Phi \left(R^{\prime },L_m^{\prime }\right)\left(T_{00}^{\prime }-\frac{1}{3}{T}^{\prime }{\eta }_{00}\right)\hfill \end{array}$$

$$\begin{array}{ccc}& & +\Lambda \left(R^{\prime },L_m^{\prime }\right){\eta }_{km}{U}^{\prime k}{U}^{{}^{\prime }m}+\frac{1}{f_R\left(R^{\prime },L_m^{\prime }\right)}{U}^{\prime k}{U}^{\prime m}\frac{{\partial }^2}{\partial x^{\prime k}\partial x^{\prime m}}{f}_R\left(R^{\prime },L_m^{\prime }\right)\hfill \end{array}$$

$$\begin{array}{ccc}& & +\varphi \left(R^{\prime },L_m^{\prime }\right)\left(T_{km}^{\prime }-\frac{1}{3}T{\eta }_{km}\right)U^{\prime k}{U}^{\prime m}\hfill \end{array}$$

$$\begin{array}{ccc}& & +\frac{\partial }{\partial x^{\prime k}}\left\{\frac{\partial }{\partial x^{\prime m}}ln\left[f_{L_m}\left(R^{\prime },L_m^{\prime }\right)\frac{d{L}_m^{\prime }\left(\rho \right)}{d\rho }\right]\left(L_m^{\prime }{\eta }^{km}-T^{\prime km}\right)\right\}\hfill \end{array}$$

(93)

Then, since in the Newtonian limit, one can omit all time derivatives, we obtain for $R_{0l0}^i$ the expression [111]:

$$R_{0l0}^i=\frac{{\partial }^2\varphi }{\partial x^i\partial x^l}$$

(94)

For simplicity, in the following, we omit the primes for the geometrical and physical quantities in the Newtonian approximation. Therefore, in $f\left(R,L_m\right)$-modified gravity with a curvature-matter coupling, we obtain for the tidal acceleration of the test particles the expression:

$$\frac{d^2{\eta }^i}{d{t}^2}=F^i\approx \frac{{\partial }^2\varphi }{\partial x^i\partial x^l}{\eta }^l+R_{jlm}^i{\eta }^l{V}^j{V}^m+{\eta }^l\frac{\partial f^i}{\partial x^l}$$

(95)

where $V^j$ and $V^m$ are the Newtonian three-dimensional velocities. In the Newtonian approximation, in modified theories of gravity with a curvature-matter coupling, the tidal force tensor is defined as:

$$\frac{\partial F^i}{\partial {\eta }^l}=\frac{{\partial }^2\varphi }{\partial x^i\partial x^l}+R_{jlm}^i{V}^j{V}^m+\frac{\partial f^i}{\partial x^l}$$

(96)

and its trace gives the generalized Poisson equation:

$$\frac{\partial F^i}{\partial {\eta }^i}=\Delta \varphi +R_{jm}{V}^j{V}^m+\frac{\partial f^i}{\partial x^i}$$

(97)

In Newtonian gravity:

$${\chi }_{il}=-\frac{{\partial }^2\varphi }{\partial x^i\partial x^l}$$

(98)

represents the Newtonian tidal tensor [111]. In the Newtonian approximation, the spherical potential of a given particle with mass M is $\varphi \left(r\right)=-M/8\pi r$. By choosing a frame of reference, so that the x-axis passes through the particle’s position, corresponding to the radial spherical coordinate, that is, $(x=r,y=0,z=0)$, the Newtonian tidal tensor is diagonal and has the only non-zero components:

$${\chi }_{ii}=\mathrm{diag}\left(\frac{2M}{8\pi r^3},-\frac{M}{8\pi r^3},-\frac{M}{8\pi r^3}\right)$$

(99)

The Newtonian tidal force ${\overrightarrow{F}}_t$ can be written as $F_{tx}=2M\Delta x/8\pi r^3$, $F_{ty}=-GM\Delta y/8\pi r^3$ and $F_{tz}=-M\Delta z/8\pi r^3$, respectively [111]. These results can be used to derive the generalization of the Roche limit in modified gravity with an arbitrary coupling between matter and curvature.

The Roche limit is the closest distance $r_{Roche}$ that a celestial object with mass m, radius $R_m$ and density ${\rho }_m$, held together only by its own gravity, can come to a massive body of mass M, radius $R_M$ and density ${\rho }_M$, respectively, without being pulled apart by the massive object’s tidal (gravitational) force [108,109,110]. For simplicity, we will consider $M\gg m$, so that the center of mass of the system nearly coincides with the geometrical center of the mass M.

The elementary Newtonian theory of this process is as follows. Consider a small mass $\Delta m$ located at the surface of the small object of mass m. There are two forces acting on $\Delta m$, the gravitational attraction of the mass m, given by:

$$F_G=\frac{m\Delta m}{8\pi R_m^2}$$

(100)

and the tidal force exerted by the massive object, which is given by:

$$F_{tr}=\frac{M\Delta mR}{8\pi r^3}$$

(101)

where r is the distance between the centers of the two celestial bodies. The Roche limit is reached at the distance $r=r_{Roche}$, when the gravitational force and the tidal force exactly balance each other, $F_G=F_{tr}$, thus giving [108,109,110]:

$$r_{Roche}=R_m{\left(\frac{M}{m}\right)}^{1/3}=2^{1/3}{R}_M{\left(\frac{{\rho }_M}{{\rho }_m}\right)}^{1/3}$$

(102)

In modified gravity with a curvature-matter coupling, the equilibrium between gravitational and tidal forces occurs at a distance $r_{Roche}$ given by the equation:

$$\left(\frac{M}{8\pi r_{Roche}^3}+R_{jrm}^r{V}^j{V}^m+\frac{\partial f^r}{\partial r}\right)R_m=\frac{m}{8\pi R_m^2}+f^r$$

(103)

where $f^r$ is the radial component of the extra-force, which modifies the Newtonian gravitational force, and the curvature tensor $R_{jrm}^r$ (no summation upon r) must be evaluated in the coordinate system in which the Newtonian tidal tensor is diagonal. Hence, we obtain the generalized Roche limit in the presence of arbitrary curvature-matter coupling as [100]:

$$\begin{array}{ccc}\hfill r_{Roche}& \approx & R_m{\left(\frac{M}{m}\right)}^{1/3}\left[1+\frac{8\pi R_m^3}{3m}\left(R_{jrm}^r{V}^j{V}^m+\frac{\partial f^r}{\partial r}\right)-\frac{8\pi R_m^2}{3m}{f}^r\right]\hfill \end{array}$$

(104)

where we have assumed that the gravitational effects due to the coupling between matter and curvature are small as compared to the Newtonian ones. Relativistic corrections to the Newtonian tidal accelerations caused by a massive rotating source, such as, for example, the Earth, could be determined experimentally, at least in principle, thus leading to the possibility of testing relativistic theories of gravitation by measuring such effects in a laboratory.

4. Extended $\mathbf{f}\left(\mathbf{R},{\mathbf{L}}_{\mathbf{m}}\right)$ Gravity with the Generalized Scalar Field and Kinetic Term Dependencies

4.1. Action and Field Equations

We generalize $f\left(R,L_m\right)$ gravity, outlined in the previous section, by considering a novel gravitational model with an action given by an arbitrary function of the Ricci scalar, the matter Lagrangian density, a scalar field and a kinetic term constructed from the gradients of the scalar field, respectively [24]. The action is given by:

$$S=\int f\left(R,L_m,\varphi ,g^{\mu \nu }{\nabla }_{\mu }\varphi {\nabla }_{\nu }\varphi \right)\sqrt{-g}{d}^{4}x$$

(105)

where $\sqrt{-g}$ is the determinant of the metric tensor $g_{\mu \nu }$, and $f\left(R,L_m,\varphi ,g^{\mu \nu }{\nabla }_{\mu }\varphi {\nabla }_{\nu }\varphi \right)$ is an arbitrary function of the Ricci scalar R, the matter Lagrangian density, $L_m$, a scalar field ϕ and the gradients constructed from the scalar field, respectively. The only restriction on the function f is to be an analytical function of R, $L_m$, ϕ, and of the scalar field kinetic energy, respectively, that is, f must possess a Taylor series expansion about any point.

Note that one may motivate the introduction of the action Equation (105) by the extensive interest in the literature between couplings between the scalar and the matter fields [112,113,114,115,116,117,118,119,120,121,122,123,124,125,126]. Indeed, such couplings generically appear in Kaluza–Klein theories with compactified dimensions [127] or in the low energy effective limit of string theories [112,113,115,116,117], where the dilaton has been proposed to be a good candidate for quintessence [112] or the inflaton [113]. The theory represented by action Equation (105) offers a generalization of the above theories in a single theoretical framework.

Here, we introduce first the “reduced” energy-momentum tensor ${\tau }_{\mu \nu }$ of the matter, which is defined as [67],

$${\tau }_{\mu \nu }=-\frac{2}{\sqrt{-g}}\frac{\delta \left(\sqrt{-g}{L}_m\right)}{\delta g^{\mu \nu }}=g_{\mu \nu }{L}_m-2\frac{\delta L_m}{\delta g^{\mu \nu }}$$

(106)

In addition to this, we assume that the scalar field ϕ is independent of the metric, i.e., $\delta \varphi /\delta g^{\mu \nu }\equiv 0$, and in the following, we define, for simplicity, ${\left(\nabla \varphi \right)}^2=g^{\mu \nu }{\nabla }_{\mu }\varphi {\nabla }_{\nu }\varphi$.

By varying the action S with respect to the metric tensor and to ϕ, we obtain the following field equation:

$$\begin{array}{c}\hfill f_R{R}_{\mu \nu }+\left(g_{\mu \nu }{\nabla }_{\lambda }{\nabla }^{\lambda }-{\nabla }_{\mu }{\nabla }_{\nu }\right)f_R-\frac{1}{2}\left(f-f_{L_m}{L}_m\right)g_{\mu \nu }=\frac{1}{2}{f}_{L_m}{\tau }_{\mu \nu }-f_{{\left(\nabla \varphi \right)}^2}{\nabla }_{\mu }\varphi {\nabla }_{\nu }\varphi \end{array}$$

(107)

and the evolution equation for the scalar field:

$${\square }_{{\left(\nabla \varphi \right)}^2}\varphi =\frac{1}{2}{f}_{\varphi }$$

(108)

respectively. The subscript of f denotes a partial derivative with respect to the arguments, i.e., $f_R=\partial f/\partial R$, $f_{L_m}=\partial f/\partial L_m$, $f_{{\left(\nabla \varphi \right)}^2}=\partial f/\partial {\left(\nabla \varphi \right)}^2$, $f_{\varphi }=\partial f/\partial \varphi$, and:

$${\square }_{{\left(\nabla \varphi \right)}^2}=\frac{1}{\sqrt{-g}}\frac{\partial }{\partial x^{\mu }}\left[f_{{\left(\nabla \varphi \right)}^2}\sqrt{-g}{g}^{\mu \nu }\frac{\partial }{\partial x^{\nu }}\right]$$

(109)

is the generalized d’Alembert operator of $f\left(R,L_m,\varphi ,{\nabla }_{\mu }\varphi {\nabla }^{\mu }\varphi \right)$ gravity.

The contraction of Equation (107) provides the following relation between the Ricci scalar R, the matter Lagrangian density $L_m$, the derivatives of the scalar field and the trace $\tau ={\tau }_{\mu }^{\mu }$ of the “reduced” energy-momentum tensor:

$$\begin{array}{c}\hfill f_{R}R+3{\nabla }_{\mu }{\nabla }^{\mu }{f}_R-2\left(f-f_{L_m}{L}_m\right)=\frac{1}{2}{f}_{L_m}\tau -f_{{\left(\nabla \varphi \right)}^2}{\nabla }_{\mu }\varphi {\nabla }^{\mu }\varphi \end{array}$$

(110)

By taking the covariant divergence of Equation (107), we obtain for the covariant divergence of the “reduced” energy-momentum tensor the following expression:

$$\begin{array}{c}\hfill \frac{1}{2}{\nabla }^{\sigma }\left(f_{L_m}{\tau }_{\mu \sigma }\right)=\frac{1}{2}\left(L_m{\nabla }_{\mu }{f}_{L_m}-f_{\varphi }{\nabla }_{\mu }\varphi \right)+f_{{\left(\nabla \varphi \right)}^2}{\nabla }_{\mu }\varphi {\nabla }_{\sigma }{\nabla }^{\sigma }\varphi +{\nabla }_{\mu }\varphi {\nabla }^{\sigma }\varphi {\nabla }_{\sigma }{f}_{{\left(\nabla \varphi \right)}^2}\end{array}$$

(111)

This relationship was deduced by taking into account the following mathematical identities:

$$\begin{array}{c}\hfill {\nabla }^{\mu }{R}_{\mu \nu }=\frac{1}{2}{\nabla }_{\nu }R,\left({\nabla }_{\nu }\square -\square {\nabla }_{\nu }\right)f_R=-\left({\nabla }^{\mu }{f}_R\right)R_{\mu \nu }\end{array}$$

(112)

and considering torsion-free space-times, such that $\left[{\nabla }_{\sigma }{\nabla }_{ϵ}-{\nabla }_{ϵ}{\nabla }_{\sigma }\right]\psi =0$, where ψ is any scalar-field. Now, using Equation (109), we get:

$$\begin{array}{c}\hfill {\nabla }^{\sigma }\left(f_{L_m}{\tau }_{\mu \sigma }\right)=L_m{\nabla }_{\mu }{f}_{L_m}\end{array}$$

(113)

For $\varphi \equiv 0$, Equation (107) reduces to the field equations of the $f\left(R,L_m\right)$ model considered in the previous section [15]. For $\varphi \ne 0$, one recovers the conservation equations for either General Relativity and Brans–Dicke-like scalar-tensor theories (with and without scalar/matter coupling [128]). For instance, the total Lagrangian of the simplest matter-scalar field-gravitational field theory, with a scalar field kinetic term and a self-interacting potential $V\left(\varphi \right)$, corresponds to the choice:

$$f=\frac{R}{2}+L_m+\frac{\lambda }{2}{g}^{\mu \nu }{\nabla }_{\mu }\varphi {\nabla }_{\nu }\varphi +V\left(\varphi \right)$$

(114)

where λ is a constant. The corresponding field equations can be immediately obtained from Equation (107) as:

$$\begin{array}{c}\hfill R_{\mu \nu }-\frac{1}{2}R{g}_{\mu \nu }={\tau }_{\mu \nu }-\lambda {\nabla }_{\mu }\varphi {\nabla }_{\nu }\varphi +\left[\frac{\lambda }{2}{g}^{\alpha \beta }{\nabla }_{\alpha }\varphi {\nabla }_{\beta }\varphi +V\left(\varphi \right)\right]g_{\mu \nu }\end{array}$$

(115)

The scalar field satisfies the evolution equation:

$$\frac{1}{\sqrt{-g}}\frac{\partial }{\partial x^{\mu }}\left[\sqrt{-g}{g}^{\mu \nu }\frac{\partial \varphi }{\partial x^{\nu }}\right]=\frac{1}{\lambda }\frac{dV\left(\varphi \right)}{d\varphi }$$

(116)

while the energy-momentum tensor obeys the conservation equation ${\nabla }^{\sigma }{\tau }_{\mu \sigma }=0$.

4.2. Models with Nonminimal Matter-Scalar Field Coupling

As an example of the application of the formalism developed in the previous section, we consider a simple phenomenological model, in which a scalar field is non-minimally coupled to pressureless matter with rest mass density ρ. For the action of the system, we consider:

$$S=\int \left[\frac{R}{2}-F\left(\varphi \right)\rho +\lambda g^{\mu \nu }{\nabla }_{\mu }\varphi {\nabla }_{\nu }\varphi \right]\sqrt{-g}{d}^{4}x$$

(117)

where $F\left(\varphi \right)$ is an arbitrary function of the scalar field that couples non-minimally to ordinary matter. The field equations for this model are given by:

$$\begin{array}{ccc}\hfill R_{\mu \nu }-\frac{1}{2}R{g}_{\mu \nu }& =& F\left(\varphi \right)\rho U_{\mu }{U}_{\nu }+2\lambda \left[{\nabla }_{\mu }\varphi {\nabla }_{\nu }\varphi -\frac{1}{2}{g}_{\mu \nu }{\nabla }_{\alpha }\varphi {\nabla }^{\alpha }\varphi \right]\hfill \end{array}$$

(118)

where $U^{\mu }$ is the four-velocity of the matter fluid. The scalar field satisfies the evolution equation:

$$\square \varphi =-\frac{1}{2\lambda }\frac{dF\left(\varphi \right)}{d\varphi }\rho$$

(119)

where □ is the usual d’Alembert operator defined in a curved space. The total energy-momentum tensor of the scalar field-matter system is given by:

$$T_{\mu \nu }=F\left(\varphi \right)\rho U_{\mu }{U}_{\nu }+2\lambda \left[{\nabla }_{\mu }\varphi {\nabla }_{\nu }\varphi -\frac{1}{2}{g}_{\mu \nu }{\nabla }_{\alpha }\varphi {\nabla }^{\alpha }\varphi \right]$$

(120)

The Bianchi identities imply that ${\nabla }_{\nu }{T}^{\mu \nu }=0$, and in the following, we assume that the mass density current is conserved, i.e., ${\nabla }_{\nu }\left(\rho U^{\nu }\right)=0$. Using the latter condition and the mathematical identity given by $\left[{\nabla }_{\sigma }{\nabla }_{ϵ}-{\nabla }_{ϵ}{\nabla }_{\sigma }\right]\psi =0$, we have:

$$F\left(\varphi \right)\rho \left(U^{\nu }{\nabla }_{\nu }{U}^{\mu }\right)+\rho U^{\mu }{U}^{\nu }\frac{dF\left(\varphi \right)}{d\varphi }{\nabla }_{\nu }\varphi +2\lambda \left({\nabla }^{\mu }\varphi \right)\square \varphi =0$$

(121)

Then, by eliminating the term $\square \varphi$ with the help of Equation (119), we obtain:

$$U^{\nu }{\nabla }_{\nu }{U}^{\mu }+\left[\frac{d}{d\varphi }lnF\left(\varphi \right)\right]\left(U^{\mu }{U}^{\nu }{\nabla }_{\nu }\varphi -{\nabla }^{\mu }\varphi \right)=0$$

(122)

Using the identity $U^{\nu }{\nabla }_{\nu }{U}^{\mu }\equiv \frac{d^2{x}^{\mu }}{d{s}^2}+{\Gamma }_{\alpha \beta }^{\mu }{u}^{\alpha }{u}^{\beta },$ where ${\Gamma }_{\alpha \beta }^{\mu }$ are the Christoffel symbols corresponding to the metric, the equation of motion of the test particles non-minimally coupled to an arbitrary scalar field takes the form:

$$\frac{d^2{x}^{\mu }}{d{s}^2}+{\Gamma }_{\alpha \beta }^{\mu }{U}^{\alpha }{U}^{\beta }+\left[\frac{d}{d\varphi }lnF\left(\varphi \right)\right]\left(U^{\mu }{U}^{\nu }{\nabla }_{\nu }\varphi -{\nabla }^{\mu }\varphi \right)=0$$

(123)

A particular model can be obtained by assuming that $F\left(\varphi \right)$ is given by a linear function:

$$F\left(\varphi \right)=\frac{\Lambda +1}{2}\left[1+\frac{1}{2}\left(\Lambda -1\right)\varphi \right]$$

(124)

where Λ is a constant. Then, the equation of motion becomes:

$$\frac{d^2{x}^{\mu }}{d{s}^2}+{\Gamma }_{\alpha \beta }^{\mu }{U}^{\alpha }{U}^{\beta }+\left(U^{\mu }{U}^{\nu }-g^{\mu \nu }\right){\nabla }_{\nu }ln\left[1+\frac{\Lambda -1}{2}\varphi \right]=0$$

(125)

In order to simplify the field equations, we adopt for λ the value $\lambda =-\left({\Lambda }^2-1\right)/8$. Then, Equation (119), determining the scalar field, takes the simple form $\square \varphi =\rho$.

The gravitational field equations take the form:

$$R_{\mu \nu }-\frac{1}{2}{g}_{\mu \nu }R=\frac{\Lambda +1}{2}{T}_{\mu \nu }$$

(126)

with the total energy-momentum tensor given by:

$$\begin{array}{c}\hfill T_{\mu \nu }=\left[1+\frac{\Lambda -1}{2}\varphi \right]\rho U_{\mu }{U}_{\nu }-\frac{\Lambda -1}{2}\left[{\nabla }_{\mu }\varphi {\nabla }_{\nu }\varphi -\frac{1}{2}{g}_{\mu \nu }{\nabla }_{\alpha }\varphi {\nabla }^{\alpha }\varphi \right]\end{array}$$

(127)

For $\Lambda =1$, we reobtain the general relativistic model for dust. Other possible choices of the function $F\left(\varphi \right)$, such as $F\left(\varphi \right)=exp\left(\varphi \right)$, can be discussed in a similar way.

A more general model can be obtained by adopting for the matter Lagrangian the general expression [66,68,128]:

$$L_m=-\left[\rho +\rho \int \frac{dp\left(\rho \right)}{\rho }-p\left(\rho \right)\right]$$

(128)

where ρ is the rest-mass energy density and p is the thermodynamic pressure, which, by assumption, satisfies a barotropic equation of state, $p=p\left(\rho \right)$. By assuming that the matter Lagrangian does not depend on the derivatives of the metric and that the particle matter fluid current is conserved (${\nabla }_{\nu }\left(\rho u^{\nu }\right)=0$), the Lagrangian given by Equation (128) is the unique matter Lagrangian that can be constructed from the thermodynamic parameters of the fluid, and it is valid for all gravitational theories satisfying the two previously mentioned conditions [128].

The gravitational field equations and the equation describing the matter-scalar field coupling are given by:

$$R_{\mu \nu }-\frac{1}{2}{g}_{\mu \nu }R=F\left(\varphi \right)ϵU_{\mu }{U}_{\nu }-p{g}_{\mu \nu }+\lambda Q_{\mu \nu }$$

(129)

$$\square \varphi =\frac{1}{2\lambda }\frac{dF\left(\varphi \right)}{d\varphi }ϵ$$

(130)

where $Q_{\mu \nu }={\nabla }_{\mu }\varphi {\nabla }_{\nu }\varphi -\frac{1}{2}{\nabla }_{\lambda }\varphi {\nabla }^{\lambda }\varphi g_{\mu \nu }$ and where the total energy density is $ϵ=\rho +\rho \int dp/\rho -p$ [68,128]. With the use of the conservation equation ${\nabla }_{\nu }\left(\rho U^{\nu }\right)=0$, one obtains the equation of motion of massive test particles as:

$$\frac{d^2{x}^{\mu }}{d{s}^2}+{\Gamma }_{\alpha \beta }^{\mu }{U}^{\alpha }{U}^{\beta }+\left(U^{\mu }{U}^{\nu }-g^{\mu \nu }\right){\nabla }_{\nu }ln\left[1+\int \frac{dp}{\rho }\right]=0$$

(131)

Models with scalar field-matter coupling were considered in the framework of the Brans–Dicke theory [124], with the action of the model given by $S=\int \left[\varphi R/2+\left(\omega /\varphi \right){\nabla }_{\mu }\varphi {\nabla }^{\mu }\varphi +F\left(\varphi \right)L_m\right]\sqrt{-g}{d}^{4}x$. Such models can give rise to a late-time accelerated expansion of the Universe for very high values of the Brans–Dicke parameter ω. Other models with interacting scalar field and matter have been considered in [120]. We emphasize that the gravitational theory considered in this work generalizes all of the above models.

In conclusion, the general formalism outlined in this work can be extremely useful in a variety of scenarios, such as, in describing the interaction between dark energy, modeled as a scalar field, and dark matter, or ordinary matter (neutrinos), with or without pressure, matter-scalar field interactions in inflation, as well as in the study of the interactions of the scalar field (representing dark matter and/or dark energy) and the electromagnetic component in the very early Universe. Moreover, they can provide a realistic description of the late expansion of the Universe, where a possible interaction between ordinary matter and dark energy cannot be excluded a priori.

5. $\mathbf{f}(\mathbf{R},\mathbf{T})$ Gravity

5.1. Action and Gravitational Field Equations

In this section, we consider another extension of standard General Relativity, namely, $f(R,T)$-modified theories of gravity, where the gravitational Lagrangian is given by an arbitrary function of the Ricci scalar R and of the trace of the energy-momentum tensor T [16]. The dependence from T may be induced by exotic imperfect fluids or quantum effects (conformal anomaly). The action takes the following form:

$$S=\frac{1}{16\pi }\int f\left(R,T\right)\sqrt{-g}{d}^{4}x+\int L_m\sqrt{-g}{d}^{4}x$$

(132)

where $f\left(R,T\right)$ is an arbitrary function of the Ricci scalar, R, and of the trace T of the energy-momentum tensor of the matter, $T_{\mu \nu }$. For the matter content, we assume that it consists of a fluid that can be characterized by two thermodynamic parameters only, the energy density and the pressure, respectively. With the use of Equation (17), it follows that the trace of the energy-momentum tensor T can be expressed as a function of the matter Lagrangian as:

$$T=g^{\mu \nu }{T}_{\mu \nu }=4L_m-2g^{\mu \nu }\frac{\delta L_m}{\delta g^{\mu \nu }}$$

(133)

Hence, $f(R,T)$ gravity theories can be interpreted as extensions of the $f\left(R,L_m\right)$-type gravity theories, with the gravitational action depending not only on the matter Lagrangian only, but also on its variation with respect to the metric.

By varying the action S of the gravitational field with respect to the metric tensor components $g^{\mu \nu }$, we obtain the field equations of the $f\left(R,T\right)$ gravity model as:

$$\begin{array}{c}\hfill f_R\left(R,T\right)R_{\mu \nu }-\frac{1}{2}f\left(R,T\right)g_{\mu \nu }+\left(g_{\mu \nu }\square -{\nabla }_{\mu }{\nabla }_{\nu }\right)f_R\left(R,T\right)\\ \hfill =8\pi T_{\mu \nu }-f_T\left(R,T\right)T_{\mu \nu }-f_T\left(R,T\right){\Theta }_{\mu \nu }\end{array}$$

(134)

We have defined the variation of T with respect to the metric tensor as:

$$\frac{\delta \left(g^{\alpha \beta }{T}_{\alpha \beta }\right)}{\delta g^{\mu \nu }}=T_{\mu \nu }+{\Theta }_{\mu \nu }$$

(135)

where:

$${\Theta }_{\mu \nu }\equiv g^{\alpha \beta }\frac{\delta T_{\alpha \beta }}{\delta g^{\mu \nu }}$$

(136)

Note that when $f(R,T)\equiv f\left(R\right)$, from Equation (134), we obtain the field equations of $f\left(R\right)$ gravity.

Contracting Equation (134) gives the following relation between the Ricci scalar R and the trace T of the energy-momentum tensor:

$$\begin{array}{ccc}& & f_R\left(R,T\right)R+3\square f_R\left(R,T\right)-2f\left(R,T\right)=8\pi T-f_T\left(R,T\right)T-f_T\left(R,T\right)\Theta \hfill \end{array}$$

(137)

where we have denoted $\Theta ={\Theta }_{\mu }^{\mu }$.

By eliminating the term $\square f_R\left(R,T\right)$ between Equation (134) and Equation (137), the gravitational field equations can be written in the form:

$$\begin{array}{ccc}\hfill f_R\left(R,T\right)\left(R_{\mu \nu }-\frac{1}{3}R{g}_{\mu \nu }\right)+\frac{1}{6}f\left(R,T\right)g_{\mu \nu }& =& 8\pi \left(T_{\mu \nu }-\frac{1}{3}T{g}_{\mu \nu }\right)-f_T\left(R,T\right)\left(T_{\mu \nu }-\frac{1}{3}T{g}_{\mu \nu }\right)\hfill \\ & & -f_T\left(R,T\right)\left({\Theta }_{\mu \nu }-\frac{1}{3}\Theta g_{\mu \nu }\right)+{\nabla }_{\mu }{\nabla }_{\nu }{f}_R\left(R,T\right)\hfill \end{array}$$

(138)

Taking into account the covariant divergence of Equation (134), with the use of the following mathematical identity [50]:

$$\begin{array}{c}\hfill {\nabla }^{\mu }\left[f_R\left(R,T\right)R_{\mu \nu }-\frac{1}{2}f\left(R,T\right)g_{\mu \nu }+\left(g_{\mu \nu }\square -{\nabla }_{\mu }{\nabla }_{\nu }\right)f_R\left(R,T\right)\right]\equiv -\frac{1}{2}{g}_{\mu \nu }{f}_T(R,T){\nabla }^{\mu }T\end{array}$$

(139)

we obtain, for the divergence of the energy-momentum tensor $T_{\mu \nu }$, the equation:

$${\nabla }^{\mu }{T}_{\mu \nu }=\frac{f_T\left(R,T\right)}{8\pi -f_T\left(R,T\right)}\left[\left(T_{\mu \nu }+{\Theta }_{\mu \nu }\right){\nabla }^{\mu }ln{f}_T\left(R,T\right)+{\nabla }^{\mu }{\Theta }_{\mu \nu }-\frac{1}{2}{g}_{\mu \nu }{\nabla }^{\mu }T\right]$$

(140)

Once the matter Lagrangian is known, the calculation of the tensor ${\Theta }_{\mu \nu }$ yields (we refer the reader to [16] for more details):

$${\Theta }_{\mu \nu }=-2T_{\mu \nu }+g_{\mu \nu }{L}_m-2g^{\alpha \beta }\frac{{\partial }^2{L}_m}{\partial g^{\mu \nu }\partial g^{\alpha \beta }}$$

(141)

More specifically, in the case of the electromagnetic field, the matter Lagrangian is given by:

$$L_m=-\frac{1}{16\pi }{F}_{\alpha \beta }{F}_{\gamma \sigma }{g}^{\alpha \gamma }{g}^{\beta \sigma }$$

(142)

where $F_{\alpha \beta }$ is the electromagnetic field tensor, we obtain ${\Theta }_{\mu \nu }=-T_{\mu \nu }$. For a massless scalar field ϕ with Lagrangian $L_m=g^{\alpha \beta }{\nabla }_{\alpha }\varphi {\nabla }_{\beta }\varphi$, we obtain ${\Theta }_{\mu \nu }=-T_{\mu \nu }+(1/2)T{g}_{\mu \nu }$.

The problem of the perfect fluids, described by an energy density ρ, pressure p and four-velocity $U^{\mu }$, is more subtle, since there is no unique definition of the matter Lagrangian. However, in the present study we assume that the energy-momentum tensor of the matter is given by $T_{\mu \nu }=\left(\rho +p\right)U_{\mu }{U}_{\nu }+p{g}_{\mu \nu }$, and the matter Lagrangian can be taken as $L_m=p$. The four-velocity $U_{\mu }$ satisfies the conditions $U_{\mu }{U}^{\mu }=-1$ and $U^{\mu }{\nabla }_{\nu }{U}_{\mu }=0$, respectively. Then, with the use of Equation (141), we obtain, for the variation of the energy-momentum of a perfect fluid, the expression:

$${\Theta }_{\mu \nu }=-2T_{\mu \nu }+p{g}_{\mu \nu }$$

(143)

The choice $L_m=-\rho$ leads to a different form for ${\Theta }_{\mu \nu }$, ${\Theta }_{\mu \nu }=-2T_{\mu \nu }+\rho g_{\mu \nu }$. With this form of the matter Lagrangian, the gravitational field equations take a different form. The nature of the qualitative and quantitative differences between the two different versions of the $f(R,T)$ gravity, corresponding to different matter Lagrangians, is still an open question.

5.2. Specific Cosmological Solution

In the present Section, we consider a particular class of $f(R,T)$-modified gravity models, obtained by explicitly specifying the functional form of f. Generally, the field equations also depend, through the tensor ${\Theta }_{\mu \nu }$, on the physical nature of the matter field. Hence, in the case of $f(R,T)$ gravity, depending on the nature of the matter source, for each choice of f, we can obtain several theoretical models, corresponding to different matter models [129].

As a specific case of a $f(R,T)$ model, consider a correction to the Einstein–Hilbert action given by $f\left(R,T\right)=R+2f\left(T\right)$, where $f\left(T\right)$ is an arbitrary function of the trace of the energy-momentum tensor of matter. The gravitational field equations immediately follow from Equation (134), and is given by:

$$R_{\mu \nu }-\frac{1}{2}R{g}_{\mu \nu }=8\pi T_{\mu \nu }-2f^{\prime }\left(T\right)T_{\mu \nu }-2f^{\prime }\left(T\right){\Theta }_{\mu \nu }+f\left(T\right)g_{\mu \nu }$$

(144)

where the prime denotes a derivative with respect to the argument.

Consider a perfect fluid, so that ${\Theta }_{\mu \nu }=-2T_{\mu \nu }-p{g}_{\mu \nu }$, and the field equations become:

$$R_{\mu \nu }-\frac{1}{2}R{g}_{\mu \nu }=8\pi T_{\mu \nu }+2f^{\prime }\left(T\right)T_{\mu \nu }+\left[2p{f}^{\prime }\left(T\right)+f\left(T\right)\right]g_{\mu \nu }$$

(145)

In the case of dust with $p=0$, the gravitational field equations are given by:

$$R_{\mu \nu }-\frac{1}{2}R{g}_{\mu \nu }=8\pi T_{\mu \nu }+2f^{\prime }\left(T\right)T_{\mu \nu }+f\left(T\right)g_{\mu \nu }$$

(146)

These field equations were proposed in [130] to solve the cosmological constant problem. The simplest cosmological model can be obtained by assuming a dust universe ($p=0$, $T=-\rho$) and by choosing the function $f\left(T\right)$, so that $f\left(T\right)=\lambda T$, where λ is a constant. By assuming that the metric of the universe is given by the flat Robertson–Walker metric:

$$d{s}^2=-d{t}^2+a^2\left(t\right)\left(d{x}^2+d{y}^2+d{z}^2\right)$$

(147)

the gravitational field equations are given by:

$$\begin{array}{ccc}\hfill 3\frac{{\dot{a}}^2}{a^2}& =& \left(8\pi +3\lambda \right)\rho \hfill \end{array}$$

(148)

$$\begin{array}{ccc}\hfill 2\frac{\ddot{a}}{a}+\frac{{\dot{a}}^2}{a^2}& =& \lambda \rho \hfill \end{array}$$

(149)

respectively. Thus, this $f(R,T)$ gravity model is equivalent to a cosmological model with an effective dark energy density ${\Lambda }_{\mathrm{eff}}\propto H^2$, where $H=\dot{a}/a$ is the Hubble function [130]. It is also interesting to note that, generally, for this choice of $f(R,T)$, the gravitational coupling becomes an effective and time-dependent coupling, of the form $G_{\mathrm{eff}}=G\pm 2f^{\prime }\left(T\right)$. Thus, the term $2f\left(T\right)$ in the gravitational action modifies the gravitational interaction between matter and curvature, replacing G by a running gravitational coupling parameter.

The field equations reduce to a single equation for H:

$$2\dot{H}+3\frac{8\pi +2\lambda }{8\pi +3\lambda }{H}^2=0$$

(150)

with the general solution given by:

$$H\left(t\right)=\frac{2\left(8\pi +3\lambda \right)}{3\left(8\pi +2\lambda \right)}\frac{1}{t}$$

(151)

The scale factor evolves according to $a\left(t\right)=t^{\alpha }$, with $\alpha =2\left(8\pi +3\lambda \right)/3\left(8\pi +2\lambda \right)$.

5.3. Further Applications

The $f(R,T)$ gravitational model analysed above has been given a great amount of recent attention. We will briefly outline a few cosmological applications. In [131], the cosmological reconstruction of $f(R,T)$ gravity describing matter-dominated and accelerated phases was analysed. Special attention was paid to the specific case of $f(R,T)=f_1\left(R\right)+f_2\left(T\right)$. The use of an auxiliary scalar field was considered with two known examples for the scale factor corresponding to an expanding Universe. In the first example, where ordinary matter is usually neglected for obtaining the unification of matter-dominated and accelerated phases with $f\left(R\right)$ gravity, it was shown that this unification can be obtained in the presence of ordinary matter. In the second example, as in $f\left(R\right)$ gravity, the model of $f(R,T)$ gravity with a transition of the matter-dominated phase to the acceleration phase was obtained.

In [132], it was shown that the dust fluid reproduces ΛCDM, he phantom-non-phantom era and phantom cosmology. Furthermore, different cosmological models were reconstructed, including the Chaplygin gas and the scalar field with some specific forms of $f(R,T)$. The numerical simulation for the Hubble parameter shows good agreement with the baryon acoustic oscillations (BAO) observational data for low redshifts $z<2$. In [133], through a numerical reconstruction, it was shown that specific $f(R,T)$ models are able to reproduce the same expansion history generated, in standard General Relativity, by dark matter and holographic dark energy. It was further shown that these theories are able to reproduce the four known types of future finite-time singularities [134]. In [135], a non-equilibrium picture of thermodynamics was discussed at the apparent horizon of an FRW universe in $f(R,T)$ gravity, where the validity of the first and second law of thermodynamics in this scenario were checked. It was shown that the Friedmann equations can be expressed in the form of the first law of thermodynamics and that the second law of thermodynamics holds both in the phantom and non-phantom phases.

The energy conditions have also been extensively explored in $f(R,T)$ gravity, for instance, in an FRW universe with perfect fluid [136], and in the context of exact power-law solutions [137], where certain constraints were found that have to be satisfied to ensure that power law solutions may be stable and match the bounds prescribed by the energy conditions. In [138], it was also shown that the energy conditions were satisfied for specific models. Furthermore, an analysis of the perturbations and stabilities of de Sitter solutions and power-law solutions was performed, and it was shown that for those models in which the energy conditions are satisfied, de Sitter solutions and power-law solutions may be stable.

Solutions of the Gödel universe in the framework of $f(R,T)$-modified theories of gravity were also obtained [139]. Algorithms were derived for constructing five-dimensional Kaluza–Klein cosmological space-times in the presence of a perfect fluid source in $f(R,T)$ gravity [140]. Spatially homogeneous and anisotropic Bianchi type-V cosmological models in a scalar-tensor theory of gravitation [141] were also explored. The Bianchi type I model with perfect fluid as a matter content in $f(R,T)$ gravity was analysed, and the physical and kinematical properties of the model were also discussed [142]. In [143], the specific case of the conservation of the energy-momentum tensor was considered, and cosmological solutions were obtained for a homogeneous and isotropic model of the Universe.

In [144], the evolution of scalar cosmological perturbations in the metric formalism were analysed. According to restrictions on the background evolution, a specific model within these theories was assumed in order to guarantee the standard continuity equation. Using a completely general procedure, the complete set of differential equations for the matter density perturbations was found. In the case of sub-Hubble modes, the density contrast evolution reduces to a second-order equation, and it was shown that for well-motivated $f(R,T)$ Lagrangians, the quasistatic approximation provide different results from the ones derived in the frame of the concordance ΛCDM model, which constrains severely the viability of such theories. In [145], the Ricci and modified Ricci dark energy models were considered in the context of $f(R,T)$ gravity, and it was found that specific models can reproduce the expansion history of the Universe in accordance with the present observational data. The Dolgov–Kawasaki stability condition were also obtained for the specific reconstructed $f(R,T)$ functions.

6. $\mathbf{f}\left(\mathbf{R},\mathbf{T},{\mathbf{R}}_{\mathbf{\mu }\mathbf{\nu }}{\mathbf{T}}^{\mathbf{\mu }\mathbf{\nu }}\right)$ Gravity

6.1. Action and Field Equations

For the specific case of a traceless energy-momentum tensor, $T=0$, for instance, when the electromagnetic field is involved, the gravitational field equations for the $f(R,T)$ theory [16] reduce to that of the field equations for $f\left(R\right)$ gravity, and all non-minimal couplings of gravity to the matter field vanish. This motivates a further generalization of $f(R,T)$ gravity that consists of including an explicit first order coupling between the matter energy-momentum $T_{\mu \nu }$ and the Ricci tensor [25,26]. In contrast to $f(R,T)$ gravity, for $T=0$, this extra coupling still has a non-minimal coupling to the electromagnetic field via the $R_{\mu \nu }{T}^{\mu \nu }$ coupling term in the action, which is non-zero in general. As in the previous section, we consider the matter content as consisting of a perfect fluid characterized by its energy density and thermodynamic pressure only.

The action, is given by:

$$S=16\pi G\int d^{4}x\sqrt{-g}f\left(R,T,R_{\mu \nu }{T}^{\mu \nu }\right)+\int d^{4}x\sqrt{-g}{L}_m$$

(152)

The only requirement imposed on the function $f\left(R,T,R_{\mu \nu }{T}^{\mu \nu }\right)$ is that it is an arbitrary analytical function in all arguments. The gravitational field equations take the following form:

$$\begin{array}{c}\hfill (f_R-f_{RT}{L}_m)G_{\mu \nu }+\left[\square f_R+\frac{1}{2}R{f}_R-\frac{1}{2}f+f_T{L}_m+\frac{1}{2}{\nabla }_{\alpha }{\nabla }_{\beta }\left(f_{RT}{T}^{\alpha \beta }\right)\right]g_{\mu \nu }\end{array}$$

$$\begin{array}{c}\hfill -{\nabla }_{\mu }{\nabla }_{\nu }{f}_R+\frac{1}{2}\square \left(f_{RT}{T}_{\mu \nu }\right)+2f_{RT}{R}_{\alpha (\mu }{T}_{\nu )}^{\alpha }-{\nabla }_{\alpha }{\nabla }_{(\mu }\left[T_{\nu )}^{\alpha }{f}_{RT}\right]\end{array}$$

$$\begin{array}{c}\hfill -\left(f_T+\frac{1}{2}{f}_{RT}R+8\pi G\right)T_{\mu \nu }-2\left(f_T{g}^{\alpha \beta }+f_{RT}{R}^{\alpha \beta }\right)\frac{{\partial }^2{L}_m}{\partial g^{\mu \nu }\partial g^{\alpha \beta }}=0\end{array}$$

(153)

The trace of the gravitational field equation, Equation (153), is obtained as:

$$\begin{array}{c}\hfill 3\square f_R+\frac{1}{2}\square \left(f_{RT}T\right)+{\nabla }_{\alpha }{\nabla }_{\beta }\left(f_{RT}{T}^{\alpha \beta }\right)+R{f}_R-T{f}_T-\frac{1}{2}RT{f}_{RT}+2R_{\alpha \beta }{T}^{\alpha \beta }{f}_{RT}\\ \hfill +R{f}_{RT}{L}_m+4f_T{L}_m-2f-8\pi GT-2g^{\mu \nu }\left(g^{\alpha \beta }{f}_T+R^{\alpha \beta }{f}_{RT}\right)\frac{{\partial }^2{L}_m}{\partial g^{\mu \nu }\partial g^{\alpha \beta }}=0\end{array}$$

(154)

The second derivative of the matter Lagrangian with respect to the metric is non-zero if the matter Lagrangian is of second or higher order in the metric. Thus, for a perfect fluid with $L_m=-\rho$, or a scalar field with $L_m=-{\partial }_{\mu }\varphi {\partial }^{\mu }\varphi /2$, this term can be dropped. However, for instance, considering the Maxwell field, we have $L_m=-F_{\mu \nu }{F}^{\mu \nu }/4$; this term results in ${\partial }^2{L}_m/\partial g^{\mu \nu }\partial g^{\alpha \beta }=-F_{\mu \alpha }{F}_{\nu \beta }/2$, thus giving a non-zero contribution to the field equations.

In analogy with the standard Einstein field equation, one can write the field Equation (153) as:

$$G_{\mu \nu }=8\pi G_{\mathrm{eff}}{T}_{\mu \nu }-{\Lambda }_{\mathrm{eff}}{g}_{\mu \nu }+T_{\mu \nu }^{\mathrm{eff}}$$

(155)

where we have defined the effective gravitational coupling $G_{\mathrm{eff}}$, the effective cosmological constant ${\Lambda }_{\mathrm{eff}}$ and an effective energy-momentum tensor $T_{\mu \nu }^{\mathrm{eff}}$ as:

$$G_{\mathrm{eff}}=\frac{G+\frac{1}{8\pi }(f_T+\frac{1}{2}{f}_{RT}R-\frac{1}{2}\square f_{RT})}{f_R-f_{RT}{L}_m}$$

(156)

$${\Lambda }_{\mathrm{eff}}=\frac{2\square f_R+R{f}_R-f+2f_T{L}_m+{\nabla }_{\alpha }{\nabla }_{\beta }\left(f_{RT}{T}^{\alpha \beta }\right)}{2(f_R-f_{RT}{L}_m)}$$

(157)

and:

$$\begin{array}{c}\hfill T_{\mu \nu }^{\mathrm{eff}}=\frac{1}{f_R-f_{RT}{L}_m}\{{\nabla }_{\mu }{\nabla }_{\nu }{f}_R-{\nabla }_{\alpha }{f}_{RT}{\nabla }^{\alpha }{T}_{\mu \nu }-\frac{1}{2}{f}_{RT}\square T_{\mu \nu }-2f_{RT}{R}_{\alpha (\mu }{T}_{\nu )}^{\alpha }\\ \hfill +{\nabla }_{\alpha }{\nabla }_{(\mu }\left[T_{\nu )}^{\alpha }{f}_{RT}\right]+2\left(f_T{g}^{\alpha \beta }+f_{RT}{R}^{\alpha \beta }\right)\frac{{\partial }^2{L}_m}{\partial g^{\mu \nu }\partial g^{\alpha \beta }}\}\end{array}$$

(158)

respectively. Note that, in general, $G_{\mathrm{eff}}$ and ${\Lambda }_{\mathrm{eff}}$ are not constants and depend on the specific model considered.

6.2. Equation of Motion of the Massive Test Particles in the $f\left(R,T,R_{\mu \nu }{T}^{\mu \nu }\right)$ Gravity Theory

The covariant divergence of the energy-momentum tensor can be obtained by taking the divergence of the gravitational field equation, Equation (153), which takes the following form:

$$\begin{array}{c}\hfill {\nabla }^{\mu }{T}_{\mu \nu }=\frac{2}{\left(1+R{f}_{TR}+2f_T\right)}\{{\nabla }_{\mu }\left(f_{RT}{R}^{\sigma \mu }{T}_{\sigma \nu }\right)+{\nabla }_{\nu }\left(L_m{f}_T\right)-\frac{1}{2}\left(f_{RT}{R}_{\rho \sigma }+f_T{g}_{\rho \sigma }\right){\nabla }_{\nu }{T}^{\rho \sigma }\\ \hfill -G_{\mu \nu }{\nabla }^{\mu }\left(f_{RT}{L}_m\right)-\frac{1}{2}\left[{\nabla }^{\mu }\left(R{f}_{RT}\right)+2{\nabla }^{\mu }{f}_T\right]T_{\mu \nu }\}\end{array}$$

(159)

where we have assumed that ${\partial }^2{L}_m/\partial g^{\mu \nu }\partial g^{\alpha \beta }=0$ and have used the mathematical identities:

$$\begin{array}{c}\hfill {\nabla }^{\mu }\left(f_R{R}_{\mu \nu }+\square f_R{g}_{\mu \nu }-\frac{1}{2}f{g}_{\mu \nu }-{\nabla }_{\mu }{\nabla }_{\nu }{f}_R\right)=-\frac{1}{2}\left[f_T{\nabla }_{\nu }T+f_{RT}{\nabla }_{\nu }\left(R_{\rho \sigma }{T}^{\rho \sigma }\right)\right]\end{array}$$

(160)

$$2T_{\mu \tau ;\delta [;\rho ;\sigma ]}=T_{\mu \tau ;\alpha }{R}_{\delta \rho \sigma }^{\alpha }+T_{\alpha \tau ;\delta }{R}_{\mu \rho \sigma }^{\alpha }+T_{\mu \alpha ;\delta }{R}_{\tau \rho \sigma }^{\alpha }$$

(161)

and $\left[\square ,{\nabla }_{\nu }\right]T=R_{\mu \nu }{\nabla }^{\mu }T$, respectively.

In order to find the equation of motion for a massive test particle, as in the previous sections, we consider the energy-momentum tensor of a perfect fluid. Following the procedure outlined above, we obtain the equation of motion for a massive test particle, considering the matter Lagrangian $L_m=p$, as:

$$\frac{d^2{x}^{\lambda }}{d{s}^2}+{\Gamma }_{\mu \nu }^{\lambda }{U}^{\mu }{U}^{\nu }=f^{\lambda }$$

(162)

where the extra force acting on the test particles is given by:

$$\begin{array}{ccc}\hfill f^{\lambda }& =& \frac{1}{\rho +p}[\left(f_T+R{f}_{RT}\right){\nabla }_{\nu }\rho -\left(1+3f_T\right){\nabla }_{\nu }p-(\rho +p)f_{RT}{R}^{\sigma \rho }\left({\nabla }_{\nu }{h}_{\sigma \rho }-2{\nabla }_{\rho }{h}_{\sigma \nu }\right)\hfill \\ & & -f_{RT}{R}_{\sigma \rho }{h}^{\sigma \rho }{\nabla }_{\nu }\left(\rho +p\right)]\frac{h^{\lambda \nu }}{1+2f_T+R{f}_{RT}}\hfill \end{array}$$

(163)

Contrary to the nonminimal coupling presented in [9] and as can be seen from the above equations, the extra force does not vanish, even with the Lagrangian $L_m=p$.

The extra-force is perpendicular to the four-velocity, satisfying the relation $f^{\mu }{U}_{\mu }=0$. In the absence of any coupling between matter and geometry, with $f_T=f_{RT}=0$, the extra-force takes the usual form of the standard general relativistic fluid motion, i.e., $f^{\lambda }=-h^{\lambda \nu }{\nabla }_{\nu }p/\left(\rho +p\right)$. In the case of $f\left(R,T,R_{\mu \nu }{T}^{\mu \nu }\right)$ gravity, there is an explicit dependence of the extra-force on the Ricci tensor $R_{\sigma \rho }$, which makes the deviation from the geodesic motion more important for regions with strong curvatures.

6.3. Cosmological Applications: Specific Case of $f=R+\alpha R_{\mu \nu }{T}^{\mu \nu }$

Let us now consider some examples of cosmological solutions of the theory. In order to obtain explicit results and as the first step, one has to fix the functional form of the function $f\left(R,T,R_{\mu \nu }{T}^{\mu \nu }\right)$. We analyze the evolution and dynamics of the Universe, assuming that the Universe is isotropic and homogeneous, with the matter content described by the energy density ρ and thermodynamic pressure p with the matter Lagrangian chosen as $L_m=-\rho$. Consider the Friedmann–Lemaitre–Robertson–Walker (FLRW) metric, the Hubble parameter $H=\dot{a}/a$ and the deceleration parameter q, defined as:

$$q=\frac{d}{dt}\frac{1}{H}-1$$

(164)

Note that if $q<0$, the expansion of the Universe is accelerating, while positive values of $q>0$ describe decelerating evolutions.

For simplicity, consider the case in which the interaction between matter and geometry takes place only via the coupling between the energy-momentum and Ricci tensors, i.e.,

$$f=R+\alpha R_{\mu \nu }{T}^{\mu \nu }$$

(165)

In order to pass the Solar System and the other astrophysical tests, the correction term in Equation (165) must be small, which implies that α is a small parameter. This simple case serves as an example to show the main differences of the present theory with $f(R,T)$ gravity, considered before [16]. The gravitational field equations for this form of f are given by:

$$\begin{array}{c}\hfill G_{\mu \nu }+\alpha [2R_{\sigma (\mu }{T}_{\nu )}^{\sigma }-\frac{1}{2}{R}_{\rho \sigma }{T}^{\rho \sigma }{g}_{\mu \nu }-\frac{1}{2}R{T}_{\mu \nu }-\frac{1}{2}\left(2{\nabla }_{\sigma }{\nabla }_{(\nu }{T}_{\mu )}^{\sigma }-\square T_{\mu \nu }-{\nabla }_{\alpha }{\nabla }_{\beta }{T}^{\alpha \beta }{g}_{\mu \nu }\right)\\ \hfill -G_{\mu \nu }{L}_m-2R^{\alpha \beta }\frac{{\partial }^2{L}_m}{\partial g^{\mu \nu }\partial g^{\alpha \beta }}]-8\pi G{T}_{\mu \nu }=0\end{array}$$

(166)

For the case of the FLRW metric, the modified Friedmann equations are:

$$3H^2=\frac{\kappa }{1-\alpha \rho }\rho +\frac{3}{2}\frac{\alpha }{1-\alpha \rho }H\left(\dot{\rho }-\dot{p}\right)$$

(167)

and:

$$2\dot{H}+3H^2=\frac{2\alpha }{1+\alpha p}H\dot{\rho }-\frac{\kappa p}{1+\alpha p}+\frac{1}{2}\frac{\alpha }{1+\alpha p}\left(\ddot{\rho }-\ddot{p}\right)$$

(168)

respectively, where we have denoted $\kappa =8\pi G$ for simplicity. When $\alpha =0$, we recover the standard Friedmann equations. To remove the under determinacy of the field equations, we impose an equation of state for the cosmological matter, $p=p\left(\rho \right)$. A standard form of the cosmological matter equation of state is $p=\omega \rho$, where $\omega =\mathrm{constant}$ and $0\le \omega \le 1$.

6.3.1. High Cosmological Density Limit of the Field Equations

First, we consider the high energy density limit of the system of modified cosmological equations Equations (167) and (168) and assume that the constant α is small, so that $\alpha \rho \ll 1$ and $\alpha p\ll 1$, respectively. In the high-energy limit, $\rho =p$, then Equation (167) and (168) take the approximate form:

$$\begin{array}{ccc}\hfill 3H^2& =& \kappa \rho \hfill \end{array}$$

(169)

$$\begin{array}{ccc}\hfill 2\dot{H}+3H^2& =& -\kappa \rho +2\alpha H\dot{\rho }\hfill \end{array}$$

(170)

The time evolution of the Hubble parameter is described by the equation:

$$\left(1-\frac{6\alpha }{\kappa }\right)H^2+3H^2=0$$

(171)

and, hence, for this model, the evolution of the Hubble parameter is given by:

$$H\left(t\right)=\frac{\sqrt{\left(C_1+3\kappa t\right){}^2-24\alpha \kappa }+C_1+3\kappa t}{12\alpha }$$

(172)

where $C_1>0$ is an arbitrary integration constant. The scale factor of the Universe is given by:

$$a\left(t\right)=C_2\frac{exp\left[\frac{\left(C_1+3\kappa t\right)\sqrt{\left(C_1+3\kappa t\right){}^2-24\alpha \kappa }+9\kappa t^2+6\kappa C_{1}t}{72\alpha \kappa }\right]}{\sqrt[3]{\sqrt{\left(C_1+3\kappa t\right){}^2-24\alpha \kappa }+C_1+3\kappa t}}$$

(173)

where $C_2>0$ is an arbitrary integration constant. In order to have a physical solution, the integration constant $C_1$ must satisfy the constraint $C_1\ge \sqrt{24\kappa \alpha }$.

The values of the integration constant can be determined from the condition $H\left(0\right)=H_0$, and $a\left(0\right)=a_0$, where $H_0$ and $a_0$ are the initial values of the Hubble parameter and of the scale factor of the Universe, respectively. This condition immediately provides for $C_1$ the following value $C_1=(6\alpha H_0^2+\kappa )/H_0$. For the integration constant $C_2$, we obtain:

$$\begin{array}{c}\hfill C_2=a_0\sqrt[3]{\sqrt{\frac{{\left(\kappa -6\alpha H_0^2\right)}^2}{H_0^2}}+6\alpha H_0+\frac{\kappa }{H_0}}exp\left[-\frac{\sqrt{\frac{{\left(\kappa -6\alpha H_0^2\right)}^2}{H_0^2}}\left(6\alpha H_0^2+\kappa \right)}{72\alpha H_0\kappa }\right]\end{array}$$

(174)

In the small time limit, the scale factor is approximated by:

$$a\left(t\right)\approx a_0\left(1+\frac{\kappa }{6H_0\alpha }t\right)$$

(175)

The deceleration parameter is obtained as:

$$q\left(t\right)=-\frac{36\alpha H_0\kappa }{\sqrt{\frac{{\left(6\alpha H_0^2+3H_0\kappa t+\kappa \right)}^2}{H_0^2}-24\alpha \kappa }\left[6\alpha H_0^2+H_0\sqrt{\frac{{\left(6\alpha H_0^2+3H_0\kappa t+\kappa \right)}^2}{H_0^2}-24\alpha \kappa }+3H_0\kappa t+\kappa \right]}-1$$

(176)

and it can be represented in a form of a power series as:

$$q\left(t\right)\approx -1-\frac{18\alpha H_0^2}{\kappa -6\alpha H_0^2}+\frac{6H_0{\kappa }^2}{{\left(\kappa -6\alpha H_0^2\right)}^3}t$$

(177)

For small values of time, if $24\alpha H_0^2\ll \kappa$, $q\approx -1$, and the Universe starts its expansion from a de Sitter-like phase, entering, after a finite time interval, into a decelerating phase. On the other hand, if $\kappa >6\alpha H_0^2$, $q<-1$, and the non-singular Universe experiences an initial super-accelerating phase.

6.3.2. The Case of Dust

Next, we consider the case of low density cosmological matter, with $p=0$. Moreover, we assume again that the condition $\alpha \rho \ll 1$ holds. Then, the gravitational field equations, Equations (167) and (168), corresponding to a FLRW universe, take the approximate form:

$$\begin{array}{ccc}\hfill 3H^2& =& \kappa \rho +\frac{3}{2}\alpha H\dot{\rho }\hfill \end{array}$$

(178)

$$\begin{array}{ccc}\hfill 2\dot{H}+3H^2& =& 2\alpha H\dot{\rho }+\frac{1}{2}\alpha \ddot{\rho }\hfill \end{array}$$

(179)

First, we consider the matter dominated phase of the model, in which the non-accelerating expansion of the Universe can be described by a power law form of the scale factor, so that $a=t^m$, $m=\mathrm{constant}$ and $H=m/t$, respectively. The deceleration parameter is given by $q=1/m-1$. Therefore, Equation (178) gives, for the time evolution of the density, the equation:

$$\frac{3\alpha m}{2t}\dot{\rho }+\kappa \rho -3\frac{m^2}{t^2}=0$$

(180)

with the general solution given by:

$$\rho \left(t\right)=\frac{e^{-\frac{\kappa t^2}{3\alpha }}\left[3{\rho }_0\alpha e^{\frac{\kappa t_0^2}{3\alpha }}+\text{Ei}\left(\frac{t^2\kappa }{3\alpha }\right)-\text{Ei}\left(\frac{t_0^2\kappa }{3\alpha }\right)\right]}{3\alpha }$$

(181)

where $\text{Ei}\left(z\right)=-{\int }_{-z}^{\infty }{e}^{-t}dt/t$ is the exponential integral function, and we have used the initial condition $\rho \left(t_0\right)={\rho }_0$. By substituting the expressions of the density and of the Hubble parameter into Equation (179), in the first order of approximation, we obtain the following constraint on m, $(9m^2-10m+1)/3t^2+O\left(t^2\right)\approx 0$, which is (approximately) satisfied if m is given by the algebraic equation $9m^2-10m+1=0$, having the solutions $m_1=1$ and $m_2=1/9$, respectively. The deceleration parameters corresponding to these solutions are $q_1=0$ and $q_2=8$, respectively. Since a value of the deceleration parameter of the order of $q=8$ seems to be ruled out by the observations, the physical solution has a scale factor $a=t$, and $q=0$. The cosmological solutions with zero value of the deceleration parameter are called marginally accelerating, and they describe the pre-accelerating phase of the cosmic expansion.

Now, we look for a de Sitter-type solution of the field equations for the pressureless matter, Equations (178) and (179), by taking $H=H_0=\mathrm{constant}$. Then, it follows that, in order to have an accelerated expansion, the matter density must satisfy the equation:

$$\ddot{\rho }-H_0\dot{\rho }+\frac{2\kappa }{\alpha }\rho =0$$

(182)

with the general solution given by:

$$\begin{array}{ccc}\hfill \rho \left(t\right)& =& e^{\frac{1}{2}{H}_0\left(t-t_0\right)}\{\frac{\sqrt{\alpha }\left(2{\rho }_{01}-H_0{\rho }_0\right)}{\sqrt{\alpha H_0^2-8\kappa }}sinh\left[\frac{\sqrt{\alpha H_0^2-8\kappa }}{2\sqrt{\alpha }}\left(t-t_0\right)\right]\hfill \\ & & +cosh\left[\frac{\sqrt{\alpha H_0^2-8\kappa }}{2\sqrt{\alpha }}\left(t-t_0\right)\right]\}\hfill \end{array}$$

(183)

where we have used the initial conditions $\rho \left(t_0\right)={\rho }_0$ and $\dot{\rho }\left(t_0\right)={\rho }_{01}$, respectively. Therefore, in the presence of a non-trivial curvature-matter coupling, once the evolution of the matter density is given by Equation (183), the time evolution of the Universe is of the de Sitter-type.

7. Dark Matter As a Curvature-Matter Coupling Effect

According to Newton’s gravitation theory, the rotation of hydrogen clouds around galaxies should show a Keplerian decrease with distance r of the orbital rotational speed $v_{tg}$ outside the luminous baryonic matter, $v_{tg}^2\propto M\left(r\right)/r$, where $M\left(r\right)$ is the dynamical mass of the galaxy. However, observations show that the rotation curves are rather flat [146,147,148,149]. The rotational velocities increase near the center of the galaxy, and then, they remain nearly constant at a value of $v_{tg\infty }\sim$ 200–300 km/s. Hence, the observations give a general mass profile $M\left(r\right)\approx r{v}_{tg\infty }^2/G$ [146,147,148,149]. Consequently, even at large distances, where very little luminous matter can be detected, the mass within a distance r from the center of the galaxy increases linearly with the radial distance r.

This unusual behavior of the galactic rotation curves is usually explained by postulating the existence of some dark (invisible) matter, distributed in a spherical halo around the galaxies. The rotation curves are obtained by measuring the frequency shifts z of the 21-cm radiation emission from the neutral hydrogen gas clouds. Usually, astronomers report the resulting z in terms of a velocity field $v_{tg}$ [146,147,148]. In the standard cosmological model, dark matter is assumed to be a cold, pressureless medium. Many possible candidates for dark matter have been proposed, the most popular ones being the weakly interacting massive particles (WIMP) (for a review of the particle physics aspects of dark matter, see [150]). The interaction cross-section with normal baryonic matter is expected to be small, but non-zero, and we expect that dark matter particles can be detected directly. However, after more than 20 years of intense experimental and observational effort, presently, no convincing non-gravitational evidence for dark matter exists.

In the present section, we will briefly analyze, following [42], the possibility that the dynamic behavior of hydrogen clouds rotating around the galactic center can be explained due to the presence of a curvature-matter coupling. We will restrict our analysis to the case of modified gravity models with linear curvature-matter coupling. The extra-terms in the gravitational field equations modify the equations of motion of test particles and induce a supplementary gravitational interaction, which can account for the observed behavior of the galactic rotation curves.

Note that the modified gravity approach has been explored, as an alternative to dark matter, in different contexts [151,152,153,154], such as in $f\left(R\right)$ gravity and the recently proposed hybrid metric-Palatini gravity [155,156,157,158].

7.1. Stable Circular Orbits and Frequency Shifts in Modified Gravity with Linear Curvature-Matter Coupling

In the following, we assume that gas clouds behave like massive test particles, moving in a static and spherically-symmetric space-time outside the galactic baryonic mass distribution. In the galactic space-time, we consider two observers $O_E$ and $O_{\infty }$, with four-velocities $u_E^{\mu }$ and $u_{\infty }^{\mu }$, respectively. Observer $O_E$ is the light emitter (i.e., the gas clouds placed at a point $P_E$ of the space-time), and $O_{\infty }$ represents the detector at point $P_{\infty }$, located far from the emitter, at “spatial infinity” [159].

Without any loss of generality, we can assume that the gas clouds move in the galactic plane $\theta =\pi /2$, so that $u_E^{\mu }={\left(\dot{t},\dot{r},0,\dot{\varphi }\right)}_E$, where the dot stands for derivation with respect to the affine parameter s. On the other hand, we suppose that the detector is static (i.e., $O_{\infty }$’s four-velocity is tangent to the static Killing field $\partial /\partial t$), and in the chosen coordinate system, its four-velocity is $u_{\infty }^{\mu }={\left(\dot{t},0,0,0\right)}_{\infty }$[159]. In this section, we use the $(+,-,-,-)$ signature.

The static spherically symmetric metric outside the galactic baryonic mass distribution is given by:

$$d{s}^2=e^{\nu \left(r\right)}d{t}^2-e^{\lambda \left(r\right)}d{r}^2-r^2\left(d{\theta }^2+{sin}^2\theta d{\varphi }^2\right)$$

(184)

where the metric coefficients ν and λ are functions of the radial coordinate r only. The motion of a test particle in the gravitational field in modified gravity with a linear curvature-matter coupling is described by the Lagrangian given by Equation (24):

$$L=Q\left[e^{\nu \left(r\right)}{\left(\frac{dt}{ds}\right)}^2-e^{\lambda \left(r\right)}{\left(\frac{dr}{ds}\right)}^2-r^2{\left(\frac{d\Omega }{ds}\right)}^2\right]$$

(185)

where $d{\Omega }^2=d{\theta }^2+{sin}^2\theta d{\varphi }^2$, and the function Q is given by:

$$Q=\left[1+\zeta f_2\left(R\right)\right]\frac{d{L}_m\left(\rho \right)}{d\rho }$$

(186)

where ζ is a constant.

For $\theta =\pi /2$, $d{\Omega }^2=d{\varphi }^2$. From the Lagrange equations, it follows that we have two constants of motion, the energy E:

$$E=Q{e}^{\nu \left(r\right)}\dot{t}$$

(187)

and the angular momentum l, given by:

$$l=Q{r}^2\dot{\varphi }$$

(188)

where a dot denotes the derivative with respect to the affine parameter s. The condition of the normalization of the four-velocity $U^{\mu }{U}_{\mu }=1$ gives:

$$1=e^{\nu \left(r\right)}{\dot{t}}^2-e^{\lambda \left(r\right)}{\dot{r}}^2-r^2{\dot{\varphi }}^2$$

(189)

from which, with the use of the constants of motion, we obtain:

$$E^2=Q^2{e}^{\nu +\lambda }{\dot{r}}^2+e^{\nu }\left(\frac{l^2}{r^2}+Q^2\right)$$

(190)

Equation (190) shows that the radial motion of the particles in modified gravity with linear matter-geometry coupling is the same as that of a particle in ordinary Newtonian mechanics, with velocity $\dot{r}$, position-dependent mass $m_{eff}=2Q^2{e}^{\nu +\lambda }$ and energy $E^2$, respectively, moving in the effective potential:

$$V_{eff}\left(r\right)=e^{\nu \left(r\right)}\left(\frac{l^2}{r^2}+Q^2\right)$$

(191)

The conditions for stable circular orbits $\partial V_{eff}/\partial r=0$ and $\dot{r}=0$ determine the energy and the angular momentum of the particle as:

$$l^2=\frac{1}{2}\frac{r^{3}Q\left({\nu }^{\prime }Q+2Q^{\prime }\right)}{1-r{\nu }^{\prime }/2}$$

(192)

and:

$$E^2=\frac{e^{\nu }Q\left(r{Q}^{\prime }+Q\right)}{1-r{\nu }^{\prime }/2}$$

(193)

respectively.

The line element given by Equation (184) can be rewritten in terms of the spatial components of the velocity, normalized with the speed of light, measured by an inertial observer far from the source, as:

$$d{s}^2=d{t}^2\left(1-v^2\right)$$

(194)

where:

$$v^2=e^{-\nu }\left[e^{\lambda }{\left(\frac{dr}{dt}\right)}^2+r^2{\left(\frac{d\Omega }{dt}\right)}^2\right]$$

(195)

For a stable circular orbit $dr/dt=0$, and the tangential velocity of the test particle is obtained as:

$$v_{tg}^2=e^{-\nu }{r}^2{\left(\frac{d\Omega }{dt}\right)}^2$$

(196)

In terms of the conserved quantities E and l, the angular velocity is given, for $\theta =\pi /2$, by:

$$v_{tg}^2=\frac{e^{\nu }}{r^2}\frac{l^2}{E^2}$$

(197)

With the use of Equations (192) and (193), we obtain:

$$v_{tg}^2=\frac{1}{2}\frac{r\left({\nu }^{\prime }Q+2Q^{\prime }\right)}{r{Q}^{\prime }+Q}$$

(198)

Thus, the rotational velocity of the test body in modified gravity with linear coupling between matter and geometry is determined by the metric coefficient $exp\left(\nu \right)$, and by the function Q and its derivative with respect to the radial coordinate r. In the standard general relativistic limit $\zeta =0$, $Q=1$, and we obtain $v_{tg}^2=r{\nu }^{\prime }/2$.

7.2. The Effect of the Curvature-Matter Coupling on the Light Shifts

The velocities of the hydrogen clouds rotating around galactic centers are determined from the red and blue shifts of the radiation emitted by the gas moving on circular orbits on both sides of the central region. The radiation travels on null geodesics with tangent $k^{\mu }$. We may restrict without any loss of generality $k^{\mu }$ to lie in the equatorial plane $\theta =\pi /2$ and evaluate the frequency shift for a light signal emitted from $O_E$, in motion, in a circular orbit, and detected by $O_{\infty }$ at infinity. The frequency shift associated with the emission and detection of the light signal is given by:

$$z=1-\frac{{\omega }_E}{{\omega }_{\infty }}$$

(199)

where ${\omega }_I=-k_{\mu }{u}_I^{\mu }$, and the index I refers to emission ($I=E$) or detection ($I=\infty$) at the corresponding space-time point [159,160]. Two frequency shifts, corresponding to the frequency shifts of a receding or approaching gas cloud, respectively, are associated with light propagation in the same and opposite direction of motion of the emitter, respectively. In terms of the tetrads $e_{\left(0\right)}=e^{-\nu /2}\partial /\partial t$, $e_{\left(1\right)}=e^{-\lambda /2}\partial /\partial r$, $e_{\left(2\right)}=r^{-1}\partial /\partial \theta$, $e_{\left(3\right)}={\left(rsin\theta \right)}^{-1}\partial /\partial \varphi$, the frequency shifts take the form [159]:

$$z_{\pm }=1-e^{\left[{\nu }_{\infty }-\nu \left(r\right)\right]/2}\left(1\mp v\right)\Gamma$$

(200)

where $v={\left[{\sum }_{i=1}^3{\left(u_{\left(i\right)}/u_{\left(0\right)}\right)}^2\right]}^{1/2}$, with $u_{\left(i\right)}$ the components of the particle’s four velocity along the tetrad (i.e., the velocity measured by an Eulerian observer, whose world-line is tangent to the static Killing field). $\Gamma ={\left(1-v^2\right)}^{-1/2}$ is the usual Lorentz factor, and $exp\left({\nu }_{\infty }\right)$ is the value of $exp\left[\nu \left(\left(r\right)\right)\right]$ for $r\to \infty$. In the case of circular orbits in the $\theta =\pi /2$ plane, we obtain:

$$z_{\pm }=1-e^{\left[{\nu }_{\infty }-\nu \left(r\right)\right]/2}\frac{1\mp \sqrt{r\left({\nu }^{\prime }Q+2Q^{\prime }\right)/2\left(r{Q}^{\prime }+Q\right)}}{\sqrt{1-r\left({\nu }^{\prime }Q+2Q^{\prime }\right)/2\left(r{Q}^{\prime }+Q\right)}}$$

(201)

It is convenient to define two other quantities, $z_D=\left(z_{+}-z_{-}\right)/2$ and $z_A=\left(z_{+}+z_{-}\right)/2$, respectively [159]. In the modified gravity model with a linear curvature-matter coupling, the redshift factors are given by:

$$z_D\left(r\right)=e^{\left[{\nu }_{\infty }-\nu \left(r\right)\right]/2}\frac{\sqrt{r\left({\nu }^{\prime }Q+2Q^{\prime }\right)/2\left(r{Q}^{\prime }+Q\right)}}{\sqrt{1-r\left({\nu }^{\prime }Q+2Q^{\prime }\right)/2\left(r{Q}^{\prime }+Q\right)}}$$

(202)

and:

$$z_A\left(r\right)=1-\frac{e^{\left[{\nu }_{\infty }-\nu \left(r\right)\right]/2}}{\sqrt{1-r\left({\nu }^{\prime }Q+2Q^{\prime }\right)/2\left(r{Q}^{\prime }+Q\right)}}$$

(203)

respectively, which can be easily related to the galactic observations [159]. $z_A$ and $z_D$ satisfy the relation:

$${\left(z_A-1\right)}^2-z_D^2=exp\left[2\left({\nu }_{\infty }-\nu \left(r\right)\right)\right]$$

(204)

and thus, in principle, by assuming that the metric tensor component $exp\left[\nu \left(\left(r\right)\right)\right]$ is known, Q and $Q^{\prime }$ can be obtained directly from the astrophysical observations. Hence, the observations of the red- and blue-shifts of the radiation emitted by hydrogen clouds rotating around the galactic center could provide a direct observational test of the galactic geometry and, implicitly, of the modified gravity models with linear coupling between matter and geometry.

7.3. Galactic Rotation Curves and the Curvature-Matter Coupling

The tangential velocities $v_{tg}$ of gas clouds, considered as massive test particles moving around the centers of the galaxies, cannot be measured directly. Instead, they can be inferred from the redshift $z_{\infty }$ of the radiation observed at spatial infinity, for which:

$$1+z_{\infty }=exp\left[\left({\nu }_{\infty }-\nu \right)/2\right]\left(1\pm v_{tg}\right)/\sqrt{1-v_{tg}^2}$$

(205)

Since the velocities of the gas clouds are non-relativistic, with $v_{tg}\le \left(4/3\right)\times {10}^{-3}$, we observe $v_{tg}\approx z_{\infty }$ (as the first term of a geometric series), with the consequence that the lapse function $exp\left(\nu \right)$ necessarily tends at infinity to unity, i.e., $e^{\nu }\approx e^{{\nu }_{\infty }}/\left(1-v_{tg}^2\right)\approx e^{{\nu }_{\infty }}\to 1$. The observations show that at distances large enough from the galactic center, $v_{tg}\approx$ constant.

In the following, we use this observational constraint to reconstruct the coupling term between curvature and matter in the “dark matter”-dominated region, far away from the baryonic matter distribution. By assuming that $v_{tg}=$ constant, Equation (198) can be written as:

$$v_{tg}^2\frac{1}{rQ}\frac{d}{dr}\left(rQ\right)=\frac{{\nu }^{\prime }}{2}+\frac{Q^{\prime }}{Q}$$

(206)

and can be immediately integrated to give:

$$Q\left(r\right)={\left(\frac{r}{r_0}\right)}^{v_{tg}^2/\left(1-v_{tg}^2\right)}exp\left[-\frac{\nu }{2\left(1-v_{tg}^2\right)}\right]$$

(207)

where $r_0$ is an arbitrary constant of integration. Since hydrogen clouds are a pressureless dust ($p=0$) that can be characterized by their density ρ only, the Lagrangian of the matter (gas cloud) is given by $L_{mat}\left(\rho \right)=\rho$. Therefore, from Equations (186) and (207), we obtain:

$$\begin{array}{ccc}\hfill \zeta f_2\left(R\right)& =& {\left(\frac{r}{r_0}\right)}^{v_{tg}^2/2\left(1-v_{tg}^2\right)}exp\left[-\frac{\nu }{4\left(1-v_{tg}^2\right)}\right]-1\approx \frac{v_{tg}^2}{2\left(1-v_{tg}^2\right)}ln\frac{r}{r_0}-\frac{\nu }{4\left(1-v_{tg}^2\right)}-\hfill \\ & & \frac{v_{tg}^2}{8\left(1-2v_{tg}^2\right)}\nu ln\frac{r}{r_0}\hfill \end{array}$$

(208)

For an arbitrary velocity profile $v_{tg}=v_{tg}\left(r\right)$, the general solution of Equation (206) is given by:

$$\sqrt{Q\left(r\right)}=1+\zeta f_2\left(R\right)=\sqrt{Q_0}exp\left[\frac{1}{2}\int \frac{v_{tg}^2\left(r\right)/r-{\nu }^{\prime }/2}{1-v_{tg}^2\left(r\right)}dr\right]$$

(209)

where $Q_0$ is an arbitrary constant of integration.

For $v_{tg}^2$, we assume the simple empirical dark halo rotational velocity law [161]:

$$v_{tg}^2=\frac{v_0^2{x}^2}{a^2+x^2}$$

(210)

where $x=r/r_{opt}$, $r_{opt}$ is the optical radius containing 83% of the galactic luminosity. The parameters a, the ratio of the halo core radius and $r_{opt}$ and the terminal velocity $v_0$ are all functions of the galactic luminosity L. For spiral galaxies $a=1.5{\left(L/L_{*}\right)}^{1/5}$ and $v_0^2=v_{opt}^2\left(1-{\beta }_{*}\right)\left(1+a^2\right)$, where $v_{opt}=v_{tg}\left(r_{opt}\right)$, and ${\beta }_{*}=0.72+0.44{log}_{10}\left(L/L_{*}\right)$, with $L_{*}={10}^{10.4}{L}_{\odot }$.

One can assume that the coupling between the neutral hydrogen clouds and the geometry is small, $\zeta f_2\left(R\right)L_{mat}\ll 1$, and consequently, the galactic geometry is not significantly modified in the vacuum outside the baryonic mass distribution with mass $M_B$, corresponding to $L_{mat}\approx 0$. Moreover, we assume, for simplicity, that outside the baryonic matter distribution, the galactic metric is given by the Schwarzschild metric (which is also a solution of the vacuum field equations of the $f\left(R\right)$-modified gravity [162] and still gives the dominant contribution to the total metric, even if the curvature-matter coupling is not small), written as:

$$e^{\nu }=e^{-\lambda }=1-\frac{2R_0}{x}$$

(211)

where $R_0=G{M}_B/r_{opt}$, from Equation (209), we obtain:

$$\begin{array}{c}\hfill 1+\zeta f_2\left(R\right)=exp\left[\alpha \times \mathrm{arctanh}\left(\frac{\sqrt{1-v_0^2}}{a}x\right)\right]\times {\left(1-\frac{2R_0}{x}\right)}^{\beta }\frac{{\left[\left(1-v_0^2\right)x^2+a^2\right]}^{\gamma }}{Q_0^{-1/2}{x}^{1/4}}\end{array}$$

(212)

where:

$$\alpha =-\frac{a{R}_0{v}_0^2}{2\sqrt{1-v_0^2}\left[a^2+4\left(1-v_0^2\right)R_0^2\right]}$$

(213)

$$\beta =\frac{a^2-4R_0^2}{4\left[a^2+4\left(1-v_0^2\right)R_0^2\right]}$$

(214)

and:

$$\gamma =-\frac{v_0^2\left[a^2+6\left(1-v_0^2\right)R_0^2\right]}{4\left(1-v_0^2\right)\left[a^2+4\left(1-v_0^2\right)R_0^2\right]}$$

(215)

respectively. Thus, the geometric part of the curvature-matter coupling can be completely reconstructed from the observational data on the galactic rotation curves.

As one can see from Equation (198), in the limit of large r, when ${\nu }^{\prime }\to 0$ (in the case of the Schwarzschild metric ${\nu }^{\prime }\approx 2M_B/r^2$), the tangential velocity of test particles at infinity is given by:

$$v_{tg\infty }^2=\frac{r}{r+Q/Q^{\prime }}$$

(216)

which, due to the presence of the curvature-matter coupling, does not decay to zero at large distances from the galactic center, a behavior that is perfectly consistent with the observational data and that is usually attributed to the existence of dark matter.

By using the simple observational fact of the constancy of the galactic rotation curves, the curvature-matter coupling function can be completely reconstructed, without any supplementary assumption. If, for simplicity, we consider again that the metric in the vacuum outside the galaxy can be approximated by the Schwarzschild metric, with $exp\left(\nu \right)=1-2G{M}_B/r$, where $M_B$ is the mass of the baryonic matter of the galaxy, then in the limit of large r, we have $\nu \to 0$. Therefore, from Equation (208), we obtain:

$$\zeta \underset{r\to \infty }{lim}{f}_2\left(R\right)\approx \frac{v_{tg}^2}{2\left(1-v_{tg}^2\right)}ln\frac{r}{r_0}$$

(217)

If the galactic rotation velocity profiles and the galactic metric are known, the coupling function can be reconstructed exactly over the entire mass distribution of the galaxy.

One can formally associate an approximate “dark matter” mass profile $M_{DM}\left(r\right)$ with the tangential velocity profile, which is determined by the non-minimal curvature-matter coupling and is given by:

$$M_{DM}\left(r\right)\approx \frac{1}{2G}\frac{r^2\left({\nu }^{\prime }Q+2Q^{\prime }\right)}{r{Q}^{\prime }+Q}$$

(218)

The corresponding “dark matter” density profile ${\rho }_{DM}\left(r\right)$ can be obtained as:

$$\begin{array}{c}\hfill {\rho }_{DM}\left(r\right)=\frac{1}{4\pi r^2}\frac{dM}{dr}=\frac{1}{4\pi G}\times \left[\frac{{\nu }^{\prime }Q+2Q^{\prime }}{r\left(r{Q}^{\prime }+Q\right)}+\frac{{\nu }^{\prime \prime }Q+\nu Q^{\prime }+2Q^{\prime \prime }}{2\left(r{Q}^{\prime }+Q\right)}\right.-\left.\frac{\left({\nu }^{\prime }Q+2Q^{\prime }\right)\left(r{Q}^{\prime \prime }+2Q^{\prime }\right)}{2{\left(r{Q}^{\prime }+Q\right)}^2}\right]\end{array}$$

(219)

7.4. Constraining the Curvature-Matter Coupling with Galactic Stellar Distributions

Other observational constraints on $M_{DM}$ and ${\rho }_{DM}$ can be obtained from the study of the galactic stellar populations. We assume that each galaxy consists of a single, pressure-supported stellar population that is in dynamic equilibrium and traces an underlying gravitational potential resulting from the non-minimal curvature-matter coupling. In spherical symmetry, the equivalent mass profile induced by the geometry-matter coupling (the mass profile of the “dark matter” halo) relates to the moments of the stellar distribution function via the Jeans equation [149]:

$$\frac{d}{dr}\left[{\rho }_s〈v_r^2〉\right]+\frac{2{\rho }_s\left(r\right)\beta }{r}=-\frac{G{\rho }_s{M}_{DM}\left(r\right)}{r^2}$$

(220)

where ${\rho }_s\left(r\right)$, $〈v_r^2〉$ and $\beta \left(r\right)=1-〈v_{\theta }^2〉/〈v_r^2〉$ describe the three-dimensional density, radial velocity dispersion and orbital anisotropy of the stellar component, where $〈v_{\theta }^2〉$ is the tangential velocity dispersion. By assuming that the anisotropy is a constant, the Jeans equation has the solution [163]:

$${\rho }_s〈v_r^2〉=G{r}^{-2\beta }{\int }_r^{\infty }{s}^{2(1-\beta )}{\rho }_s\left(s\right)M_{DM}\left(s\right)ds$$

(221)

With the use of Equation (218), we obtain for the stellar velocity dispersion equation:

$${\rho }_s〈v_r^2〉\approx \frac{1}{2}{r}^{-2\beta }{\int }_r^{\infty }{s}^{2(2-\beta )}{\rho }_s\left(s\right)\frac{{\nu }^{\prime }\left(s\right)Q\left(s\right)+2Q^{\prime }\left(s\right)}{s{Q}^{\prime }\left(s\right)+Q\left(s\right)}ds$$

(222)

After projection along the line of sight, the “dark matter” mass profile can be related to two observable profiles, the projected stellar density $I\left(R\right)$, and to the stellar velocity dispersion ${\sigma }_p\left(R\right)$, according to the relation [149]:

$${\sigma }_P^2\left(R\right)=\frac{2}{I\left(R\right)}{\int }_R^{\infty }\left(1-\beta \frac{R^2}{r^2}\right)\frac{{\rho }_s〈v_r^2〉r}{\sqrt{r^2-R^2}}dr$$

(223)

Given a projected stellar density model $I\left(R\right)$, one recovers the three-dimensional stellar density from [149]:

$${\rho }_s\left(r\right)=-(1/\pi ){\int }_r^{\infty }\left(\frac{dI}{dR}\right){\left(R^2-r^2\right)}^{-1/2}dR$$

(224)

Therefore, once the stellar density profile $I\left(R\right)$, the stellar velocity dispersion $〈v_r^2〉$ and the galactic metric are known, with the use of the integral Equation (222), one can obtain the explicit form of the curvature-matter coupling function Q and the equivalent mass profile induced by the non-minimal coupling between matter and curvature. The simplest analytic projected density profile is the Plummer profile [149], given by:

$$I\left(R\right)=\mathcal{L}{\left(\pi r_{half}^2\right)}^{-1}{\left(1+R^2/r_{half}^2\right)}^{-2}$$

(225)

where $\mathcal{L}$ is the total luminosity and $r_{half}$ is the projected half-light radius (the radius of the cylinder that encloses half of the total luminosity).

7.5. Stability of the Stable Circular Orbits in Modified Gravity with Curvature-Matter Coupling

An important physical requirement for the circular orbits of the test particles moving around galaxies is that they must be stable. Let $r_0$ be a circular orbit and consider a perturbation of it of the form $r=r_0+\delta$, where $\delta \ll r_0$ [160]. Taking expansions of $V_{eff}\left(r\right)$, $exp\left(\nu +\lambda \right)$ and $Q^2\left(r\right)$ about $r=r_0$, it follows from Equation (190) that:

$$\ddot{\delta }+\frac{1}{2}{Q}^2\left(r_0\right)e^{\nu \left(r_0\right)+\lambda \left(r_0\right)}{V}_{eff}^{\prime \prime }\left(r_0\right)\delta =0$$

(226)

The condition for the stability of the simple circular orbits requires $V_{eff}^{\prime \prime }\left(r_0\right)>0$[160]. This gives for the coupling function Q the constraint:

$$\begin{array}{c}\hfill {\left.\frac{d^{2}Q}{d{r}^2}\right|}_{r=r_0}<{\left.\left[{\nu }^{\prime \prime }\left(\frac{l^2}{r^2}+Q^2\right)+{\nu }^{\prime }\left(-\frac{2l^2}{r^2}+2Q{Q}^{\prime }\right)-\frac{6l^2}{r^4}\right]\right|}_{r=r_0}\end{array}$$

(227)

a condition that must be satisfied at any point $r_0$ of the galactic space-time.

From the observational, as well as from the theoretical point of view, an important problem is to estimate an upper bound for the cutoff of the constancy of the tangential velocities. If in the large r limit, the coupling function satisfies the condition $Q^{\prime }/Q\to 0$, then $Q/Q^{\prime }\to \infty$, and in this limit, the tangential velocity decays to zero. If the exact functional form of Q is known, the value of r at which the rotational velocity becomes zero can be accurately estimated.

8. Discussion and Conclusions

In the present paper, we have reviewed the cosmological and astrophysical applications of generalized curvature-matter couplings in $f\left(R\right)$-type gravity models, and we have extensively analysed some of their applications. Specific models were explored and presented in detail. The gravitational field equations in the metric formalism, in the presence of a nonminimal coupling between curvature and matter, were presented, as well as the equations of motion for test particles, which follow from the covariant divergence of the energy-momentum tensor. Generally, the motion is non-geodesic and takes place in the presence of an extra force orthogonal to the four-velocity. The Newtonian limit of the equation of motion was also described, and a procedure for obtaining the energy-momentum tensor of the matter in the framework of these gravity models was presented. On the other hand, it was shown that the gravitational field equations are equivalent to the effective Einstein equations of the $f\left(R\right)$ model in empty spacetime, but differ from them, as well as from standard General Relativity, in the presence of matter. Therefore, the predictions of these gravity theories could lead to some major differences, as compared to the predictions of standard General Relativity or to its extensions that ignore the role of matter, in several problems of current interest, such as cosmology, dark matter, dark energy, gravitational collapse or the generation of gravitational waves.

We have also reviewed, based on the already existing literature, the implications of the curvature-matter coupling models on the galactic dynamics, and we have shown that the behavior of the neutral hydrogen gas clouds outside of the galaxies and their flat rotation curves can be explained in terms of a non-minimal coupling between matter and curvature. We have also shown that in the curvature-matter coupling “dark matter” model, all of the relevant physical quantities, including the “dark mass” associated with the coupling, and which plays the role of dark matter, its corresponding density profile, as well as the curvature-matter coupling function can be expressed in terms of observable parameters: the tangential velocity, the baryonic (luminous) mass and the Doppler frequency shifts of test particles moving around the galaxy. Therefore, this opens the possibility of directly testing the modified gravity models with non-minimal coupling between matter and geometry by using direct astronomical and astrophysical observations at the galactic or extra-galactic scale. Since the observations on the galactic rotation curves are obtained from the Doppler frequency shifts, we have generalized the expression of the frequency shifts by including the effect of the curvature-matter coupling. Thus, at least in principle, the coupling function can be obtained directly from astronomical observations. We have also presented two other possibilities for testing these classes of gravity theories, at the level of the Solar System, by using the perihelion precession of Mercury, and at astrophysical scales, from possible observations of the gravitational tidal effects.

As future avenues of research, one should aim to characterize as much as possible the phenomenology predicted by these theories with curvature-matter couplings in order to find constraints arising from observations. The study of these phenomena may also provide some specific signatures and effects, which could distinguish and discriminate between the various theories of modified gravity. We also propose to use a background metric to analyse the dynamic system for specific curvature-matter coupling models and to use the data of SNIa (Type Ia supernovae), BAO and CMB shift parameter to obtain restrictions for the respective models and to explore in detail the analysis of structure formation.