f(R,L_m) Gravity Insights

• f(R,L_m) gravity is a framework where the Lagrangian is an arbitrary function of the Ricci scalar and matter, modifying general relativity through curvature–matter couplings.
• The theory yields altered field equations, extra forces causing non-geodesic motion, and modified conservation laws that address cosmic acceleration and structure formation.
• It applies to diverse phenomena from late-time cosmology and compact star equilibrium to wormhole constructions, while raising challenges in stability and energy–momentum conservation.

$f(R,L_m)$ gravity is a metric modified-gravity framework in which the gravitational Lagrangian is taken to be an arbitrary function of the Ricci scalar $R$ and the matter Lagrangian density $L_m$, rather than a sum of a purely geometric Einstein–Hilbert term and a minimally coupled matter term. In its basic form, the theory is defined by

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

and reduces to standard General Relativity (GR) for $f(R,L_m)=R/2+L_m$ in units $8\pi G=c=1$ (Harko et al., 2010). The explicit curvature–matter dependence generically modifies the field equations, changes the covariant balance laws, and can induce non-geodesic motion through an extra force orthogonal to the four-velocity. In the literature represented here, the framework has been applied to late-time cosmology, dynamical-systems analyses, perturbative growth of structure, bouncing and baryogenic early-universe scenarios, white-dwarf equilibrium, and several wormhole constructions (Harko et al., 2010).

1. General action, field equations, and formal structure

The general metric variation of the action yields the field equations

$$f_R\,R_{\mu\nu}+\bigl(g_{\mu\nu}\Box-\nabla_\mu\nabla_\nu\bigr)f_R -\tfrac12\bigl(f-f_{L_m}L_m\bigr)g_{\mu\nu} =\tfrac12\,f_{L_m}\,T_{\mu\nu},$$

with

$$f_R\equiv \frac{\partial f}{\partial R},\qquad f_{L_m}\equiv \frac{\partial f}{\partial L_m},$$

and

$$T_{\mu\nu}=-\frac{2}{\sqrt{-g}}\, \frac{\delta(\sqrt{-g}\,L_m)}{\delta g^{\mu\nu}}.$$

The corresponding trace equation is

$$R\,f_R+3\Box f_R-2\bigl(f-f_{L_m}L_m\bigr)=\tfrac12\,f_{L_m}\,T,$$

where $R$0 [(Harko et al., 2010); (Shukla et al., 2023)].

This structure interpolates between several familiar limits. The linear choice $R$1 reproduces the Einstein equations exactly, and in that case the generalized non-conservation law collapses to the standard $R$2 (Shukla et al., 2023). By contrast, nonlinear matter sectors such as $R$3, mixed couplings such as $R$4, and cosmological forms such as $R$5 introduce effective source terms that cannot be reabsorbed into standard GR without redefining the matter sector (Kavya et al., 2023, Jaybhaye et al., 2024, Sahlu et al., 2024).

A persistent technical point is that the theory depends not only on the functional choice of $R$6, but also on the prescription for $R$7. Across the studies summarized here, both $R$8 and $R$9 are used, with materially different consequences for conservation laws and effective dynamics (Lobato et al., 2022, Sahlu et al., 2024).

2. Energy–momentum balance, extra force, and the role of $L_m$0

A defining property of generic $L_m$1 models is the modified covariant balance equation

$L_m$2

or, equivalently in another common presentation,

$L_m$3

Hence, unless $L_m$4 is constant or a special matter prescription is chosen, matter and geometry exchange energy–momentum [(Harko et al., 2010); (Shukla et al., 2023)].

For a perfect fluid with $L_m$5, the non-conservation law implies non-geodesic motion. The equation of motion may be written as

$L_m$6

with an extra force $L_m$7 orthogonal to the four-velocity, $L_m$8 (Harko et al., 2010). In the weak-field, slow-motion limit, this becomes an Euler-type equation with an additional acceleration

$L_m$9

beyond the Newtonian gravitational acceleration and the ordinary pressure-gradient term (Harko et al., 2010). A related Newtonian-limit treatment gives an effective Poisson relation with $S=\int d^4x\,\sqrt{-g}\,f(R,L_m),$0 under the assumptions used there, together with an extra acceleration $S=\int d^4x\,\sqrt{-g}\,f(R,L_m),$1 (Haghani et al., 2021).

A common misconception is that $S=\int d^4x\,\sqrt{-g}\,f(R,L_m),$2 is unavoidable in every concrete realization. The white-dwarf study based on

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

instead chooses $S=\int d^4x\,\sqrt{-g}\,f(R,L_m),$4, specifically so that $S=\int d^4x\,\sqrt{-g}\,f(R,L_m),$5 remains valid on shell. In that setting, the modified dynamics enter through altered field equations and hydrostatic balance rather than through explicit matter non-conservation (Lobato et al., 2022). This suggests that the choice of matter Lagrangian is not a secondary convention but a central structural ingredient of the theory.

3. FLRW cosmology, autonomous systems, and late-time acceleration

In a spatially flat FLRW background,

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

the modified Friedmann equations can be written as

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

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

with model-dependent matter evolution (Sahlu et al., 2024).

Several late-time cosmological realizations have been analyzed. For

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

a dynamical-systems reformulation yields three critical points: a matter-dominated decelerating point $f(R,L_m)=R/2+L_m$0 with $f(R,L_m)=R/2+L_m$1, $f(R,L_m)=R/2+L_m$2, and $f(R,L_m)=R/2+L_m$3; a de Sitter attractor $f(R,L_m)=R/2+L_m$4 with $f(R,L_m)=R/2+L_m$5 and $f(R,L_m)=R/2+L_m$6; and a scaling solution $f(R,L_m)=R/2+L_m$7 whose stability depends on parameter ranges. In that analysis, $f(R,L_m)=R/2+L_m$8 is stable for

$f(R,L_m)=R/2+L_m$9

and the phase-plane trajectories begin near $8\pi G=c=1$0, pass close to $8\pi G=c=1$1, and end in $8\pi G=c=1$2 (Shukla et al., 2023).

A broader nonlinear matter class,

$8\pi G=c=1$3

has been split into two cosmological cases with $8\pi G=c=1$4 and an uncoupled radiation sector. In Case A,

$8\pi G=c=1$5

the standard radiation $8\pi G=c=1$6 matter $8\pi G=c=1$7 de Sitter sequence occurs only for $8\pi G=c=1$8, with acceleration essentially enforced by the vacuum term. In Case B,

$8\pi G=c=1$9

the interval $$f_R\,R_{\mu\nu}+\bigl(g_{\mu\nu}\Box-\nabla_\mu\nabla_\nu\bigr)f_R -\tfrac12\bigl(f-f_{L_m}L_m\bigr)g_{\mu\nu} =\tfrac12\,f_{L_m}\,T_{\mu\nu},$$0 produces a physical scaling de Sitter future attractor inside the bounded simplex, with $$f_R\,R_{\mu\nu}+\bigl(g_{\mu\nu}\Box-\nabla_\mu\nabla_\nu\bigr)f_R -\tfrac12\bigl(f-f_{L_m}L_m\bigr)g_{\mu\nu} =\tfrac12\,f_{L_m}\,T_{\mu\nu},$$1 and $$f_R\,R_{\mu\nu}+\bigl(g_{\mu\nu}\Box-\nabla_\mu\nabla_\nu\bigr)f_R -\tfrac12\bigl(f-f_{L_m}L_m\bigr)g_{\mu\nu} =\tfrac12\,f_{L_m}\,T_{\mu\nu},$$2 and without introducing $$f_R\,R_{\mu\nu}+\bigl(g_{\mu\nu}\Box-\nabla_\mu\nabla_\nu\bigr)f_R -\tfrac12\bigl(f-f_{L_m}L_m\bigr)g_{\mu\nu} =\tfrac12\,f_{L_m}\,T_{\mu\nu},$$3 (Navarro-Coydán et al., 15 Jan 2026).

Other background constructions emphasize phenomenological reconstruction rather than phase-space completeness. A transit-dark-energy model based on

$$f_R\,R_{\mu\nu}+\bigl(g_{\mu\nu}\Box-\nabla_\mu\nabla_\nu\bigr)f_R -\tfrac12\bigl(f-f_{L_m}L_m\bigr)g_{\mu\nu} =\tfrac12\,f_{L_m}\,T_{\mu\nu},$$4

finds that as $$f_R\,R_{\mu\nu}+\bigl(g_{\mu\nu}\Box-\nabla_\mu\nabla_\nu\bigr)f_R -\tfrac12\bigl(f-f_{L_m}L_m\bigr)g_{\mu\nu} =\tfrac12\,f_{L_m}\,T_{\mu\nu},$$5 (or $$f_R\,R_{\mu\nu}+\bigl(g_{\mu\nu}\Box-\nabla_\mu\nabla_\nu\bigr)f_R -\tfrac12\bigl(f-f_{L_m}L_m\bigr)g_{\mu\nu} =\tfrac12\,f_{L_m}\,T_{\mu\nu},$$6) the cosmographic set approaches

$$f_R\,R_{\mu\nu}+\bigl(g_{\mu\nu}\Box-\nabla_\mu\nabla_\nu\bigr)f_R -\tfrac12\bigl(f-f_{L_m}L_m\bigr)g_{\mu\nu} =\tfrac12\,f_{L_m}\,T_{\mu\nu},$$7

with a transition redshift reported in the range $$f_R\,R_{\mu\nu}+\bigl(g_{\mu\nu}\Box-\nabla_\mu\nabla_\nu\bigr)f_R -\tfrac12\bigl(f-f_{L_m}L_m\bigr)g_{\mu\nu} =\tfrac12\,f_{L_m}\,T_{\mu\nu},$$8 depending on the dataset (Pradhan et al., 2022). A constrained quintessence model with

$$f_R\,R_{\mu\nu}+\bigl(g_{\mu\nu}\Box-\nabla_\mu\nabla_\nu\bigr)f_R -\tfrac12\bigl(f-f_{L_m}L_m\bigr)g_{\mu\nu} =\tfrac12\,f_{L_m}\,T_{\mu\nu},$$9

and a parametrized deceleration function yields late-time quintessence behavior, with joint best-fit values

$$f_R\equiv \frac{\partial f}{\partial R},\qquad f_{L_m}\equiv \frac{\partial f}{\partial L_m},$$0

for the four parameters of that construction (Singh et al., 2022). A separate freezing-quintessence analysis of

$$f_R\equiv \frac{\partial f}{\partial R},\qquad f_{L_m}\equiv \frac{\partial f}{\partial L_m},$$1

with a sinh ansatz for $$f_R\equiv \frac{\partial f}{\partial R},\qquad f_{L_m}\equiv \frac{\partial f}{\partial L_m},$$2 reports that the model approaches a de Sitter epoch identical to $$f_R\equiv \frac{\partial f}{\partial R},\qquad f_{L_m}\equiv \frac{\partial f}{\partial L_m},$$3CDM as $$f_R\equiv \frac{\partial f}{\partial R},\qquad f_{L_m}\equiv \frac{\partial f}{\partial L_m},$$4, and for the quoted best-fit values the adiabatic sound speed remains in $$f_R\equiv \frac{\partial f}{\partial R},\qquad f_{L_m}\equiv \frac{\partial f}{\partial L_m},$$5 for all $$f_R\equiv \frac{\partial f}{\partial R},\qquad f_{L_m}\equiv \frac{\partial f}{\partial L_m},$$6 (Myrzakulov et al., 2024).

4. Perturbations, growth observables, and observational status

Beyond the background, $$f_R\equiv \frac{\partial f}{\partial R},\qquad f_{L_m}\equiv \frac{\partial f}{\partial L_m},$$7 gravity modifies the evolution of matter perturbations through altered effective couplings and matter-dilution laws. In the model

$$f_R\equiv \frac{\partial f}{\partial R},\qquad f_{L_m}\equiv \frac{\partial f}{\partial L_m},$$8

the background was constrained using 57 OHD points, 1048 SNIa distance moduli, and their combination, while perturbations were analyzed with 14 growth-rate data points and 30 $$f_R\equiv \frac{\partial f}{\partial R},\qquad f_{L_m}\equiv \frac{\partial f}{\partial L_m},$$9 measurements. The reported OHD+SNIa best fit is

$$T_{\mu\nu}=-\frac{2}{\sqrt{-g}}\, \frac{\delta(\sqrt{-g}\,L_m)}{\delta g^{\mu\nu}}.$$0

while the $$T_{\mu\nu}=-\frac{2}{\sqrt{-g}}\, \frac{\delta(\sqrt{-g}\,L_m)}{\delta g^{\mu\nu}}.$$1 fit gives

$$T_{\mu\nu}=-\frac{2}{\sqrt{-g}}\, \frac{\delta(\sqrt{-g}\,L_m)}{\delta g^{\mu\nu}}.$$2

Using AIC and BIC, that study finds substantial support from OHD alone but weaker support once SNIa or combined datasets are included (Sahlu et al., 2024).

A more direct growth-rate study considered

$$T_{\mu\nu}=-\frac{2}{\sqrt{-g}}\, \frac{\delta(\sqrt{-g}\,L_m)}{\delta g^{\mu\nu}}.$$3

derived a modified Poisson law with

$$T_{\mu\nu}=-\frac{2}{\sqrt{-g}}\, \frac{\delta(\sqrt{-g}\,L_m)}{\delta g^{\mu\nu}}.$$4

and used 23 $$T_{\mu\nu}=-\frac{2}{\sqrt{-g}}\, \frac{\delta(\sqrt{-g}\,L_m)}{\delta g^{\mu\nu}}.$$5 measurements over $$T_{\mu\nu}=-\frac{2}{\sqrt{-g}}\, \frac{\delta(\sqrt{-g}\,L_m)}{\delta g^{\mu\nu}}.$$6. The quoted best-fit values are

$$T_{\mu\nu}=-\frac{2}{\sqrt{-g}}\, \frac{\delta(\sqrt{-g}\,L_m)}{\delta g^{\mu\nu}}.$$7

with equivalent parameters $$T_{\mu\nu}=-\frac{2}{\sqrt{-g}}\, \frac{\delta(\sqrt{-g}\,L_m)}{\delta g^{\mu\nu}}.$$8, $$T_{\mu\nu}=-\frac{2}{\sqrt{-g}}\, \frac{\delta(\sqrt{-g}\,L_m)}{\delta g^{\mu\nu}}.$$9, and $$R\,f_R+3\Box f_R-2\bigl(f-f_{L_m}L_m\bigr)=\tfrac12\,f_{L_m}\,T,$$0. Relative to $$R\,f_R+3\Box f_R-2\bigl(f-f_{L_m}L_m\bigr)=\tfrac12\,f_{L_m}\,T,$$1CDM, the model produces small $$R\,f_R+3\Box f_R-2\bigl(f-f_{L_m}L_m\bigr)=\tfrac12\,f_{L_m}\,T,$$2 deviations in $$R\,f_R+3\Box f_R-2\bigl(f-f_{L_m}L_m\bigr)=\tfrac12\,f_{L_m}\,T,$$3 at intermediate and high redshifts, a difference of $$R\,f_R+3\Box f_R-2\bigl(f-f_{L_m}L_m\bigr)=\tfrac12\,f_{L_m}\,T,$$4 near $$R\,f_R+3\Box f_R-2\bigl(f-f_{L_m}L_m\bigr)=\tfrac12\,f_{L_m}\,T,$$5, and up to $$R\,f_R+3\Box f_R-2\bigl(f-f_{L_m}L_m\bigr)=\tfrac12\,f_{L_m}\,T,$$6 suppression of growth at $$R\,f_R+3\Box f_R-2\bigl(f-f_{L_m}L_m\bigr)=\tfrac12\,f_{L_m}\,T,$$7 (Goswami et al., 5 Jul 2025).

At the purely background level, multiple fits cluster around $$R\,f_R+3\Box f_R-2\bigl(f-f_{L_m}L_m\bigr)=\tfrac12\,f_{L_m}\,T,$$8CDM-like kinematics but with model-dependent interpretations. The dynamical-systems study reports

$$R\,f_R+3\Box f_R-2\bigl(f-f_{L_m}L_m\bigr)=\tfrac12\,f_{L_m}\,T,$$9

together with $R$00, and interprets the resulting $R$01 as phantom dark energy (Shukla et al., 2023). By contrast, the 2026 autonomous-systems analysis finds Case A consistent with $R$02CDM and Case B within the accelerating window, both at the background level, with $R$03 and $R$04, respectively (Navarro-Coydán et al., 15 Jan 2026). A plausible implication is that background observables alone do not isolate a unique $R$05 phenomenology; perturbative probes are structurally more discriminating.

5. Compact stars, wormholes, and other strong-field realizations

The compact-object literature shows that curvature–matter couplings can be important in regimes where pressure gradients or effective exotic stresses are large. In the white-dwarf application of

$R$06

the static spherically symmetric system was recast into TOV-like equations using the on-shell choice $R$07 and the zero-temperature Hamada–Salpeter equation of state. The resulting equilibrium sequences predict maximum masses above the Chandrasekhar mass limit. The most important reported effect is a significant increase of the mass for stars with radius $R$08 km, while GR is recovered for stars with radii larger than 3000 km independently of $R$09. For a helium core, the sequence goes from $R$10 at $R$11 to $R$12 at $R$13, and the quoted white-dwarf bound is

$R$14

(Lobato et al., 2022).

Wormhole studies have largely focused on nonlinear matter couplings of the form

$R$15

With non-commutative Gaussian and Lorentzian smearing, tideless Morris–Thorne wormholes are obtained that are traversable, horizon-free, and asymptotically flat. In both distributions, the reported energy-condition pattern is qualitatively the same: $R$16, the radial NEC is violated, the tangential NEC holds, the SEC holds, and both radial and tangential DEC are violated. A stability analysis based on the TOV equation and the sound-speed criterion yields

$R$17

as the dynamically stable range for the solutions discussed there (Kavya et al., 2023).

A charged Casimir generalization in the same gravity sector adds a GUP correction and an electromagnetic contribution through

$R$18

That analysis reports that electric charges, GUP effects, and higher model parameter values increase the throat length, and that the NEC remains violated despite a positive contribution from the electromagnetic energy density. The deflection angle was studied via the Gauss–Bonnet theorem, with explicit charge- and GUP-dependent corrections in the large-impact-parameter expansion (Rizwan et al., 2024).

Wormholes sourced by galactic dark-matter profiles have also been constructed in both

$R$19

and

$R$20

Using URC, NFW model-I, and NFW model-II halos, the reported solutions satisfy the standard wormhole requirements with appropriate parameter choices, while violating the null energy constraints in the radial sector. In that sense, dark-matter halos are found to support wormholes within the effective stress structure of these models (Jaybhaye et al., 2024).

6. Early-universe sectors, extensions, and unresolved issues

The framework has also been used in early-universe model building. In gravitational baryogenesis with

$R$21

the effective CP-violating interaction proportional to $R$22 generates a non-zero baryon-to-entropy ratio during radiation dominance. For the numerical choice

$R$23

the paper reports

$R$24

for $R$25 and $R$26, and states that acceptable values lie roughly in

$R$27

(Jaybhaye et al., 2023).

A bouncing cosmology based on

$R$28

finds a non-singular bounce for the ansatz $R$29, with $R$30 at the bounce point and an equation-of-state parameter that crosses the phantom divide line $R$31 when $R$32, i.e. at $R$33. In that treatment, a successful bounce requires NEC violation near $R$34, the SEC is likewise violated around the bounce, and for $R$35 the model can satisfy $R$36 in the immediate vicinity of the bounce, although late-time ghost-like behavior is reported as a possible instability (Koussour et al., 2024).

The direct generalization to $R$37 gravity retains $R$38 as a special case and extends the explicit coupling to the trace $R$39 of the energy–momentum tensor. In that larger framework, the Newtonian limit, generalized Poisson equation, and Dolgov–Kawasaki instability were re-examined, while the $R$40 sector continues to exhibit non-geodesic motion and matter–geometry exchange in generic circumstances (Haghani et al., 2021).

Several open issues recur across the literature. One is the physical status of the matter Lagrangian: $R$41 and $R$42 are both employed, and the choice can determine whether $R$43 vanishes or not in a given model (Lobato et al., 2022, Sahlu et al., 2024). Another is stability: the original formulation explicitly notes the need to reanalyze ghost modes, Dolgov–Kawasaki instability, and singularity formation in the generalized setting (Harko et al., 2010). A further unresolved direction is perturbation-level phenomenology. The 2026 phase-space analysis identifies tensor no-ghost and luminal propagation as minimal consistency conditions, but also argues that growth, lensing, and scalar-sector stability must be studied explicitly through observables such as $R$44 and weak lensing before the viability of nonlinear matter branches can be decisively assessed (Navarro-Coydán et al., 15 Jan 2026).

In aggregate, $R$45 gravity is best understood not as a single model but as a formal umbrella for curvature–matter couplings whose phenomenology depends crucially on the choice of $R$46, the prescription for $R$47, and the physical regime under study. The supplied literature shows that these theories can reproduce GR in controlled limits, generate late-time acceleration, alter structure growth, support super-Chandrasekhar white dwarfs, and sustain wormhole geometries, while simultaneously raising unresolved questions about conservation laws, stability, and precision observational discrimination.