Anisotropic quark stars in ${\boldsymbol {f(R,L_m,T)} }$ gravity

Einstein's General Relativity (GR) theory [1] has stood as the cornerstone of gravitational theory for over a century, successfully describing the dynamics of gravity on both cosmological and astrophysical scales. However, the discovery of the accelerated expansion of the universe and discrepancies in galactic dynamics have motivated the exploration of alternative theories of gravity such as $ f(R,L_m) $ and $ f(R,T) $ [2−5]. The product of the unification of latter theories is $ f(R,L_m,T) $ gravity theory [6, 7], whose novelty is that it modifies the gravitational action to account for observed phenomena beyond the scope of GR.

Quark stars (QSs) [8−10], also known as strange stars or strange QSs, are theoretical entities composed of ultra-dense quark matter, proposed as an alternative to neutron stars (NSs) formed from the gravitational collapse of massive stars [11]. The conjecture is that these hypothetical objects emerge in the aftermath of supernova explosions, where intense pressure and density facilitate the conversion of protons and electrons into quarks [12]. Recent works [13−19] explore several aspects of QSs, including their structure, formation, and observable properties. Investigations into formation mechanisms scrutinize conditions conducive to QS genesis, including the role of quark deconfinement within NS cores during supernova events, elucidating the transition from neutron matter to quark matter under extreme pressure and temperature conditions. Crucially, recent researchers focus on elucidating the equation of state (EOS) of quark matter, delineating the intricate relationship among pressure, energy density, and temperature within these exotic objects [20−23]. Theoretical models and computational simulations are leveraged to unravel the behavior of quark matter under extreme astrophysical environments. In this context, the theory of $ f(R,L_m,T) $ gravity, proposed by Haghani and Harko [6], offers a promising framework for investigating these celestial objects. This theory generalizes and unifies the $ f(R,T) $ and $ f(R,L_m) $ gravity models, where R is the Ricci scalar, T is the trace of the energy-momentum tensor, $ T_{\mu\nu} $, and $ L_m $ is the matter Lagrangian [7]. In $ f(R,L_m,T) $ gravity, the gravitational Lagrangian is given by an arbitrary function of R, T, and $ L_m $, such that $ L_{grav} = f(R,L_m,T) $. The full action in $ f(R,L_m,T) $ gravity theories is expressed as [6]

where g is the determinant of the metric tensor $ g_{\mu\nu} $.

As highlighted previously, $ f(R,L_m,T) $ gravity theories refer to modifications of Einstein's GR theory. These theories are part of the broader framework of modified gravity, which addresses various cosmological and astrophysical phenomena beyond what GR can explain. Notably, $ f(R,L_m,T) $ gravity theories have garnered attention for their potential to elucidate the accelerated expansion of the universe [24, 25], gravitational lensing effects [26−33], and the dynamics of galaxies and galaxy clusters [34−37], all without invoking the existence of dark matter or dark energy. They have been proposed as alternatives to the conventional concept of dark energy. Furthermore, within the realm of modified gravity theory, investigations into dilaton fields [38−43], quasinormal modes [44−46], and inflationary scenarios [47−49] further expand extant understanding of gravitational physics and its implications for the cosmos. These interconnected topics underscore the interdisciplinary nature of modern gravitational research and its quest to uncover the fundamental mechanisms governing the universe's behavior.

This study was conducted to investigate the impact of the $ f(R,L_m,T) $ gravity theory on the internal structure of QSs. We focused on the algebraic function originally proposed in Refs. [6, 7], i.e., $ f(R,L_m,T) = R + \alpha T L_m $, with α being a matter-geometry coupling constant. We examined the matter Lagrangian density, expressed as $ L_m = -\rho $, following which we analyzed its effects on the principal macroscopic properties of compact stars, including mass and radius. We derived the modified Tolman-Oppenheimer-Volkoff (TOV) equations. We conducted numerical analysis on different parameter variations, and present the results in this paper. Further, we examined the static stability criterion, adiabatic index, and sound velocity. Our findings will provide valuable insights into the behavior of compact stars in the context of $ f(R,L_m,T) $ gravity. It should also be noted that although Ref. [50] uses a similar gravitational framework, we used a distinct equation of state (EoS) for quark matter and focused on different aspects of QS properties. The referenced study [50] explores color superconducting quark matter with the CFL phase, whereas our investigation centered around the effects of anisotropy and different EoS models to reveal new insights about the influence of the $ \alpha T L_m $ term on stellar structure.

Although our manuscript correctly implements the first matching condition by ensuring the continuity of the metrics on the surface of the star, one can examine the necessity to seek the additional junction conditions required for a comprehensive matching [51, 52] because it is known that mere metric continuity does not guarantee a smooth matching between interior and exterior metrics. In the context of our chosen theory, $ f(R,L_m,T) $ gravity, it is imperative to consider the conditions imposed by the presence of $ L_m $ in addition to R and T. As outlined in the literature, the matching in $ f(R,L_m,T) $ gravity necessitates the continuity of extrinsic curvature at the surface. The choice of $ f(R,L_m,T) = R + \alpha T L_m $ is motivated by its ability to introduce meaningful deviations from GR while maintaining mathematical simplicity. This form incorporates a coupling between matter and geometry through the trace of the energy-momentum tensor T and the matter Lagrangian $ L_m $, providing an insightful framework to explore high-density environments like QSs. Recent studies have shown that this model leads to observable effects in astrophysical systems, making it a suitable candidate to examine the structure and stability of QSs [6, 7]. Importantly, it reduces to GR when $ \alpha = 0 $, ensuring consistency with known results [2, 3]. After examining this, we discovered that $ f(L_{m}) $ contributes to the thin-shell's stress-energy tensor only if the extrinsic curvature is continuous at the surface. Therefore, our analysis shows that similar to $ f(T) $ [53], the inclusion of $ f(L_m) $ does not introduce further junction conditions. Therefore, our study aligns with the requisite smooth matching conditions established within the framework of $ f(R,L_m,T) $ gravity, ensuring the robustness of our results. However, this matter is purely informational and outside the scope of this article. This article's primary goal is to examine how $ f(R,L_m,T) $ gravity affects compact stars' interior structure. It is anticipated that this hypothesis will be most evident in these stars' high-density cores.

The paper is structured as follows: In Section II, we review $ f(R, L_m, T) $ gravity theory and serve the gravitational field equations of the theory. Section III discusses the gravitational field equations in $ f(R,L_m,T) $ theories, while in Sec. IV, we prescribe a quark matter EoS and a specific quasi-local EoS that describe local anisotropy for QSs. Numerical results and discussions for different parameter variations are presented in Sec. V, followed by an examination of the static stability criterion, adiabatic index, and sound velocity in Sec. VI. Finally, we draw our conclusions in Sec. VII. (We adopt geometrized units throughout this exposition, setting $ c = G = 1 $, while retaining physical units for clarity in presentation.)

In this section, we review the modifications to classical gravity theories in the context of generalized $ f(R, L_m, T) $ gravity [6, 7], where arbitrary functions of the Ricci scalar R, the matter Lagrangian $ L_m $, and the trace of the stress-energy tensor T are included. The gravitational Lagrangian in this modified framework is $ f(R, L_m, T) $, allowing for a more comprehensive description of gravitational interactions in the presence of matter.

The gravitational Lagrangian in this unified framework encompasses arbitrary functions of the Ricci scalar R, the trace T of the energy-momentum tensor $ T_{\mu\nu} $, and the matter Lagrangian $ L_m $, yielding $ L_{grav} = f(R,L_m,T) $. Consequently, the complete action in $ f(R,L_m,T) $ gravity theory results in Eq. (1). The variation of action (1) with respect to the inverse metric $ g^{\mu\nu} $ yields the following field equations in $ f(R,L_m,T) $ gravity:

where $ \Box \equiv \partial_{\mu}(\sqrt{-g}g^{\mu\nu}\partial_{\nu})/\sqrt{-g} $, and $ f_{R} $, $ f_{T} $, $ f_{L} $ denote the partial derivatives of f with respect to R, T, and $ L_{m} $, respectively. We adopted $ L_m = -\rho $ primarily for its physical relevance in modeling compact stars. The energy density ρ directly contributes to the gravitational source term, aligning with conventional understanding in GR, where matter density is central to the formation of gravitational fields. In contrast, using $ L_m = p $ would introduce significant differences in the mass-radius relations and stability, as shown in [7].On the other hand, alternative forms of $ L_m $ could be explored in future work. Although $ L_m = p $ has been explored in some studies, as in Ref. [2], we chose $ L_m = -\rho $ to maintain consistency with the standard treatment of compact objects and simplify numerical comparisons with GR. This choice ensures that the matter's gravitational effects remain physically intuitive, with ρ acting as the primary source. $ R_{\mu\nu} $ stands for the Ricci tensor, $ \nabla_{\mu} $ represents the covariant derivative concerning the symmetric connection associated with $ g_{\mu\nu} $, and $ \tau_{\mu\nu} $ is a new tensor defined as [6]

The form of $ f(R,L_m,T) $ dictates the nature of the gravitational dynamics. Specifically, when $ f(R,L_m,T) = f(R) $, Eq. (2) reduces to the field equations of metric $ f(R) $ gravity [54−56]. Similarly, for $ f(R,L_m,T) = f(R,T) $, the $ f(R,T) $ gravity model is recovered, while $ f(R,L_m,T) = f(R,L_m) $ yields the field equations of the $ f(R,L_m) $ theory. Furthermore, $ f(R,L_m,T) = R $ retrieves the standard field equations of GR: $ R_{\mu\nu}-\dfrac{1}{2}g_{\mu\nu}R = 8\pi T_{\mu\nu} $.

The covariant divergence of the field equations (2) leads to the non-conservation equation of the energy-momentum tensor $ T_{\mu\nu} $:

where $ f_m = f_T + \dfrac{1}{2}f_L $ and $ A_\nu = \nabla^\mu(f_T \tau_{\mu\nu}) $.

The trace of the field equations yields a second-order differential equation:

where τ represents the trace of the tensor $ \tau_{\mu\nu} $. For the specific functional form $ f(R,L_m,T) = f(R) $, this equation reduces to the well-known dynamical equation for the Ricci scalar in pure $ f(R) $ gravity theories [57, 58].

For simplicity, we focus on the algebraic function originally proposed in Ref. [6], i.e., $ f(R,L_m,T) = R + \alpha T L_m $, with α denoting a matter-geometry coupling constant. In this context, Eqs. (2) and (5) simplify to

and

respectively, where $ G_{\mu\nu} $ denotes the Einstein tensor. Remarkably, the Einstein field equations $ G_{\mu\nu} = 8\pi T_{\mu\nu} $ and the conservation equation $ \nabla^{\mu}T_{\mu\nu} = 0 $ are retrieved when $ \alpha = 0 $.

In this section, we will explore anisotropic QSs within the framework of the $ f(R,L_m,T) = R + \alpha TL_m $ gravity model. To begin our analysis, we will concentrate on the static and spherically symmetric spacetime. To do so, we initially adopt the metric ansatz described by the following line element:

where $ \nu(r) $ and $ \lambda(r) $, both functions of the radial coordinate r, are the two unknown functions. Within a spherically symmetric star, the matter source is modeled as an anisotropic fluid, indicating that the radial pressure $ p_r $ differs from the transverse pressure $ p_{\perp} $. The energy-momentum tensor for the static configuration is given by:

where ρ represents the energy density and $ \Delta \equiv p_{\perp} - p_r $. Here, $ u_{\mu} $ denotes the fluid 4-velocity, satisfying $ u_{\mu}u^{\mu} = -1 $, which can be expressed as $ u^\mu ={\rm e}^{-\nu/2}\delta_0^\mu $. Consequently, $ T_\mu^\nu = \text{diag} (-\rho, p_r, p_{\perp}, p_{\perp}) $ and $ T = -\rho + p_r + 2 p_{\perp} $.

Given that the Lagrangian density $ L_m $ corresponding to the matter source is not unique, we have the option to choose either $ L_m = p $ or $ L_m = -\rho $ (see [59] for further discussion). Although consensus on which Lagrangian to consider is lacking, we opted to assume $ L_m = -\rho $ and derive the equations of hydrostatic equilibrium for QSs within the context of $ f(R,L_m,T) $ gravity. This assumption allows us to rewrite Eq. (6) as:

Utilizing the spherically symmetric metric (8) and Eq. (10), the non-zero components of the field equations become

where the prime denotes a derivative with respect to r. On the other hand, the covariant divergence of Eq. (7) yields the modified TOV equations:

where

According to Eq. (14), the standard conservation equation does not hold for this theory. It can be readily shown that when $ \alpha = 0 $, the standard conservation equation is recovered.

Figures 1 and 2 depict various physical quantities as functions of the radial coordinate r, illustrating the sensitivity of the system to parameters α and β, respectively.

Figure 1.  (color online) From top to bottom, energy density ρ, radial $ p_r $ and transverse $ p_t $ pressures as functions of the radial coordinate r. We start by varying $ \alpha \in [-0.2,0.2] \mu_1 $ where $ \mu_1 = 10^{-79} $ s⁴/kg². We define the other parameters, $ B = 60 $ MeV/fm³ and $ \beta = 1.5 $. The dashed black line represents the GR solution $ (\alpha = 0) $ for an anisotropic fluid distribution.

Figure 2.  (color online) From top to bottom, energy density ρ, radial $ p_r $ and transverse $ p_t $ pressures as functions of the radial coordinate r. We start by varying the anisotropy parameter $ \beta \in [-1.5, 1.5] $, whereas the other parameters are $ B = 60 $ MeV/fm³ and $ \alpha = -0.2 \mu_1 $ with $ \mu_1 = 10^{-79} $ s⁴/kg², respectively. The black dashed line denotes an isotropic arrangement in $ f(R,L_m,T) $ gravity.

Before exploring the internal composition of compact stars in $ f(R,L_m,T) = R + \alpha T L_m $ gravity, we introduced a mass function $ m(r) $ through the relation:

where $ m(r) $ represents the total mass enclosed within the radius r. Subsequently, by substituting this expression into Equations (11)−(14), we ultimately derive

In this framework, we encountered five unknown functions with three differential equations. However, introducing two appropriate ad hoc assumptions enabled us to close the system of equations. In this context, we considered an EoS relating radial pressure to the energy density of the fluid and a quasi-local EoS proposed by Horvat et al. in Ref. [60], which offers reliable solutions and aligns effectively with the GW170817 data. Based on these considerations, we numerically solved the structure equations with the appropriate boundary conditions, as follows:

ensuring regularity at the center of the star, where $ \rho_c $ denotes the central energy density. Finally, we solved these equations by integrating outward until the pressure vanished, i.e., $ p_r(R) = p_{\perp}(R) = 0 $, where R is identified as the radius of the star. Then, the interior solution was matched with the exterior Schwarzschild vacuum solution:

where M signifies the total mass of the star.

Observations of compact objects via gravitational wave (GW) signals or other recent astrophysical data have provided strong and direct evidence in favor of more massive compact stars. Interestingly, these observations also impose constraints on the NS EoS, suggesting that their internal composition likely includes free quarks in their core, rather than only neutrons and other non-fundamental particles. The central density is expected to reach several times higher than the nuclear saturation density ($ n_0 \approx 0.148\ \text{fm}^{-3} $), leading to the possible formation of quark matter [61, 62]. This density is typical within atomic nuclei and provides a reference point for understanding the extreme conditions that may lead to the deconfinement of quarks. This highly compressed nuclear matter is anticipated to undergo a phase transition to deconfined quark matter, releasing its constituent quarks and gluons. Therefore, the EoS of quark matter plays a pivotal role in shaping the structure of the star beyond nuclear saturation density. From Sec. IV onward, we use "NS" units ($ \hbar = c = 1 $) consistently in discussing the numerical results and the value of α. To ensure that all physical quantities are presented coherently throughout the manuscript, the value of α mediated by $ \mu_1 $ is expressed using units appropriate for NS contexts.

In this work, we focused on the applications of strange matter, which are described by the MIT bag model EoS for modeling compact stellar objects. The MIT bag model is a simple phenomenological model for quark matter proposed in the 1970s to explain hadrons in terms of quarks [63]. In the model, it is assumed that quarks are asymptotically free and confined to a spherical region of space by the bag constant B. Consequently, the MIT bag model establishes a relation between energy density and pressure, given by

This expression indicates that at $ \rho = 4B $, the external pressure acting on a bag vanishes. The constant B, referred to as the bag pressure, is typically constrained within a range of $ 57 \leq B \leq 92 $ MeV/fm³. This range is supported by various studies in the field [64, 65]. Following [66], we used $ B = 60\ {\rm{MeV}}/{\rm{fm}}^3 $ to align the model with observational astrophysical data. Noteworthily, this chosen value also satisfies the $ 2M_\odot $ constraints within the framework of GR. For local anisotropy, we employed the quasilocal EoS suggested by Horvat et al. [60]. The choice of this model was driven by two key factors: (a) The simplicity and physical consistency of the Horvat model, which connects the anisotropy directly to the local compactness parameter $ \mu = \dfrac{2m(r)}{r} $, makes it an intuitive choice for compact stars. (b) Furthermore, as shown in [7], the Horvat model emerges naturally from the modified TOV equation for ${\rm d}p_r/{\rm d}r$, making it a suitable and physically grounded approach for studying anisotropic compact stars. Although other anisotropic models could have been considered, such as those based on pressure or density gradients, the Horvat model offers a straightforward and analytically tractable way to explore anisotropic effects. This model's ability to smoothly transition between isotropic and anisotropic cases by varying the parameter β further supports its use in our analysis. Within this framework, various solutions have been obtained in both GR and modified gravity theory. Thus, we express the quasi-local equation as

where the free parameter β plays a crucial role in measuring the deviation from isotropy, taking positive or negative values, see Refs. [67−73] for detailed discussions. In our calculations, we assumed $ \beta \in [-1.5,1.5] $. The quantity $ \mu \equiv 2m(r)/r $ represents the local measure of compactness. Furthermore, the choice of Eq. (23) ensures that $ \Delta = 0 $ at the center, thereby recovering the isotropic case. In Figs. 1 and 2, we depict the energy density ρ, radial pressure $ p_r $, and transverse pressure $ p_t $ as functions of the radial coordinate r for two different parameter sets (see Tables 1 and 2 for more details). It is evident from these figures that the energy density of the QS is non-zero at the surface, and the pressure components vanish simultaneously at the surface of the star, i.e., $ p_r \left( r \rightarrow R \right) = p_{\perp} \left( r \rightarrow R \right) = 0 $. Clearly, when $ \beta = 0 $, the anisotropy factor vanishes at the center, ensuring regularity in the interior.

$\alpha \times 10^{-78} \dfrac{s^4}{kg^2}$	M/$ M_{\odot} $	R/km	$ \rho_c $/(MeV/fm3)	$ M/R $
−0.2	2.72	11.68	847	0.345
−0.1	2.63	11.57	847	0.337
0.0	2.49	11.42	847	0.323
0.1	2.29	11.34	791	0.299
0.2	1.96	11.15	735	0.260

Table 1.  Summary of the resulting properties of anisotropic QSs for the variation of α in $ f(R,L_m,T) $ gravity, as briefly discussed in Sec. V.A.

β	M/$ M_{\odot} $	R/km	$ \rho_c $/(MeV/fm3)	$ M/R $
−1.5	1.91	10.39	1,415	0.272
−1.0	2.03	10.59	1,334	0.284
−0.5	2.16	10.78	1,253	0.298
0.0	2.30	10.96	1,173	0.311
0.5	2.44	11.20	1,052	0.323
1.0	2.58	11.47	931	0.334
1.5	2.72	11.67	851	0.346

Table 2.  Summary of the resulting properties of anisotropic QSs for the variation of β in $ f(R,L_m,T) $ gravity, as briefly discussed in Sec. V.B.

To investigate the influence of pressure anisotropy on the stellar structure, we first numerically integrated Equations (18), (19), (22) and (23) with the boundary conditions (20). By varying the central density, we obtained mass as a function of the radial coordinate r, the mass-radius ($ M-R $) and the compactness $ (M-M/R) $ relations. Here, we used several values of $ \alpha \in [-0.2,0.2] \mu_1 $ where $ \mu_1 = 10^{-79} $ s⁴/kg² [7] and $ \beta \in [-1.5, 1.5] $ to demonstrate the results for anisotropic stars in $ f(R,L_m,T) = R+ \alpha T L_m $ gravity.

From top to bottom, we plotted $ M-r $, $ M-R $ and $ (M/R-M) $ curves for the variation of $ \alpha \in [-0.2,0.2] \mu_1 $ where $ \mu_1 = 10^{-79} $ s⁴/kg² in Fig. 3. The other parameters were $ B = 60 $ MeV/fm³ and $ \beta = 1.5 $, respectively. For this particular model and dataset, we observed that the mass of the QS was inversely proportional to the values of α. This observation is consistent with the results reported in Ref. [7]. Indeed, the reported value of maximum mass for QS is 2.72 $ M_{\odot} $ with radius $ R = 11.68 $ km for $ \alpha = -0.2 \mu_1 $. For completeness, we list the maximum masses and corresponding radii of QSs obtained from the numerical calculation in Table 1 based on which we can quantify how coupling constant α contributes to $ M_{\rm{max}} $ estimates. Moreover, the maximum mass in GR (for anisotropic fluid) is 2.49 $ M_{\odot} $ with radius $ R = 11.42 $ km. In other words, we can say that sufficiently negative values of α increase the maximum mass and corresponding radius, see Table 1. It is worth noting that the deviation between GR and $ f(R,L_m,T) $ gravity is significant. Putting the observational constraints into consideration, we considered several observations, in particular, PSR J0952-0607 with mass $ M = 2.35 \pm 0.17 M_{\odot} $ (Yellow) [74], PSR J0740+6620 with the pulsar mass $ M = 2.08 \pm 0.07 $ $ M_{\odot} $ (Orange) [75] and PSR J0348+0432 with the mass of $ M = 2.01 \pm 0.04 M_{\odot} $ (Pink) [76]. We can also discuss our results in light of constraints obtained from the detection of GW190814 event with a mass of $ 2.59^{+0.08}_{-0.09} M_{\odot} $ (Green) [77] and GW170817 event [78]; see middle panel of Fig. 3. More importantly, our results can be described as consistent with the main constraint from recent observational astrophysical data. Finally, we plotted $ (M-M/R) $ diagram in the lower panel of Fig. 3. We observed that for different values of α the maximum compactness differs according to the maximum masses. It can also be noted from Table 1 that the value of maximum compactness increases with decreasing values of α, which is consistent with the $ M-R $ diagram. Moreover, the value of maximum compactness lies within the range of $ 0.260 < M/R < 0.345 $. It is worth noting that theoretical considerations [79, 80] suggest α should be small to avoid large deviations from GR predictions, especially for stable, high-density QS configurations. Observationally, fitting to specific QS data could provide indirect constraints. Large negative α enhances mass, while positive α tends to reduce it, as shown in Table 1. These values influence compactness and stability within observed astrophysical ranges.

Figure 3.  (color online) Mass as a function of radial coordinate (top), the mass-radius (middle), the mass-compactness (bottom) relations for MIT bag model EoS used in this study with parameters mentioned in Table 1. The horizontal bands correspond to the observed masses of pulsars: PSR J0952-0607 [74], PSR J0740+6620 [75], and PSR J0348+0432 [76] with different colour bars. Furthermore, we include the constraint from GW190814 [77] and GW170817 event [78]. In the relevant graphs, the black dashed line represents the GR solution $ (\alpha = 0) $ for an anisotropic fluid distribution.

We further analyzed the effect of anisotropic parameter $ \beta \in [-1.5, 1.5] $ on the stellar structure in Fig. 4. The sequence of curves that describe the behavior of $ M-r $, $ M-R $ and $ (M/R-M) $ for the parameter set: $ B = 60 $ MeV/fm³ and $\alpha = -0.2 \times 10^{-78}~\rm s^4/kg^2$. As can be observed from the figure, the mass of the QS increases for increasing values β. The result of the main aspects are summarized in Table 2, and Fig. 4 (middle panel) shows some associated $ M-R $ relations for QSs. Depending on the model parameter, the value of maximum mass is $ M_{\rm{max}} = $ 2.72 $ M_{\odot} $ and the corresponding radius is $ 11.67 $ km at $ \beta = 1.5 $. As for the values of $ \beta = 0 $ (the isotropic fluid) the maximum mass goes upto $ M_{\rm{max}} = $ 2.30 $ M_{\odot} $ and the corresponding radius is $ 10.96 $ km. As expected, the presence of pressure anisotropy leads to a considerable changes in the $ M-R $ relations. Concerning the constraints, we see that our model is compatible with maximum QS masses. The $ (M-M/R) $ diagram is plotted in Fig. 4 (lower panel). The maximum compactness of QSs monotonically increases with increasing β values. Finally, the results are presented in Table 2, where, as can be observed, the maximum compactness corresponding to the maximum mass lies within the range of $0.272 < M/R < 0.346$.

Figure 4.  (color online) The mass as a function of radial coordinate (top), the mass-radius (middle), the mass-compactness (bottom) relations for MIT bag model EoS used in this study with parameters mentioned in Table 2. Imposing constraints from astrophysical observations are the same as of Fig. 2. In the relevant graphs, the black dashed line represents an isotropic configuration in $ f(R,L_m,T) $ gravity.

Apart from the mass-radius relations, the most important issue is related to stability of the configuration. A detailed study is presented subsequently in this section.

Let us now focus on the stability of the equilibrium configuration through static stability criterion (SSC) [81, 82]. This criteria has been presented in $ M-\rho_c $ plane, where M is the mass and $ \rho_c $ is the central density of the star. In a similar vein, SSC has also been applied in modified gravity theory , see Refs. [83−85] for more. We should also note that this is a necessary but insufficient condition. The mathematical form of these inequalities is as follows

In Fig. 5, we show the $ M-\rho_c $ curves for the consider cases mentioned above. For plotting, we use the same parameters as shown in Fig. 1 and 2, respectively. Here, we see that the total mass is an increasing function of the central density, reaching a point where $ (M_{\text{max}}, R_{M_{\text{max}}}) $ exists. This point is known as a boundary point (indicated in the figure by pink points) which can separate the stable configuration region indicated by ${\rm d}M/{\rm d}\rho_c > 0$ from the unstable one.

Figure 5.  (color online) This profile is for $ M-\rho_c $ curves. The same parameter set as used for Figs. 1 to 2 is considered. The pink points are representing as a boundary wall that separates the stable configuration region indicated by ${\rm d}M/{\rm d}\rho_c > 0$ from the unstable one.

An additional test was performed to evaluate the stability of the configuration via adiabatic index, γ. The stability of compact objects was extensively studied, with reference to notable contributions by Chandrasekhar [86, 87], who laid the groundwork for understanding dynamical instability in compact stars. His work provides key insights into the limits of stability for stellar objects, including NSs and white dwarfs. Maintaining the spherical symmetry of the background, Chandrasekhar performed the dynamical stability using the theory of infinitesimal radial perturbations [86]. The adiabatic index (γ) is given by

The expression (26) is associated with the sound speed, and the subscript S indicates the derivation at constant entropy. γ has some restrictions related with the dynamical instability of the spherical static object. For polytropic stellar model, authors in [88] have shown that the adiabatic index γ is greater than $ \gamma> \gamma_{cr} = 4/3 $, below which configurations are unstable against radial perturbations. For more details about the role of adiabatic index on stellar structure, see Ref. [89]. In Fig. 6, we display the dependence of γ for two considered cases as a function of radial coordinate r. Based on the figures, it can be concluded that $ \gamma > \gamma_{cr} $ holds for both cases, ensuring the stability of the configuration under consideration.

Figure 6.  (color online) The adiabatic index γ is plotted as a function of radial coordinate r. The same parameter set used for Figs. 1 and 2 was considered. The value of $ \gamma_{cr} = 4/3 $ is indicated by the brown dot-dash line.

The propagation of sound speed is another important criteria for checking the stability of QSs. The sound speed within the star is given by $v^2_s = {\rm d}p_{r, \perp}/ {\rm d}\rho$, and lies within the range of $ 0<v^2_s <1 $because the speed of sound does not exceed the speed of light. Furthermore, the sound speed along radial direction is constant inside the star, whereas the tangential velocity has been plotted (both cases) in Fig. 7. As can be observed from those figures, the sound speed along transverse direction lies within the specified range. Thus, we can conjecture that QSs with anisotropic pressure can exist in $ f(R,L_m,T) $ gravity.

Figure 7.  (color online) Sound speed along transverse direction within the stellar interior. The same parameter set as used for Figs. 1 and 2 is considered.

In this study, we investigated the behavior of QSs within the framework of $ f(R,L_m,T) $ gravity, focusing on the specific model $ f(R,L_m,T) = R + \alpha T L_m $. This model offers a novel perspective on the internal structure of compact stars, allowing us to explore the impact of modifications to the gravitational action on the properties of QSs.

Our analysis revealed several key findings. First, we derived the modified TOV equations for QSs in $ f(R,L_m,T) $ gravity, considering a specific choice of the matter Lagrangian density $ L_m = -\rho $. These equations govern the hydrostatic equilibrium of QSs and provide valuable insights into their macroscopic properties.

Next, the modified TOV equations are solved numerically with proper boundary conditions to demonstrate the results for anisotropic stars in $ f(R,L_m,T) = R+ \alpha T L_m $ gravity. Our results showed the sensitivity of QS properties to the MIT bag parameter, highlighting the importance of considering various scenarios when studying compact stars in modified gravity theories. The model was implemented such as to vary the coupling parameter α and anisotropic parameter β to determine the maximum mass of QSs. The results are interesting because of the inverse relationship between the maximum mass for sufficiently negative values of α, whereas the situation is reverse for increasing positive values of β. Furthermore, we analyzed the static stability criterion, adiabatic index, and sound velocity of QSs to assess their dynamical stability. Our calculations showed that QS models within the considered parameter ranges satisfy stability criteria, indicating the robustness of these configurations against radial oscillations and sound propagation.

Our modified QS model showed consistency with observed high-mass compact stars, such as PSR J0740+ 6620, which had an estimated mass of $ M = 2.08 \pm 0.07 \, M_{\odot} $ and radius constraints provided by NICER data [90]. The influence of the coupling parameter α in the $ f(R, L_m, T) $ gravity model supported the attainment of higher maximum masses in QS configurations, aligning with this observational range. Additionally, GW data, particularly from events GW190814 and GW170817, imposed stringent constraints on compact star models. The high secondary mass in GW190814 ($ 2.50-2.67 \, M_{\odot} $) presented a challenge for traditional neutron star EoS but aligned with the higher mass predictions of our modified model. By adjusting the coupling parameter, we also achieved QS configurations compatible with the tidal deformability limits from GW170817 [91], indicating that our model can satisfy these observational constraints. This agreement with astrophysical data underscored the potential of $ f(R, L_m, T) $ gravity as a framework for understanding the properties of ultra-dense compact objects.

Overall, our study contributes to the growing body of research on compact stars in alternative gravity theories. By exploring the behavior of QSs in $ f(R,L_m,T) $ gravity, we deepen our understanding of the gravitational effects on extreme astrophysical objects. In conclusion, the nonlinear term $ \alpha T L_m $ might influence additional QS characteristics, including compactness and rotation, through adjustments to the stress-energy tensor and equilibrium equations. This term might suggest surface effects, such as tension, and impact the moment of inertia during rapid rotation, which we intend to explore further in a separate study.

Future investigations could extend to other modified gravity models and explore additional observational consequences for QSs such as GW signatures and electromagnetic emissions. Moreover, the implications of our findings for the formation and evolution of quark stars in diverse astrophysical environments are worth exploring. Additionally, incorporating magnetic fields or studying rotating quark stars could further illumine the complex dynamics of these exotic objects.

The authors would like to thank the editor and reviewer for their insightful comments and constructive suggestions, which have greatly enhanced the clarity and quality of this manuscript.