Generalizing the coupling between geometry and matter: $f(R,L_m,T)$ gravity

We generalize and unify the $f(R,T)$ and $f(R,L_m)$ type gravity models by assuming that the gravitational Lagrangian is given by an arbitrary function of the Ricci scalar $R$, of the trace of the energy-momentum tensor $T$, and of the matter Lagrangian $L_m$, so that $L_{grav}=f(R,L_m,T)$. We obtain the gravitational field equations in the metric formalism, the equations of motion for test particles, and the energy and momentum balance equations, 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 equations of motion is also investigated, and the expression of the extra acceleration is obtained for small velocities and weak gravitational fields. The generalized Poisson equation is also obtained in the Newtonian limit, and the Dolgov-Kawasaki instability is also investigated. The cosmological implications of the theory are investigated for a homogeneous, isotropic and flat Universe for two particular choices of the Lagrangian density $f(R,L_m,T)$ of the gravitational field, with a multiplicative and additive algebraic structure in the matter couplings, respectively, and for two choices of the matter Lagrangian, by using both analytical and numerical methods.


I. INTRODUCTION
The unexpected discovery of the accelerating behavior of the Universe [1][2][3][4][5], that began at a redshift z of around z ≈ 0.6, brought into question the very foundations of the most successful gravitational theory developed so far, general relativity. The simplest answer to explain the late time acceleration is to resort again to the old cosmological constant Λ, which was introduced by Einstein to build the first general relativistic cosmological model [6]. Based on this ansatz, and strongly supported by the excellent fits a cosmological constant provide to the observations, the standard model of cosmology, called the ΛCDM paradigm, was established, which also requires the existence of another intriguing (and yet undetected) constituent of the Universe, called dark matter [7]. The ΛCDM model describes well observations [8][9][10][11], however, it is confronted with a major problem: no fundamental physical theory can give a reason for it. The main problem arises from the difficulties encountered in explaining the nature and origin of the cosmological constant [12][13][14]. A popular explanation for Λ as a Planck-scale vacuum energy density ρ vac leads to the "worst prediction in physics" [15], since ρ vac ≈ ( /c) k P l k dS k 2 + (mc/ ) 2 d 3 k ≈ ρ P l = c 5 / G 2 = 10 93 g/cm 3 , a predictions that differs by a factor of around 10 −120 from the observed value of the cosmological constant, ρ Λ = Λc 2 /8πG ≈ 10 −30 g/cm 3 [10].
Actually, the ΛCDM model can fit the cosmological observational data at a very high level of precision. Additionally, it is a very simple (empirical in some sense) general relativistic theoretical approach to cosmology. However, even at the observational level it faces some important problems, the most important being the "Hubble tension", representing the significant differences between the values of the Hubble constant, H 0 , as determined from the CMB measurement [11], and those obtained from direct measurements in the local Universe [16][17][18]. For example the SH0ES measurement of H 0 give H 0 = 74.03 ± 1.42 km/s/Mpc [16], while from the early Universe measurement of the Planck satellite one finds H 0 = 67.4 ± 0.5 km/s/Mpc [10], results which differ by ∼ 5σ.
On the other hand, important theoretical questions cannot be also answered in the framework of the ΛCDM paradigm. Why it happens that the numerical value of the cosmological constant is extremely small? Why is its value adjusted so fine? What determined the Universe to begin to accelerate only recently, after a long decelerating phase? And, after all, is a cosmological constant really necessary for theoretical and observational cosmology? Hence, the present day observational situation of cosmology strongly requires the exploration of alternative pathways for the understanding of the gravitational interaction, and its influence on the dynamics of the Universe. In fact, one can distinguish (at least) three possible theoretical approaches as substitutes of the ΛCDM paradigm.
The first avenue leading to an alternative description of cosmological evolution may be called the dark constituents model. It is based on the generalization of the Einstein gravitational field equations through the addition of two new terms in the total energy momentum tensor of the Universe, which describe dark energy and dark matter, respectively. In this approach gravitational phenomena are described by the generalized field equation [19] G µν = κ 2 T bar µν + κ 2 T DM µν (φ, ψ µ , ...) + κ 2 T DE µν (φ, ψ µ , ...), where G µν is the Einstein tensor, T bar µν , T DM µν (φ, ψ µ , ...), and T DE µν (φ, ψ µ , ...) denote the energy-momentum ten-sors of baryonic matter, dark matter and dark energy, respectively. In this approach dark energy is understood as a form of matter, an interpretation based on the equivalence between mass and energy in the theory of relativity. Both the energy-momentum tensors of dark matter and dark energy can be constructed from some scalar or vector fields. The simplest dark constituent model can be obtained by assuming that dark energy is a scalar field φ with a self-interaction potential V (φ). Hence the action of the gravitational theory be- 4 x, leading to so-called quintessence models [20][21][22][23][24][25]. Other dark constituent models are k-essence models [26][27][28], tachyon field models [29,30], models considering phantom fields [31][32][33], or quintom field models [34][35][36], respectively. In [25] it was pointed out that quintessence always lower H 0 relative to the ΛCDM model, thus making the Hubble tension worse. Chameleon field [37][38][39][40][41], Chaplygin gas [42,43] and vector field [44][45][46] dark energy models have also been investigated. Notwithstanding its impressive success, the dark constituents avenue is still confronted with its own intrinsic problems, the most important being the lack of any direct experimental evidence for the existence of dark matter. Furthermore, the existence of so many distinct field type models for dark energy, each one facing its own theoretical problems, leads to the problem of the existence of a single, theoretically consistent field theoretical picture of the major components of the Universe. For reviews of the dark energy models see [47][48][49][50].
A third way for understanding gravitation has also been developed, which we may call the dark coupling approach. The basic idea of this approach is the substitution of the standard Hilbert-Einstein Lagrangian density, which has a simple additive structure in curvature and matter Lagrangian, with a more general algebraic system. In this way the maximal extension of the standard Hilbert-Einstein gravitational Lagrangian can be established by considering that the gravitational action is an arbitrary analytical function of the curvature scalar R, of the matter Lagrangian L m , and, possibly, of other thermodynamic parameters. This approach leads naturally to the existence of a nonminimal coupling between geometry and matter. Many types of coupling between matter and geometry are possible, including the coupling between the Ricci scalar and the trace of the matter energymomentum tensor.
Therefore, in the dark coupling approach, the Einstein gravitational equations can be extended to the form with the effective energy-momentum tensor of the theory T (coup) µν (g µν , R, L m , T, R, T, ...) is obtained by considering generally the maximal extension of the Hilbert-Einstein Lagrangian, and by introducing a non-additive geometry-matter algebraic structure that describes the couplings between matter in all its forms, and spacetime curvature, respectively. It is important to mention that in both dark gravity and dark coupling theories the gravitational constant κ 2 turns into an effective quantity, generally a function of the couplings, and of the field parameters.
An alternative possibility of coupling between matter and geometry is described by the f (R, T ) gravity theory [105], in which the trace of the matter energy-momentum tensor is nonminimally coupled to geometry, with the action given by S = [f (R, T ) + L m ] √ −gd 4 x. The astrophysical and cosmological implications of f (R, T ) theory were investigated in detail in .
The viability of the f (R, T ) cosmological models were discussed in [117] for a Lagrangian density of the form f (R, T ) = R + λT + γ n T n , where γ n and n are constants. Several models were considered, corresponding to different values of the model parameters, and confronted with the H(z) data and Supernovae data. The H(z) comparison indicates consistency with ΛCDM at low redshifts, but a significant deviation of the models from the standard scenario does appear when going to higher redshifts. Similar pathological behaviors also show up in the dynamics of the deceleration parameter of the models, leading to universes where a transition from a decelerating phase to an accelerated one does not exist. The existence of such a transition is a basic requirement for any cosmological model compatible with the observational data. However, one can find models with transitions from deceleration to acceleration transition, but the predictions of these models are not in agreement with both the H(z) and Supernovae data at higher redshift regimes. On the other hand, at lower redshifts there is a good concordance between model and data. These results indicate the difficulties faced by the considered version of the f (R, T ) theory in correctly representing the cosmological evolution in its entirety.
A detailed analysis of the cosmological model based on the field equations R µν − (1/2)Rg µν = [(8π + λ)/λ]T µν + (p + T /2)g µν was performed in [121]. The maximum likelihood analysis was used to obtain the constraints on the Hubble parameter H 0 and the model parameter by taking into account the observational Hubble data set H(z), the Union 2.1 compilation data set SNeIa, the Baryon Acoustic Oscillation data BAO, and the joint data set H(z) + SNeIa and H(z) + SNeIa + BAO, respectively. It was found that the model is in good agreement with these observations.
An extension of the f (R, T ) gravity theory by introducing higher derivatives matter fields was proposed in [132], based on an action of the form S = (1/16π) f (R, T, T ) √ −gd 4 x + ǫ L m √ −gd 4 x, where ǫ = ±1. Accelerated evolution does appear in this scenario in the dust-filled Universe without any additional matter sources. The inflationary model is also in quite good agreement with observations. Gravitational theories in which the Ricci scalar is coupled to the square of the energy-momentum tensor T 2 = T µν T µν , with action of the form S = f R, T 2 √ −gd 4 x were analyzed in [133][134][135][136][137]. Theoretical models involving couplings between the Ricci tensor R µν and the energy momentum tensor T µν , extending the action of the standard f (R, T ) theory to an action of the form √ −gd 4 x were introduced and discussed in [138][139][140]. Derivative matter couplings are also considered in [141], where the Galileon interactions are promoted to contain matter Lagrangians. Nonminimal couplings between nonmetricity and matter were investigated in [142][143][144][145]. For reviews and detailed discussions on the gravitational theories with geometrymatter coupling see [146] and [147], respectively.
It is the goal of the present paper to extend, generalize, and at the same time unify two classes of gravita-tional theories with geometry-matter coupling, namely, the f (R, L m ) and the f (R, T ) theories, respectively. Hence we will consider a gravitational theory in which matter is non-minimally coupled with geometry at both the levels of the matter Lagrangian L m , and of the trace of the energy-momentum tensor T , with Lagrangian density given by f (R, L m , T ). This theory has as limiting cases the f (R), f (R, L m ) and f (R, T ) gravity theories, respectively.
After introducing the gravitational action, as a first step in our study, we obtain the gravitational field equations of the theory. The energy and momentum balance equations are also obtained, and, as usual in theories with geometry-matter coupling, it turns out the matter energy-momentum tensor is not conserved anymore. The Newtonian limit of the equations of motion is investigated in detail, and the expression of the extraacceleration is obtained for small particle velocities and weak gravitational fields. The Dolgov-Kawasaki instability is investigated in detail, and the condition for the existence of stable perturbations is obtained. The cosmological applications of the theory are also investigated. We consider two particular choices of the Lagrangian density f (R, L m , T ) of the field, having a multiplicative and additive algebraic structure in the matter couplings. We also investigate the effects of the two possible different choices of the matter Lagrangian (L m = p and L m = −ρ, respectively) on the description of the cosmological evolution.
The present paper is organized as follows. The gravitational field equations of the f (R, L m , T ) theory are obtained in Section II. The divergence of the matter energymomentum tensor is derived in Section III. The energy and momentum balance equations are also presented, and it is shown that the equation of motion of test particles is not geodetic, but it takes place in the presence of an extra-force. The Newtonian limit of the equation of motion is investigated by using a variational approach, and the general expression of the extra-acceleration is obtained in the limit of small velocities. The generalized Poisson equation for the Newtonian gravitational potential is also derived. The Dolgov-Kavasaki instability of the theory is analyzed in Section IV, and the condition for obtaining stable perturbations is found. The cosmological implications of the theory are examined in Section V, for two specific choices of the function f (R, L m , T ), which assumes a multiplicative and an additive form in the matter couplings. The cosmological evolution equations are studied for the two forms of the matter Lagrangian, L m = −ρ, and L m = p, respectively. The possibility of the existence of a de Sitter type phase for the models is explored. We discuss and conclude our results in Section VI.
We assume that the action describing the modified gravity theory with strong geometry-matter coupling takes the form where f (R, L m , T ) is an arbitrary function of the Ricci scalar R, of the trace T of the matter energy-momentum tensor T µν , T = g µν T µν , and of the Lagrangian density of the ordinary matter, L m , respectively. The main physical motivation for considering this function is that it will ensure that the modification is non-trivial for all kinds of matter fields, including radiation with T = 0. The energy-momentum tensor of the ordinary (baryonic) matter is defined in the standard way according to the relation [148] In the following we introduce the important assumption that the Lagrangian density L m of the matter depends only on the metric tensor components g µν , and not on its derivatives. Therefore for T µν we find the simple expression The variation of the action S with respect to the metric tensor components g µν can be obtained as where we have denoted f R = ∂f /∂R, f T = ∂f /∂T , and f L = ∂f /∂L m , respectively. The variation of the Ricci scalar can be obtained from the Palatini identity δR σν = ∇ ρ (δΓ ρ νσ ) − ∇ ν (δΓ ρ ρσ ) [149] as where ∇ λ represents the covariant derivative with respect to the connection Γ, associated to the metric g, while Γ λ µν denotes the Christoffel symbols associated to the metric. We also assume the metric compatibility of the covariant derivative, ∇ σ g µν = 0. The variation of the trace T = T µ µ of the matter energy-momentum tensor T µν with respect to the metric tensor is defined as where and we have denoted After substituting the above expressions into Eq. (7), and integrating by parts, one can obtain the field equation of the f (R, T, L m ) gravity as For perfect fluid matter Lagrangian with L m = −ρ, and also for the scalar field theory with L m = − 1 2 ∂ µ φ∂ µ φ+V (φ), one can immediately see that τ µν vanishes. However, in the case of vector field theory, with If f = R (the Hilbert-Einstein Lagrangian), from Eq. (12) we recover the standard field equations of general relativity, R µν − (1/2)g µν R = 8πT µν [148]. If f = f (R, L m ), we reobtain the field equations of the f (R, L m ) gravity model [91], while f = f (R, T ) gives the field equations of the f (R, T ) theory [105].
By taking the trace of the gravitational field equations (12), in the absence of electromagnetic type fields we obtain With the use of Eq. (14) the gravitational field equation of the f (R, T, L m ) gravity theory can be written in an alternative form as We call this form of the gravitational field equations of the f (R, T, L m ) the traceless representation of the theory. It is interesting to note that for f (R) = R, Eq. (15), the basic field equation of the traceless representation of the f (R, L m , T )) reduces to the geometry-matter symmetric Einstein equations [150,151], which have been proposed as a possible solution to the cosmological constant problem. For the applications and further generalizations of the geometry-matter symmetric Einstein equations see [46].

III. ENERGY AND MOMENTUM BALANCE EQUATIONS, EQUATIONS OF MOTION, AND EXTRA-FORCE IN f (R, Lm, T ) GRAVITY
It is obvious that due of the non-minimal coupling between matter and geometry, in the present theory the energy-momentum tensor is not conserved. In this Section we will first obtain the non-conservation equation of the matter energy-momentum tensor (the energy and momentum balance equations), and then we will consider the equations of motion of particles, thus obtaining the expression of the extra-force generated by the curvaturematter coupling.
A. The divergence of the matter energy-momentum tensor Let us start with taking the covariant derivation of the gravitational field equations Eq. (12) where we have defined One should note that in the case of perfect fluid and scalar field theory we have A ν = 0. By taking into account that (19) and after using the geometric identities ∇ µ G µν = 0, and [152] the non-conservation equation of the energy-momentum tensor is reduced to where we have defined Eq. (21) is direct a consequence of the presence of the matter fields in the expression of the gravitational Lagrange density, given by the function f (R, L m , T ). One can immediately see that in the case of f T = 0 = f L , the matter content of the Universe is conserved. By using the traceless representation of the field equations, one can obtain the non-conservation of the energymomentum tensor as where we have denoted

B. The energy and momentum balance equations
In the following we adopt the simplifying assumption according to which the matter content of the Universe consists of a perfect fluid that can be described by only two basic thermodynamics parameters, the energy density ρ, and the thermodynamic pressure p of the fluid, respectively. In this case the matter energy-momentum tensor is given by where the four-velocity u µ is normalized as u µ u µ = −1.
Next, we introduce the projection operator h ν λ , defined as with the basic property u ν h ν λ ≡ 0. By taking the covariant divergence of Eq. (25), we first find Therefore Eq. (21) takes the form where we have taken into account that A ν = 0 in the case of a perfect fluid. By multiplying Eq. (28) with u ν , and by considering the identity u ν ∇ µ u ν = 0, we arrive at the energy balance equation in the f (R, L m , T ) gravity theory as given bẏ where we have denoted H = (1/3)∇ µ u µ , and˙= u µ ∇ µ = d/ds, respectively. Here ds denotes the line element constructed with the help of the metric g µν , ds 2 = g µν dx µ dx ν . We note that the equality u µ ∇ µ = d/ds is only true when we apply it on a scalar function. In the case of tensors, we will substitute the definition with the absolute derivative u µ ∇ µ = D/Ds, as we will see in the following.
The multiplication of Eq. (28) with h ν λ gives the momentum balance equation for a perfect fluid in the alternative f (R, L m , T ) gravity theory as where we have defined D λ ≡ h νλ ∇ ν , with D λ orthogonal to the vector field u ν , and lying in the spatial hypersurfaces of constant time.
Alternatively, we can write down the equation of motion of particles as where or, equivalently, is the extra-force generated by the geometry-matter coupling in the f (R, L m , T ) theory. It is straightforward to see that the extra-force is perpendicular to the fourvelocity, u λ f λ = 0.
C. The variational principle for the equation of motion of a test particle In the following we assume that the extra-force in the f (R, L m , T ) theory, as given by Eq. (32), can be represented in a formal mathematical way as it happens that the extra-force admits such a mathematical representation, the equation of motion of the massive particles can be derived from a variational principle, which is given by where S p and L p = √ Q √ −g µν u µ u ν are the action and the Lagrangian density for massive test bodies moving in a curved geometry described by the metric g µν . As usual, by s we have denoted the affine parameter along the curve.
To prove this result we begin by writing down the Euler-Lagrange equations obtained from the variation of the action (35), and which are given by d ds By considering the mathematical relations and respectively, after a simple calculation we obtain the equations of motion of the particle in the form By the simple identification of the terms in the equation of motion given by Eq. (39), and the expression of the extra-force generated by the geometry-matter coupling, as given by Eq. (34), we can immediately obtain the explicit functional form of √ Q. In the limit √ Q → 1, we reobtain the usual general relativistic equation d 2 x µ ds 2 + Γ µ νλ u ν u λ = 0, describing geodesic motion in the gravitational field.

D. The Newtonian limit of the equation of motion
The variational principle (35) provides a very simple and efficient way for the study of the Newtonian limit of the equations of motion of the test particles in f (R, T, L m ) gravity. In the Newtonian limit of the weak gravitational fields and slow motions, one can approximate the metric according to where by φ we have denoted the Newtonian potential, and v is the usual tridimensional velocity of the fluid, satisfying the condition v 2 /c 2 ≪ 1. By assuming that in the Newtonian limit of weak gravitational fields the function √ Q can be represented in the f (R, L m , T ) theory as, where we assume that in the first order of approximation U (φ, ρ, p, T, L m , f m , ∇ ν T, ...) ≪ 1, the equations of motion of the fluid do follow from the variational principle and are obtained as where a denotes the total acceleration of the system, the term a N = −∇φ gives the Newtonian gravitational acceleration, while represents the supplementary acceleration induced by the inclusion of the different forms of the geometry-matter coupling in the action of the gravitational field of the f (R, L m , T ) theory.
In the case of the extra-force generated in the f (R, L m , T ) gravity, the expression of Q can be obtained from the relation The function √ Q can always be represented in a closed form obtained by integrating the left-hand side of Eq. (45). But we must point out that generally Q cannot be expressed in a simple mathematical form, and to find its expression some approximate methods must be used.
In order to obtain the explicit expressions of the extraacceleration induced by the geometry-matter coupling, we analyze the simple case in which the pressure of the gravitating fluid is given as a function of the density by a linear barotropic equation of state having the standard form p = wρ, where w is a constant satisfying the condition w ≪ 1. Hence it follows that we can approximate ρ + p ≈ ρ, and T = −ρ + 3p ≈ −ρ, respectively. Hence in this approximation T is a function of the density only. Furthermore, we assume that the function f (R, T, L m ) can be represented as f (R, T, L m ) = R + ǫΨ (L m , T ), where ǫ is a small parameter. We also adopt for the matter Lagrangian the expression L m = −ρ. Hence the gravitational Lagrangian density becomes a function of the Ricci scalar and of the matter density only, f (R, T, L m ) = R + ǫΨ (ρ). Now we can expand Ψ(ρ) near a fixed value ρ 0 of the matter density, thus obtaining in the first order of approximation where we have denoted a 0 = Ψ (ρ 0 ) and b 0 = [dΨ(ρ)/dρ] | ρ=ρ0 , respectively. Hence we immediately find f T = ǫb 0 , f L = −ǫb 0 , and f m = f T + f L /2 = ǫb 0 /2, respectively. Then, with the use of the above approximations, Eq. (45) giving the expression of √ Q takes the form giving where α = w − (ǫb 0 2) / (8π + ǫb 0 /2). Eq. (48) also maintains its validity for a fluid satisfying a linear barotropic equation of state of the form p = (γ − 1) ρ, γ = constant, with (γ − 1) not necessarily very small. Hence an equation of the form (48) can be applied for the description of both non-relativistic and the extreme relativistic regimes. By using the mathematical relation x α = exp(α ln x) ≈ 1 + α ln x, we can approximate √ Q (ρ) given by Eq. (48) as where Hence in the first order of approximation in f (R, L m , T ) gravity the equations of motion of the gravitating hydrodynamic type system can be obtained from the variational principle and they are given by with the extra-acceleration of the fluid given by In the present approach and within the framework of the adopted approximations the supplementary acceleration induced by the modification of the action of the gravitational field due to the geometry-matter coupling is a function of the matter density only.

E. The modified Poisson equation
In this Subsection we will obtain the modified Poisson equation in f (R, L m , T ) gravity. In order to find it we consider the first order approximation of the function f (R, L m , T ), which can be written with the help of a Taylor series expansion in the vicinity of a fixed point where we have denoted Therefore, in this approximation, the field equations (12) for a perfect fluid with energy-momentum tensor given by Eq. (25) become where we have denoted Equivalently, Eq. (58) can be reformulated as In the Newtonian limit one can write the metric as where φ is the Newtonian potential. In the following we assume that the macroscopic motion takes place at small velocities, and consequently in the four-velocity we can keep only the time component, and neglect all the space components with a very good approximation. Therefore we have u 0 ≈ −1, and u i ≈ 0, i = 1, 2, 3. We will apply the same approximation in the gravitational field equations, thus keeping only their time component. A similar approximation is used for the components of the energy-momentum tensor. Therefore, we can approximate T µν as T 0 0 ≈ −ρ, and T = −ρ, respectively. By taking into account the explicit expression of the metric tensor, we obtain Γ α 00 ≈ −(1/2)g αβ ∂g 00 /∂x β = ∂φ/∂x α , giving R 0 0 = −∆φ. Then Eq. (60) immediately gives ∆φ = 4πG N ew ρ + Λ, (62) with G N ew = G ef f − β/4π. Hence, we have obtained the generalized Poisson equation in f (R, L m , T ) gravity theory, which, in the first order of approximation, determines a slight modification of the gravitational constant, and induces a cosmological term.

IV. THE DOLGOV-KAWASAKI INSTABILITY
In this Section we will obtain the Dolgov-Kawasaki criterion [153] in f (R, L m , T ) gravity theory. Let's consider the dynamics of a small celestial body like, for example, the Sun. In this case, the gravitational force is very small, and the metric is very close to the flat one. However, the curvature of the space-time is still non-zero, and the geometry is curved. Suppose that there exist a solution of the field equations corresponding to constant values of the variables of the theory, R 0 , T 0 and L 0 , respectively. Let us now perturb the curvature, the matter energymomentum tensor, and the matter Lagrangian around their background values as Also, we can expand the tensor τ µν as τ µν = τ 0µν + τ 1µν . In order to be consistent with the Solar System observations, in addition to the above relations, we assume that the function f can be expressed as where ǫ is a small constant, with its value set by Solar System tests, and Φ (R, L m , T ) is an arbitrary function. Let us start with the trace of Einstein's field equation (14). To the first order in perturbations, we have f (R, L m , T ) = R 0 + ǫΦ(0) + (1 + ǫΦ R (0))R 1 + H 1 , (67) where H 1 contains all first order contributions of the matter field, i.e. T 1 and L 1 . One can also obtain The higher order derivatives of the function f can be easily written as above. We also have (71) By using the above approximations, it follows that the d'Alembertian can be expanded in the Minkowski space as The resulting trace equation (14) at first order can be written asR where we have defined the effective mass, and effective potential, respectively, as and In order to have a stable perturbation, and to avoid tachyonic instability, one must impose the condition m 2 ef f ≥ 0. However, the mass term (74) is dominated by its first term, which gives us the constraint Φ RR (0) ≥ 0. Noting that Φ RR (0) = f RR (0), we finally obtain the Dolgov-Kawasaki criterion in f (R, L m , T ) gravity as The same condition applies in the special cases of f (R), f (R, T ) and f (R, L m ) theories, which were proven elsewhere.

V. COSMOLOGICAL CONSEQUENCES OF f (R, Lm, T ) GRAVITY
Let us assume that the Universe can be described by the flat, homogenous and isotropic Friedmann-Lemaitre-Robertson-Walker (FLRW) metric of the form where a(t) is the scale factor. In the following we denote by H =ȧ/a the Hubble function, describing the rate of expansion of the Universe. The accelerating/decelerating nature of the cosmological expansion can be characterized with the help of the deceleration parameter q, defined as Also, let us assume that the Universe is filled with a perfect fluid. We adopt for the Lagrangian density of the cosmic (baryonic) matter the expression L m = −ρ.
In the comoving frame the non-zero components of the energy-momentum tensor are given by Moreover, we impose on the function f the condition f R = 0, a condition that is assumed to be valid for all times.

A. Generalized Friedmann equations
In the case of the traceless representation of the f (R, L m , T ) gravitational theory, the traceless equation (15) and the trace equation (14) reduce for the FLRW metric to the evolution equations (generalized Friedmann equations) andḢ In order to obtain a consistent solution it is usually necessary to solve both cosmological equations above. This is due to the fact that the original set of Einstein equations are decomposed into a traceless equation, and the trace equation, with L m appearing only in the second equation. To close the system of field equations we must add to cosmological system the energy balance equation (23), which gives the time evolution of the density aṡ The deceleration parameter can be obtained in a general form as Eq. (80) immediately gives the main criterion for the existence of an asymptotic de Sitter type evolution with H = constant = H 0 in f (R, L m , T ) gravity as Assuming that matter and geometry are minimally coupled f (R, L m , T ) = R + g (L m , T ), the condition for accelerated de Sitter type expansion reduces to representing a linear first order partial differential equation, with the general solution given by where α is an arbitrary constant and C is an arbitrary function of the variable (2L m − T ) /2. In the simple case of an additive structure of the gravitational Lagrangian of the form f (R, L m , T ) = R + αL m + βT , the condition for de Sitter type expansion is given by 8π + α/2 + β = 0.
The general condition for an accelerating expansion of the Universe, with q < 0, can be formulated as The condition given by Eq. (87) for the Lagrangian density.

The radiation dominated Universe
In the case of radiation dominated Universe, we have p = ρ/3, and the evolution equations (80), (81) and (82) can be simplified to and 2πρ − (αρ − 8π)Hρ = 0, substituting the above equation into the field equations, we have where A = 128π 2 − αβ. Also, in the case α > 0, we have an accelerating solution with the Hubble parameter with energy density ρ = 8π/α.

The dust Universe
In the case of dust dominated Universe, one has the evolution equations and respectively. By using Eqs. (94) and (95), we obtain Substituting the above relation for the Hubble function in the field equation (95) we find where we have defined This model has also a de Sitter type expansionary phase with H = H 0 = constant, for H 0 given by Hence in this model the de Sitter type phase is triggered by the generalized geometry-matter coupling, and takes place in the presence of ordinary matter only.

Numerical results
Now, we consider the cosmological evolution of the Universe in the present model by assuming it is filled with a mixture of dust and radiation, with the total energy density given by The equation of state of the radiation is p = p r = ρ r /3. In order to simplify the mathematical formalism we introduce a set of dimensionless parameters defined as where H 0 is the current value of the Hubble function.
In terms of the redshift the equations of motion (80) and (81) take the form and respectively, where the prime denotes the derivative with respect to the redshift. One can write the energy balance equation (82) as two coupled generalized conservation equations for the dust and the radiation separately as and respectively. It should be mentioned that the sum of above equations is equal to the original conservation equation (82), and are written in a way that in the limitᾱ = 0, they reduce to the standard conservation equations for dust and radiation.
For the initial (present day) values of the cosmological parameters we adopt the values h(0) = 1, Ω m0 = 0.305 and Ω r0 = 0.53 × 10 −4 , respectively. Evaluating Eqs. (102) and (103) at the present time (z = 0) one can find the numerical values of the constantβ in terms ofᾱ asβ The redshift dependence of the deceleration parameter is given by Using Eqs. (102) and (103) one can obtain the present value of the deceleration parameter q 0 as The behaviors of the Hubble function and of the deceleration parameter are shown in Fig. 1 for three different values ofᾱ. As one can see from the left panel of Fig. 1, the Hubble function is a monotonically increasing function of z, and up to the redshift z ≈ 2, its evolution is weakly dependent of the numerical values of the model parameterᾱ. For larger values of the redshift, significant differences between the model prediction and the observations do appear for someᾱ values. However,ᾱ = 0.0007 reproduces well the observed behavior of H(z) up to at least a redshift of z ≈ 4. In order to compare the present model with the ΛCDM theory, let us see the differences (in %) between the two models at some specific redshifts. Forᾱ = −0.002, h differs from the ΛCDM model at redshifts z = 1 and z = 3 by 0.35%, and 3.8% respectively. Also forᾱ = 0.007 for the deviations of h from the ΛCDM theory at z = 1 and z = 3 one obtains 0.12% and 1.30%, respectively.
The variation of the deceleration parameter of the model, is represented in the right panel of Fig. 1. The Universe enters into an accelerating expansionary phase at z ≈ 0.5, and in the range z ∈ (0, 1) the cosmological evolution is weakly dependent of the numerical values of the model parameters. However, for larger redshifts, the variation of q shows large deviations from the predictions of the ΛCDM model, withᾱ = 0.0007 giving the closest description to the predictions of standard general relativistic cosmology. It is well-known that there are two choices for the perfect fluid matter Lagrangian, i.e. L m = −ρ and L m = p. These two choices lead to the same matter energy-momentum tensor in the case of Einstein's gravity. However, in the case of modified theories of gravity, these two distinct possibilities lead to different dynamical evolutions of the Universe.
As in the case of L m = −ρ, the above set of equations admit the radiation dominated solution with and where we have defined D = 8αρ 2 + 144πρ − 9β.

The dust Universe
In the case of matter dominated Universe with p = 0, the cosmological evolution equations can be written as and respectively. By subtracting Eqs. (112) and (113) we obtain H(t) as Now, by using the expression of the Hubble function in terms of the energy density, and using the field equations, we obtain For the Hubble parameter, we have also plotted the observational data together with their errors [170]. The ΛCDM curve is depicted as a red solid curve.

Numerical results
To describe the general evolution of the Universe we assume again that ρ = ρ m + ρ r . The generalized Friedmann equations in terms of redshift are and respectively. The generalized conservation Eq.(82) for dust and radiation takes the form and respectively. For this case, with the use of Eqs. (117) and (118), the constant parameterβ can be obtained as The current value for the deceleration parameter is given by The variations of the Hubble function and of the deceleration parameter are shown, as functions of the redshift, in Figs. 2.
As one can see from the left panel of Fig. 2, up to a redshift of z ≈ 1, the present model gives a good description of the observational data for h(z), and its predictions almost coincide with the ΛCDM model. The dependence on the model parameterᾱ is weak in the low redshift range. However, for z > 1, the model behavior is influenced by the numerical values ofᾱ, but a good description of the observational data for the Hubble function up to z ≈ 4 can be obtained forᾱ = −0.13. Now, we can obtain the deviations (in %) of the dimensionless Hubble parameter h from the ΛCDM values. At z = 1 forᾱ = −0.13 andᾱ = −0.45, we have −1.34%, and −1.91%, respectively. However, at z = 3 the deviations become more significant, and can be obtained as −4.57% and −6.54% respectively.
A comparable situation can be seen in the description of the deceleration parameter, shown in the right panel of Fig. 2, whose evolution is strongly dependent on the model parameters even at low redshifts. However,ᾱ = −0.13 gives a very good description of the ΛCDM model values, with the differences slightly increasing in the low redshift region. Hence, forᾱ = −0.13, the present model gives a good description of the observational data for the expansion rate, and of the deceleration parameter up to a redshift of z ≈ 4.  [170]. The ΛCDM curve is depicted as a red solid curve.
algebraic structure of the form where k(R), g(T ) and h (L m ) are arbitrary functions of their arguments. With this choice of f the cosmological evolution equations can be written as and (8π + βg ′ + γh ′ )(ρ +ṗ + 4H(ρ + p)) respectively. In the following we will adopt for the additive gravitational Lagrange density the simple form corresponding to k(R) = 0, g(T ) = T , and h (L m ) = L m , so that As for the choice of the matter Lagrangian, we will restrict our analysis to the case L m = −ρ only.

The radiation dominated Universe
By assuming that the matter content of the Universe is radiation with p = ρ/3, the generalized Friedmann equations take the form while the evolution of the energy density can be obtained from the energy balance equation as (γ + 8π)ρ + 4(β + γ + 8π)Hρ = 0.
One should note that the energy balance equation is not independent of the field equations (128) and (129). From Eqs. (128) and (129) one can obtain the time evolution of the Hubble function and of the energy density as follows and where and c 1 is a constant of integration. By using Eq. (131) one can obtain the scale factor as where c 2 is a constant of integration.
These equations have the following exact solutions for the Hubble parameter and matter energy density, and ρ = λ β + 2γ + 16π sec 2 (g 2 (t)), respectively, where and c 1 is an integration constant. Now one can obtain the scale factor as a = c 2 [cos(g 2 (t))] β+2(γ+8π) where c 2 is an integration constant.

Numerical solutions
In this Section we consider the numerical solutions for the additive choice of the Lagrangian in the presence of dust and radiation, as the matter constituents of the Universe. We rescale the parameters that enter in the gravitational Lagrangian and the cosmological equations according toβ In terms of the redshift the generalized cosmological equations are and respectively, while the conservation equations for dust and radiation are One should note that the value of the parameterλ can be obtained by the use of Eqs. (143) and (144) at z = 0 in terms ofβ andγ as λ = −2 + 2 1 +β + 2γ Ω m0 + 2 (1 + 2γ) Ω r0 . (146) Also the current value of the deceleration parameter in this case is The results obtained by the numerical integration of the generalized Friedmann equations for this case are shown, for four different choices ofβ andγ, in Figs. 3.
As one can see from the left panel of Fig. 3, for some model parameter values the simple additive algebraic structure of the f (R, L m , T ) gravity theory provides a good description of the behavior of the Hubble function in the redshift range z < 1. For higher redshifts a strong dependence on the model parameters does appear, and the differences with respect to the standard ΛCDM model become important.
For the case withβ = −0.12 andγ = 0.15, the percent deviation of the dimensionless Hubble parameter h with respect to the ΛCDM values at z = 1 and z = 3 are 2.66%, and 11.80%, respectively. One can also compare these values with the corresponding quantities for the parametersβ = 0.11 andγ = −0.32, which are 6.80%, and −11.80%, respectively.
An analogues strong dependence on the numerical values ofβ andγ can be seen in the behavior of the deceleration parameter, as shown in the left panel of Fig. 3, with some combinations of the parameters providing a good description of the ΛCDM data at low redshifts, and with other combinations describing better the standard cosmological model at higher redshifts.

VI. DISCUSSIONS AND FINAL REMARKS
In the present paper we have considered a generalized gravity model with the gravitational Lagrangian density given by an arbitrary function of the Ricci scalar, of the trace of the energy-momentum tensor T , and of the Lagrange density of the matter L m , respectively. The proposed action represents a non-trivial extension of the standard Hilbert-Einstein action for the standard general relativistic gravitational field, S = [R/2 + L m ] √ −g d 4 x, and it unifies two already welldeveloped modified gravity theories, the f (R, L m ) theory, and the f (R, T ) theory, respectively. It is important to note that the classical Hilbert-Einstein action, as  [170]. The ΛCDM curve is depicted as a red solid curve.
well as most of its generalizations, has an intrinsic additive algebraic structure, and is constructed as the sum of independent geometrical and physical terms, without involving any coupling between matter and geometry.
On the other hand, the f (R, L m , T ) modified theory of gravity considered in the present paper offers the possibility of going beyond this simple algebraic structure, and also brings new degrees of freedom into the gravitational theories. Moreover, besides its theoretical interest, such a generalization may cure some of the cosmological pathologies of either f (R, T ) or f (R, L m ) theories, like, for example, those pointed out in [117], by providing a description of the cosmological evolution valid from high to low redshifts. In our study we have restricted our numerical analysis of the cosmological evolution for a redshift range 0 ≤ z ≤ 3, and at higher redshifts some differences between the considered models do appear. To clarify the nature of these deviations it is necessary to perform a fit of the model predictions and observational data. But our preliminary results already indicate that the theory can provide an acceptable description to the observational data, also representing an attractive alternative to the ΛCDM cosmology.
When formulating the f (R, L m , T ) gravity theory an interesting analogy with the geometry-matter symmetric Einstein theory is found. It is known that the field equations (16) contain the cosmological constant as a constant of integration [154][155][156][157][158]. This can be shown immediately by taking the divergence of Eqs. (16), together with the postulates of the conservation of the matter energy-momentum tensor, ∇ µ T µν = 0, which gives ∇ µ R = −8π∇ µ T , ∇ µ (R + 8π) = 0, and R + 8πT = Λ, where Λ is an integration constant. After substitution of R = Λ − 8πT into Eq. (16) we reobtain the standard Einstein equations in the presence of a cosmological constant, which thus becomes a simple integration con-stant. This representation of the Einstein field equations is also called unimodular gravity. The trace-free Einstein equations can be obtained from the variational principle [159,160] The field equations of the traceless representation of the f (R, L m , T ) gravity theory represent an extension of the geometry-matter symmetric Einstein theory, with the field equations taking the form and respectively. In the absence of the electromagnetic type interactions, the supplementary term U µν generated by the geometry-matter coupling takes the simple form As one can see from the field equations (149), the matter energy-momentum tensor T µν is generally not conserved, and ∇ µ T µν = 0. Of course one can construct conservative models in which the energy-momentum conservation is imposed a priori, ∇ µ T µν = 0, but the possible non-conservation of T µν has far reaching physical implications. First of all, the motion of free particles is not geodesic anymore, and it takes place in the presence of an extra-force generated by the geometry-matter coupling. The extra-force generates an extra-acceleration, and hence, in the Newtonian limit, the total acceleration a of a massive gravitating objects can be written as a = a N + a E . If a E = 0, then the acceleration of the object is a = a N , with the Newtonian acceleration of an object moving in the gravitational field created by a mass M given by a N = −GM r/r 3 . From the acceleration equation we obtain which allows to express the vector a N as where C is an arbitrary vector, which for simplicity will be taken as zero in the following. From Eq. (153) it follows that in order to obtain a consistent mathematical description a E · a = 0, that is, the vectors a E and a cannot be orthogonal. Without any loss of generality we will assume that they are parallel. By denoting we obtain that is, the total acceleration of an object in the gravitational field of a mass M is proportional to the Newtonian acceleration. Such a relation has been already inferred from the study of the galactic rotation curves of massive test particles gravitating around the galactic center, and it is called the radial acceleration relation (RAR) [161][162][163][164]. The radial acceleration relation is an empirical relation that tries to establish a connection between the centripetal acceleration a obs (R) = V 2 rot (R)/R observed in galaxies, and usually attributed to the presence of dark matter, and the Newtonian acceleration a bar (R) = V 2 bar /R of the observed baryonic matter distribution. The relation can be formulated as where ν(x) denotes a fitting function whose mathematical form must be determined observationally, while a + denotes an acceleration scale. A similar relation, indicating the proportionality of the total acceleration to the Newtonian one does naturally appear in the low velocity/low density limit of the f (R, L m , T ) gravity, as shown by Eq. (155), where such a relation is a direct consequence of the geometry-matter coupling. More interestingly, the extra-acceleration predicted by the theory in the Newtonian limit is dependent on the density of the baryonic matter, a E = a E (ρ), a result that indicates a possible relation with the chameleon field models [165][166][167]. The chameleon scalar field is a light scalar field, with a mass depending on the background baryonic matter density. However, such environmental effects can also be generated by the geometry-matter coupling, indicating a baryonic density dependence of the astrophysical processes on large scales. We have also obtained the generalized Poisson equation in the Newtonian limit of the theory, by using a series expansion near a fixed point of the theory, and restricting our analysis to the first order terms. The series expansion also generates a constant term, that can be interpreted as the cosmological constant, and which in the present approach is nothing but the difference of the Lagrangian density at a fixed point at its first order approxima- Hence in this approximation the cosmological constant appears as a second order effect in the series expansion of the full Lagrangian density f (R, L m , T ), and hence this interpretation may justify its smallness. In the same order of approximation there is a shift of the gravitational constant, determined by the derivatives of f with respect to L m and T estimated for some background values of the model parameters. However, the constraints on the numerical value of the gravitational constant [168,169] could lead to the determination of some limits of the gravitational Lagrangian in modified gravity with geometrymatter coupling.
The Dolgov-Kawasaki instability is a general feature of modified theories of gravity, in which higher order equations of motion do emerge that may exhibit unphysical behaviors. For example, unitarity may be broken, and ghosts or tachyons could develop. Moreover, the breaking of the laws of gravity in the weak field regime and for small curvatures leads to contradiction with basic observational facts. We have investigated in detail the Dolgov-Kawasaki instability for the f (R, L m , T ) theory, and we have obtained the criterion for the stability of the theory, f RR (0) ≥ 0. This is the standard stability condition as obtained in, for example, f (R) gravity, a condition that assures the instability of the perturbations, and the absence of the tachyonic instabilities.
An important field of application of f (R, L m , T ) is cosmology. We have obtained the generalized Friedmann equations (80) and (81), which must be usually considered together with the energy-momentum balance equation (82). The condition for the accelerated expansion is given by Eq. (87), and it can be satisfied by a large variety of functional forms of the Lagrange density of the gravitational field in the f (R, L m , T ) gravity. We have investigated in detail two distinct classes of cosmological models, corresponding to a multiplicative and an additive structure of the terms involving the geometrymatter coupling. The evolutionary cosmological models have a large variety of behaviors, and different types of accelerating or decelerating models can be constructed. Overall, for some values of the model parameters, the considered model give a good description of the cosmo-logical observational data, and they can reproduce the predictions of the ΛCDM paradigm at both low (z < 1) and higher (z ≈ 4) redshifts.
It should be noted here that in the plots of the Hubble functions H(z), we have normalized them to the current value of the Hubble function H 0 . However, for every choice of the model parameters, the current value of the Hubble function is different, and this value should be obtained by fitting the theory with the observational data, thus obtaining some observational constraints on the model parameters. This will be done in future works. As for the Hubble tension problem, we would like to point out that all modified gravity theories that make the Hubble tension better, will make the σ 0 8 tension worsen [171][172][173]. However, as a possible solution to these problems, in [174] it was shown that a decaying dark matter model can be considered as a solution to both tensions simultaneously. Such an approach can also be extended to gravitational theories with geometry-matter coupling.
The present gravitational theory still faces the problem of the degeneracy of the matter Lagrangian, leading to two different classes of models corresponding to the choices L m = −ρ, and L m = p, respectively.
However, a possible solution of this problem can be obtained as follows.
For simplicity we assume that the Lagrangian density L m of the ordinary matter is a function of the metric tensor components g µν only, and not of its derivatives. Then for the energymomentum tensor we obtain the expression T µν = L m g µν − 2∂L m /∂g µν . Moreover, we assume that the matter satisfies a barotropic equation of state p = p(ρ), and thus the matter Lagrangian, that generally can be considered a function of both ρ and p, L m = L m (ρ, p) becomes a function of the matter density only, L m = L m (ρ). By assuming that the matter current is conserved, ∇ ν (ρu ν ) = 0, and by taking into account the relation δρ = − (1/2) ρ (g µν + u µ u ν ) δg µν [94], it follows that the energy-momentum tensor can be written in the form With the use of the mathematical identity u ν ∇ ν u µ = d 2 x µ /ds 2 + Γ µ νλ u ν u λ , from the condition ∇ µ T µν = 0 we obtain the equation of motion of the fluid as d 2 x µ ds 2 + Γ µ νλ u ν u λ + h µν ∇ ν ln K dL m (ρ) dρ = 0, (158) where K is a constant of integration. A comparison of Eqs. (158) and (39) gives the relation K dL m (ρ) dρ = Q ≈ 1 + U (φ, ρ, p, T, L m , f m , ∇ ν T, ...) , where we have also used Eq. (41). The above equation leads us to the general result that in modified gravity theories with geometry-matter coupling the matter energy-momentum tensor must be constructed dynamically for each individual model, by solving the differential equation (159). In the first order of approximation, when U (φ, ρ, p, T, L m , f m , ∇ ν T, ...) ≪ 1, Eq. (159) gives L m ≈ ρ/K, which is the (approximate) expression we have adopted in our investigations of the cosmological models, and which coincides with the standard general relativistic matter Lagrangian. This limiting procedure also fixes the value of the constant K as K = −1.
To conclude, the current investigation of the unified f (R, T, L m ) gravity theory opens new avenues for the study of the theoretical, observational and even experimental aspects of the large class of alternative theories of gravity in the presence of geometry-matter coupling. In addition, the basic results obtained in the present work may also contribute to the better understanding of the physical and mathematical structure of other modified gravity theories.