Effect of pressure anisotropy on quark star structure in the Starobinsky model

The study of compact objects, such as QSs, has driven extensive research efforts owing to their unique role as astrophysical laboratories for high-density matter and potential indicators of physics beyond general relativity (GR). Numerous modified theories of gravity have been proposed to address fundamental cosmological questions, including dark energy and inflation. Among these theories, the Starobinsky model, which augments the Einstein-Hilbert action with a term proportional to the Ricci scalar squared, $ f(R) = R + aR^2 $, offers a simple yet compelling framework. This model addresses early-universe inflation and is compatible with cosmic microwave background data, making it a promising candidate for investigating gravitational interactions at astrophysical scales [1, 2].

Astrophysical observations of supermassive neutron stars [3, 4] and recent gravitational wave detections, such as GW190814 [5], suggest the need for theories that extend GR in strong-field regimes [6]. These detections challenge traditional models, particularly in describing stellar structures under high densities and extreme conditions. In this context, the $f(R) $ models [7−12] and their particular realization in Starobinsky gravity offer insight into stellar stability and compactness properties that may deviate from GR predictions, which makes them attractive for modeling compact stars like QSs [13−21]. Gravitational wave signals from GW 190814 have suggested the presence of a secondary object that is a neutron star with a mass of 2.6 $ M_{\odot} $, and modified theories of gravity have been proposed to describe this phenomenon [20, 21]

Another essential aspect of QS structure is pressure anisotropy, where the radial pressure, $ P_r $, differs from the tangential pressure, $ P_\perp $. Pressure anisotropy in compact stars is a natural consequence of various physical phenomena, including strong magnetic fields, rotation, and particle interactions at ultra-high densities [22, 23]. In quark matter, pressure anisotropy affects the star's mass-radius relation, stability, and compactness, making it crucial to include such anisotropic features when studying QSs under modified gravity [14, 24]. Indeed, pressure anisotropy has been shown to significantly alter the maximum mass and radius of compact stars, highlighting its relevance in the structure of QSs modeled within the Starobinsky framework [25].

This study aims to investigate the effect of pressure anisotropy on the structural properties of QSs within the Starobinsky model. We employ a quark matter equation of state (EoS) tailored for anisotropic conditions, combined with modified TOV [26] equations in Starobinsky gravity, to analyze the impact on mass-radius relationships and compactness. Our results extend existing studies on compact stars in $f(R) $ theories by specifically addressing how anisotropy modifies the predictions of Starobinsky gravity regarding QS structure.

This paper is structured as follows: In Sec. II, we formulate the problem and derive the field equations for Starobinsky's model under spherical symmetry. Section III presents the anisotropic quark matter EoS used in this study. In Sec. IV, we discuss our numerical results, focusing on the mass-radius relationship and stability conditions under varying anisotropic parameters. Section VII provides concluding remarks and future perspectives. Throughout the paper, we adopt the signature $ (-, +, +, +) $ and set $ c = 1 $.

We begin with the most-studied modified gravity theory, so-called $ f(R) $ gravity, in which the Lagrangian is replaced by a general function of the Ricci scalar, and the action is given by

where $ f(R) $ is a function of the Ricci scalar only, and g denotes the determinant of the metric $ g_{\mu\nu} $. $ S_{\rm matter} $ denotes the action of the matter field depending on the metric tensor and matter field $ \psi_i $. To prevent pathological scenarios such as tachyonic instabilities and ghosts, viable $ f(R) $ theories must meet the requirements [7, 8]

respectively. For our purposes, we shall consider a specific form of Starobinsky's model, that is, $ f(R)=R+aR^2 $, where $ a \geq 0 $, which is in agreement with the inequalities given in Eq. (2) [27]. Here, the free parameter a has dimensions of $ [mass]^{-2} $, and we may write it in the form $ a = 1/M^2 $, stating that the mass scale M is now the theory's free parameter. Selecting $ a = 0 $allows us to obtain the GR solution.

In general, solving fourth order differential equations in $ 4D $ spacetime is problematic. Hence, we adopt a conformal transformation by introducing the scalar field φ and metric $ \tilde{g}_{\mu \nu} $ [28−30]:

where $ A (\varphi = \exp(-\phi/\sqrt{3}) $. Using this transformation, one can rewrite the action in Eq. (1) in the Einstein frame as

where the scalar-field potential $ V(\varphi) $ takes the form [28, 29]

Taking variation with respect to the metric $ \tilde{g}_{\mu \nu} $ and scalar field φ, one arrives at the following modified field equation in the Einstein frame

where $ T^{\varphi}_{\mu \nu} $ is the the energy-momentum tensor corresponding to the scalar field and $V_{,\varphi} \equiv \dfrac{{\rm d}V(\phi)}{{\rm d}\phi}$. Moreover, the coupling constant $ \alpha (\varphi) $ is related with [30, 31]

Through the formula $ \tilde{T}_{\mu \nu} = A(\phi)^2 T_{\mu \nu} $, the energy-momentum tensor of the Einstein frame $ \tilde{T}_{\mu \nu} $ is connected to that of the Jordan frame $ T_{\mu \nu} $. We consider the matter source to be an anisotropic fluid described by

where ρ is the energy density, and $ P_{\perp} $ and $ P_r $ are the pressures perpendicular and parallel to the spacelike vector $ \chi_\nu $, respectively. Then, the only nonzero diagonal components in the energy-momentum tensor are as follows: $ T^\nu_i = \left( -\rho, P_r, P_{\perp}, P_{\perp} \right) $.

For an anisotropic fluid, the energy density, pressure components, and 4-velocity in the two frames are interconnected through the equations [30, 31]

where the tilde indicates the Einstein frame.

Now we look at a static spherically symmetric metric that describes non-rotating solutions of these field equations inside a compact star of $ f(R) $ gravity,

where ${\rm d} \Omega^2 = {\rm d}\theta^2 + \sin^2\theta {\rm d}\vartheta^2$ is the line element on the unit 2-sphere. Moreover, all the metric functions depend only on the radial coordinate r.

Under these assumptions, we reach the following set non-zero components of the field equations (for a more detailed derivation of the following field equations, see Appendix A),

At the stellar surface ($ r=R_S $), where $ P_r=0 $, the anisotropy term

vanishes because our quasi-local ansatz, $ \Delta = P_{\perp}-P_r = \beta P_r \zeta $, ensures that $ P_{\perp}=0 $when $ P_r=0 $. This guarantees that the standard TOV boundary condition is preserved at the surface, while anisotropy continues to influence the star's interior structure. For a comprehensive description of the examined configuration, we must add an EoS for the fluid to the aforementioned equations. Note that the absence of a scalar field allows one to recover the standard conservation equation.

To close the modified TOV equations (14)−(17) numerically, one needs boundary conditions for solving the interior and exterior problems simultaneously. Thus, we require boundary conditions in the natural Einstein frame to maintain regularity at the origin:

where $ \rho_c $ and $ \phi_c $ are the central values of energy density and scalar field, respectively, while at infinity

The radius of the star can be identified using the condition

where the radial pressure at the surface of the star vanishes. Referring to [30], the condition $\dfrac{{\rm d}\varphi}{{\rm d}r}(0) = 0$ guarantees the regularity of the scalar field ϕ, which consequently ensures the regularity of Φ at the stellar interior $ r = 0 $. Moreover, $ \Lambda(0)=0 $ is necessary to ensure that the Einstein frame geometry at the center is regular. The Einstein and Jordan frame metrics are interconnected through a nonsingular conformal factor, which guarantees the regularity of the Jordan frame geometry at the star's center. In contrast, we choose boundary conditions that can fulfill the asymptotic flatness criterion at infinity, which requires $ \lim_{r\to\infty} V(\varphi(r)) = 0 $ and gives $ \lim_{r\to\infty}\varphi(r) = 0 $. The criteria expressed in Eq. (20) guarantee asymptotic flatness in both the Einstein and Jordan frames.

The coordinate radius of the star is determined through Eq. (21), while the physical radius of the star, as measured in the physical Jordan frame, is expressed by

Next, we show the structure equation for a QCD-driven EoS and an EoS with an anisotropic profile that acts as a matter source for QSs in the Starobinsky model.

The latest astronomical findings of massive pulsar masses support the possibility of a quark matter or QCD matter (quantum chromodynamic) core. Therefore, both the theoretical and experimental communities have focused their efforts on QCD matter in extreme environments, such as high temperatures and/or density at the core of compact objects, for decades. Interestingly, the perturbative quantum chromodynamics (pQCD) formalism becomes more reliable as the quark mass increases [32]. This means that this method could be useful for studying stars that are densely packed.

The current study focuses on investigating the properties of compact stars made up of interacting quark matter (IQM) EoS, taking into account the superconducting effect [33] and quantum chromodynamics (pQCD) corrections. Based on the correction terms, we rewrite the thermodynamic potential Ω as [34, 35]

where μ and $ \mu_e $ are the respective average quark and electron chemical potentials. The unpaired free quark gas contribution is represented by the first term. The pQCD contribution from one-gluon exchange for gluon interaction to the $ O(\alpha_s^2) $ order is expressed by the second term with $ (1-a_4) $. The quartic coefficient $ a_4 $ is a constant accounting for the strong interactions between quarks [34]. In the above expression, the third term with accounts for the correction from the finite mass of the strange quark, if relevant. Moreover, the term with the gap parameter Δ represents the contribution from color superconductivity. Different states of color-superconducting phases are represented by $ (\xi_4, \xi_{2a}, \xi_{2b}) $. Additionally, the effective bag constant, $ B_{\rm eff} $, accounts for the nonperturbative contribution from the QCD vacuum. Generally, the range of $ B_{\rm eff} $ lies between $ 57 \geq B_{\rm eff} \geq 92 $ MeV/fm³ [36, 37].

Finally, the analytic expression related to energy density and radial pressure is [35, 38]

where the constant coefficients λ are described by

Note that $ \text{sgn}(\lambda) $ represents the sign of λ. The focus of this investigation is solely on the positive λ space [35].

As demonstrated in Ref. [35], the $ B_{\rm eff} $ parameter can be avoided using the following dimensionless rescaling:

and

so that the EoS (24) reduces to the dimensionless form

The parameter $ \bar{\lambda} $ in the above rescaled EoS directly influences the stiffness of the EoS. Larger values of $ \bar{\lambda} $ lead to a stiffer EoS by increasing the radial pressure at a given energy density. This increased stiffness enables the QSs to support more mass and, as demonstrated in Fig. 2 and Table 2, increases both the maximum mass (from about $ 2.301\,M_\odot $ to $ 2.545\,M_\odot $) and radius. Consequently, the compactness parameter $ M/R $ is slightly enhanced (rising from 0.287 to 0.294) while remaining within the Buchdahl limit. This behavior reflects how the interaction strength encoded in $ \bar{\lambda} $ plays a crucial role in determining the macroscopic properties of QSs.

Figure 2.  (color online) Mass-radius $ (M-R) $ and mass-compactness $ (M-M/R) $relations obtained from the interacting quark matter (IQM) EoS (28) admixed with quasi-local EoS (30). For the given central density, listed in Table 2, we plot all diagrams for the given parameter set: $ \beta = 0.35 $, $ B = 60 $ MeV/fm³, and $ a = 50.0 $ km² with $ \bar{\lambda} \in \{0.05, 0.10, 0.20 \} $. The presented observational constraints are the same as those of Fig. 1. The $ (M-R) $ diagram also includes mass measurements from the GW190425 event.

$ \bar{\lambda} $	$M_{\max}$ $/\rm M_{\odot}$	R/km	$ \rho_c $/(MeV/fm3)	$ M/R $
0.05	2.301	11.871	1,024	0.287
0.10	2.402	12.329	923	0.289
0.20	2.545	12.834	873	0.294

Table 2.  Fundamental features of QSs in $ R^2 $ gravity are presented. Parameters of the considered QSs are $ B = 60 $ MeV/fm³, $ \beta = 0.35 $, and $ a = 50.0 $ km² with $\bar{\lambda} \in \{0.05,~ 0.10, $$ 0.20 \}$.

Notice that, as $ \bar{\lambda} \to 0 $, the aforementioned expression represents the rescaled conventional noninteracting quark matter EoS, where $\bar{P}_r=\dfrac{1}{3}(\bar{\rho}-1)$. Moreover, for sufficiently large values of $ \bar{\lambda} $, that is, when $ \bar{\lambda}>>0 $, we obtain

or, equivalently, $ P_r={\rho}-2B_{\rm eff} $, using Eq. (26). In Refs. [35, 38, 39], authors have shown that the EoS (26) and a much wider range of $ \lambda>0 $ lead to more massive QSs that meet the 2$ M_{\odot} $ restriction.

Here, we examine the anisotropic QSs and particularly explore how we can incorporate the anisotropy as an extra matter component in the modified TOV equations in addition to the interacting quark matter (IQM) EoS. Starting with the quasi-local EoS [40], which has a broad range of applications in both modified gravity theory and GR, the explicit form of the EoS becomes

where β denotes the dimensionless free parameter that measures anisotropy. This Ansatz is motivated by the physical conditions expected in ultra-dense compact stars. In this expression, the factor $ P_r $ ensures that anisotropy scales with the local pressure and vanishes at the surface (where $ P_r=0 $), while the compactness parameter $ \zeta=2m(r)/r $ captures the enhanced gravitational effects toward the star's core. Positive values of β provide additional outward pressure that can counteract gravitational collapse, thereby enhancing stability; in contrast, negative values reduce the tangential pressure, potentially affecting stability in alternative ways. Our subsequent stability analysis confirms that this form of anisotropy keeps both the adiabatic indices above the critical threshold of $ 4/3 $ and preserves causality throughout the stellar configuration.

Although the theoretical range for β can extend from $ -2 $ to $ 2 $, we restrict our study to $ \beta \in [-0.35, 0.35] $ based on both physical and numerical considerations [41−48]. Our investigations indicate that values with $ |\beta|>0.35 $tend to push the solutions toward the limits of physical acceptability-positive values and move the configuration extremely close to the Buchdahl limit, while negative values may compromise stability by reducing the adiabatic indices below the critical threshold. This moderate range of β is also more consistent with current microphysical models and observational constraints, ensuring that our quark star models remain realistic.

Additionally, the compactness is denoted as $ \zeta \equiv 2m(r)/r $, where $ m(r) $ is the gravitational mass bounded by radius r. As $ r\to 0 $, the anisotropy of pressures approaches zero, that is, $ \Delta \to 0 $, resulting in the recovery of the isotropic solution. In [40], it was shown that two pressure components vanish at the surface of the star, that is, $ P_r \left( r \rightarrow R \right) = P_{\perp} \left( r \rightarrow R \right) = 0 $. Note that $ \beta = 0 $ maintains the regularity condition in the stellar interior as well.

Finally, we can now solve the TOV equations (Eqs. (14)−(17)) within the scenarios of interacting quark matter (IQM) EoS (28) admixed with quasi-local EoS (30). We perform a series of numerical calculation with different combinations of three parameters $ (a,\; \beta, \; , \bar{\lambda}) $ and examine the impact of these parameters on the mass-radius $ (M-R) $ relationships. The conventional $ (M, R) $ relation has been reconstructed by introducing $ (M, R) = (\bar{M}/\sqrt{4B_{\rm eff}}, \bar{R}/\sqrt{4B_{\rm eff}}) $. In the following, the mass of the star is measured in solar mass $ (M_{\odot}) $, radius is measured in $\rm km$, energy density and pressures are measured in $\rm MeV/fm^3$, and parameter a is measured in$ \text{ km}^2 $. The scalar field and compactness factor are dimensionless.

To describe mass-radius $ (M-R) $ and mass-compactness $ (M-M/R) $ relations in Fig. 1, we vary the parameter $ a \in [5.0, 50,500] $ in km². To cover the parameter space of our model, we choose the other parameters as $ B = 60 $ MeV/fm³, $ \bar{\lambda} = 0.2 $, and $ \beta = -0.35 $. As one can see, the $ (M-R) $ curves are almost indistinguishable for a ranging from $ a = 5 $ to $ a = 500 $ km². This situation is similar to that observed in Ref. [49], which supports the existence of anisotropic QSs with perturbative QCD corrections. Table 1 also includes a list of the properties of QSs, including their maximum masses and their corresponding radii. The table shows that the maximum mass for QSs can be larger than 2$ M_{\odot} $ and goes up to 2.25 $ M_\odot $ at $ a = 5.0 $ km². The orange and light purple areas correspond to the constraints on the mass and radius measurement derived from the astrophysical observation data provided by the NICER measurements of the pulsars PSR J0030+0451 [50] and PSR J0740+ 6620 [51]. Moreover, the solid blue contour regions represent the GW170817 event [52], whereas the estimation for the mass and radius of the HESS J1731-347 [53] are shown in dark and light pink. Additionally, we demonstrate the impact of model parameter a on the parameters of maximum compactness in the lower panel of Fig. 1. Notice that the trend of the $ (M-M/R) $ curves is similar to that of $ (M-R) $. This means that $ (M-M/R) $ curves are almost indistinguishable from each other, and its value goes up to $ 0.266 $, as evident from Table 1. We note that the minimal differences observed in the mass-radius relation for varying values of a (from 5 to 500 km²) arise because the $ R^2 $ correction mainly influences the high-curvature core region. In particular, the corresponding scalar field decays rapidly outside this core, and the effective potential,

Figure 1.  (color online) Mass-radius $ (M-R) $ and mass-compactness $ (M-M/R) $ relations obtained within the scenarios of interacting quark matter (IQM) EoS (28) admixed with quasi-local EoS (30). For the given central density, listed in Table 1, we plot all diagrams for the given parameter set: $ B = 60 $ MeV/fm³, $ \bar{\lambda} = 0.2 $, and $ \beta = -0.35 $ with $a \in [5,~ 50, ~ 500]$ km². The shaded areas in the diagram represent observation constraints from analysis of NICER observations the pulsars PSR J0030+0451 [50] and PSR J0740+6620 [51]. The solid blue contour regions stand for the GW170817 event [52], whereas the estimation for the mass and radius of the HESS J1731-347 [53] are shown in dark and light pink.

a/km2	$M_{\max}$ $/M_{\odot}$	R/km	$ \rho_c $/(MeV/fm3)	$ M/R $
5	2.247	12.470	923	0.267
50	2.239	12.456	974	0.266
500	2.237	12.463	974	0.266

Table 1.  Fundamental features of QSs in $ R^2 $ gravity are presented. Parameters of the considered QSs are $ B = 60 $ MeV/fm³, $ \bar{\lambda} = 0.2 $, and $ \beta = -0.35 $ with $ a \in [5, 50,500] $ km².

becomes flatter for larger a, thereby limiting its impact on the overall stellar structure. Moreover, the imposed boundary conditions further constrain the scalar field's influence, making the macroscopic properties, such as mass and radius, largely insensitive to a. This behavior appears to be a generic feature of $ R^2 $ gravity in the context of compact stars, suggesting that more sensitive observables (for example, tidal deformability) may be needed to distinguish between different values of a. Moreover, the table also indicates that the Buchdahl limit remains intact, that is, $ M/R < 4/9 $ [54].

We further conduct the analysis by varying $ \bar{\lambda} $with $ \beta = 0.35 $, $ B = 60 $ MeV/fm³, and $ a = 50.0 $ km². In Fig. 2, we show the $ (M-R) $ curves obtained for different values of $ \bar{\lambda} $ and the main QS properties summarized in Table 2. The maximum masses and their corresponding radii are increasing in conjunction with the increasing values of $ \bar{\lambda} $. In this case, the maximum gravitational mass and its radius are found to be 2.30−2.55 and 11.87−12.83 km, respectively, for the considered parameter sets (refer Table 2), which is consistent with GW190814 data. Moreover, in the lower panel of Fig. 2, we also illustrate the impact of $ \bar{\lambda} $ on the properties of maximal compactness. Observe that the maximal compactness increases as $ \bar{\lambda} $increases, typically falling within the range of 0.287 to 0.294, as also summarized in Table 1. The table further demonstrates that the Buchdahl limit is preserved, that is, $ M/R < 4/9 $ [54].

Finally, we examine the physical effect of anisotropy parameter β on the properties of QSs. In Fig. 3, we plot the $ (M-R) $ and $ (M-M/R) $ curves for$ \beta \in [-1, 1] $, $ B = 60 $ MeV/fm³, $ \bar{\lambda} = 0.2 $, and $ a = 50.0 $ km². The figure shows that the $ (M-R) $curves vary significantly when the anisotropic pressure is varied, that is, when we vary β. According to our calculations, the maximum mass and radius of QSs is more than 2$ M_{\odot} $and goes up to 2.55 $ M_\odot $ at $ \beta = 0.35 $. This illustrates that the maximum mass of the QS increases for increasing values β. Simultaneously, we measured the QS's maximum mass for the isotropic scenario, or $ \beta = 0 $, and observed a decrease to 2.34 $ M_\odot $; see Table 3. This is a clear indication that the presence of the anisotropic pressure leads to the existence of a massive QS in $ R^2 $ gravity. Moreover, the theoretically attainable maximum mass is compatible with the mass-radius measurements obtained from different astrophysical observations. We further illustrate the $ (M-M/R) $ diagram in the lower panel of Fig. 3. Notice that the $ (M-M/R) $ curves follows the same pattern as the $ (M-R) $ relations, that is, the maximum value of compactness increases with β. This value lies within the range of $ 0.266 < M/R < 0.294 $. We summarize our results for various values of β in Table 3. Interestingly, the compactness does not exceed the Buchdahl bound $ M/R < 4/9 $ [54].

Figure 3.  (color online) Mass-radius $ (M-R) $ and mass-compactness $ (M-M/R) $ relations obtained within the scenarios of interacting quark matter (IQM) EoS (28) admixed with quasi-local EoS (30). For the given central density, listed in Table 3, we plot all diagrams for the given parameter set: $ B = 60 $ MeV/fm³, $ \bar{\lambda} = 0.2 $, and $ a = 50.0 $ km² with $\beta \in \{-0.35,~ 0.00,~ 0.35\}$. The presented observational constraints are the same as those of Fig. 1. The green horizontal bar represents the maximum mass constraints from the GW190814 event.

β	$M_{\max}$ $/M_{\odot}$	R/km	$ \rho_c $/(MeV/fm3)	$ M/R $
-0.35	2.239	12.456	973	0.266
0.00	2.389	12.653	923	0.280
0.35	2.545	12.834	873	0.294

Table 3.  Fundamental features of QSs in $ R^2 $ gravity are presented. Parameters of the considered QSs are $ B = 60 $ MeV/fm³, $ \bar{\lambda} = 0.2 $, and $ a = 50.0 $ km² with $\beta \in \{-0.35,~ 0.00, $$ 0.35\}$.

To address the astrophysical feasibility of our proposed model, we assess the stability of anisotropic QSs utilizing the static stability criterion, adiabatic index, and speed of sound. We examine each of these properties with a visual representation.

We now examine the stability of the equilibrium configuration using the static stability criterion (SSC) [55, 56], which is a necessary but insufficient condition. The results are depicted in the $ M-\rho_c $ plane, with M representing the total mass and $ \rho_c $ denoting the central density of the star. Mathematically, these inequalities are as follows:

In Fig. 4, we plot the variation of $ M-\rho_c $ curves for the cases mentioned above. We utilize the parameters used in Figs. 1 and 3. The figures demonstrate that the total gravitational mass increases with growing central energy density, reaching a point where $ (M_{\text{max}}, R_{M_{\text{max}}}) $ is defined. Besides, the opposite inequality ${\rm d}M/{\rm d}\rho_c < 0$ implies instability with respect to small deformations. Moreover, Fig. 4 shows that increasing the anisotropy parameter β significantly modifies the stability boundary in the M-$ \rho_c $ plane. In particular, the critical central density where ${\rm d}M/{\rm d}\rho_c=0$ shifts to lower values and the maximum stable mass increases (e.g., from $ M_{\max}\approx2.24\,M_\odot $ at $ \beta=-0.35 $ to $ M_{\max}\approx2.55\,M_\odot $ at $ \beta=0.35 $). This behavior arises because positive β enhances tangential pressure, providing additional support against gravitational collapse, and leading to a steeper transition between stable and unstable configurations. These observations are consistent with our adiabatic index and sound speed analyses, confirming that anisotropy plays a crucial role in stabilizing the QS configurations.

Figure 4.  (color online) Relationships between mass and central density. The parameters used are the same as in Figs. 1 to 3.

The adiabatic index γ is a key parameter in determining the dynamical stability of compact stars, which was first proposed by Chandrasekhar [57]. For stability, the adiabatic index must satisfy the condition $\gamma > {4}/{3}$, which is critical for preventing gravitational collapse [58]. The adiabatic indices in the radial ($ \Gamma_r $) and tangential ($ \Gamma_t $) directions are crucial in the case of anisotropic stars. They are defined as follows [59−63]:

where $ \sigma = P_{\perp} - P_r $ is the anisotropy factor, and $ p_r' $ is the derivative of the radial pressure. The stability condition depends on the anisotropy, as pressure behaves differently in radial and tangential directions. Therefore, we also define the following adiabatic indices for each direction:

where $ v_r $ and $ v_t $ are the speeds of sound in the radial and tangential directions, respectively, ρ is the energy density, $ p_r $ is the radial pressure, and $ p_t $ is the tangential pressure.

For isotropic stars, the adiabatic index γ typically satisfies the stability criterion $\gamma > {4}/{3}$. However, in anisotropic stars, the condition $\gamma > {4}/{3}$ must be satisfied separately in both radial and tangential directions, which is why $ \Gamma_r $ and $ \Gamma_t $ become essential parameters in determining stability.

According to Chandrasekhar's approach, the condition for stability is that both adiabatic indices $ \Gamma_r $ and $ \Gamma_t $ must exceed $ {4}/{3} $ to ensure that the star is dynamically stable in both directions. In particular, anisotropy can enhance stability by increasing $ \Gamma_t $, thereby providing additional resistance to collapse in the tangential direction.

Our results for the adiabatic indices in radial and tangential directions are shown in Fig. 5, which indicate that the considered model is dynamically stable. Our analysis indicates that anisotropic effects on dynamical stability are parameter dependent. Moderate positive anisotropy (e.g., $ 0<\beta<0.35 $) generally enhances stability by keeping both $ \Gamma_r $ and $ \Gamma_t $ safely above the critical value of $ 4/3 $, whereas strong negative anisotropy (e.g., $ \beta < -0.35 $) may lead to instability in certain regions, particularly in the outer layers of the star.

Figure 5.  (color online) Adiabatic indices in radial ($ \Gamma_r $) and tangential ($ \Gamma_{\perp} $) directions with the same parameters as those of Figs. 1 to 3.

The speed of sound $ v_s $ is an important criterion for assessing both the stability and physical validity of stellar configurations. In general relativity, the sound speed must respect causality, meaning it must not exceed the speed of light (c). This imposes the condition

where $ v_s^2 $ is the square of the sound speed, P is the pressure, and ρ is the energy density. The inequality ensures that perturbations in the star propagate at subluminal speeds, preserving causality, and that the pressure is a monotonically increasing function of the energy density.

In our analysis, we compute the sound speed in both the radial ($ v_r^2 $) and tangential ($ v_t^2 $) directions to account for possible anisotropy in the stellar model. For the radial direction, the sound speed is defined as

Similarly, the tangential sound speed is given by

For anisotropic stars, both $ v_r^2 $ and $ v_{\perp}^2 $ must independently satisfy the causality condition.

Figures 6, 7, and 8 illustrate the radial and tangential speeds of sound for different parameter sets. Our results confirm that the speeds of sound in both directions remain within the physically acceptable range, $ 0 \leq v_s^2 \leq 1 $, thus satisfying the causality condition across all models considered. Additionally, we find that the variation in the speed of sound is smooth, indicating that the stellar configurations are stable and consistent with the physical requirements for compact stars.

Figure 6.  (color online) Square of the speed of sound plotted in both radial and tangential directions for the parameter sets defined in Fig. 1.

Figure 7.  (color online) Square of the speed of sound plotted in both radial and tangential directions for the parameter sets defined in Fig. 2.

Figure 8.  (color online) Square of the speed of sound plotted in both radial and tangential directions for the parameter sets defined in Fig. 3.

Therefore, the computed speeds of sound validate the stability and causal nature of our anisotropic QS models within the framework of $ R^2 $ gravity. These results also support the conclusion that the presence of anisotropy, characterized by different pressures in the radial and tangential directions, does not lead to violations of causality or instabilities in the considered QS configurations.

Due to the study of gravitational waves, physicists can investigate binary neutron or QSs via tidal effects. We show the relation between the tidal Love number ($ k_2 $) and tidal deformability (λ) for the anisotropic QSs in the Starobinsky model. QSs are massive within a limited volume, leading to one of the highly densed objects in our Universe. They form themselves against the tidal deformability λ defined as

where R is the radius of the stars. $ k_2 $ is the dimensionless parameter, which is sensitive to the EoS and plays a major role in tidal forces. To obtain the Love number $ k_2 $, we can solve the following equation written in terms of compactness $ C = M/R $

where $ y_R = y(R) $ is at the radius R and determined by

where

We also introduce the dimensionless tidal deformability

As Λ decreases, the stars become more resilient against tidal deformation.

We also plot the dimensionless tidal deformability Λ against the total mass $ M_{\odot} $ in Fig. 9 for the parameter set defined in Figs. 1 to 3. We also compare our results to the observational result of GW170817 as the constraint on the tidal deformability. Anisotropy is expected to imprint observable signatures on the gravitational wave spectrum of QSs. In particular, changes in the tidal deformability parameter Λ (as illustrated in Fig. 9) can alter the inspiral waveform during binary mergers, while shifts in the oscillation mode frequencies (e.g., f-modes and p-modes, see [64, 65]) due to increased compactness for positive β may lead to higher-frequency post-merger signals.

Figure 9.  (color online) Tidal deformation of the anisotropic QSs in the Starobinsky model plotted for the parameter sets defined in Figs. 1, 2, and 3, respectively.

In this study, we investigated the structure and stability of QSs within the framework of $ R^2 $ gravity, specifically using the Starobinsky model, where $ f(R) = R + a R^2 $. We explored the impact of anisotropic pressures and the quark matter EoS with $ \mathcal{O}(m_s^4) $ corrections on the mass-radius ($ M-R $) and mass-central density ($ M-\epsilon_c $) relations for QSs. By solving the modified TOV equations, we examined the role of the $ R^2 $ gravity parameter a, as well as the anisotropy parameters $ B_{\perp} $ and $ a_4^{\perp} $ derived from the EoS, on the properties and stability of these compact objects.

We found that the presence of anisotropy in the pressure significantly affected the mass and radius of QSs. In particular, our results indicated that increasing the anisotropy parameter β enhanced the maximum mass of the QSs, suggesting the possibility of forming super-massive pulsars, consistent with recent astronomical observations, such as PSR J0952-0607 and GW170817. The calculated mass-radius relationships and compactness of the QSs adhered to the Buchdahl limit, ensuring that the solutions remained physically viable. We also confirmed that the inclusion of anisotropic pressures does not violate causality, as the speeds of sound in both radial and tangential directions satisfied the condition $ 0 \leq v_s^2 \leq 1 $.

Moreover, we analyzed the dynamical stability of the QSs by computing the adiabatic indices in the radial and tangential directions. Our results showed that the adiabatic indices $ \Gamma_r $ and $ \Gamma_t $ exceeded the critical value of $ 4/3 $, thereby ensuring the stability of the QSs against radial perturbations. We also applied the static stability criterion, verifying that the gravitational mass increased monotonically with central density, which further supported the stability of the configurations.

In summary, our study provided a detailed analysis of QSs in $ R^2 $ gravity, incorporating both pressure anisotropy and interacting quark matter EoS. The results demonstrated that QSs in $ R^2 $ gravity could attain masses greater than $ 2M_{\odot} $, with stability and causality preserved. These findings contribute to the growing body of evidence supporting the existence of super-massive QSs, as inferred from recent astrophysical data.

While the results of this study are promising, several questions remain open. One potential direction for future research is to explore the impact of rotation on the structure and stability of anisotropic QSs in modified gravity, as rotation could significantly alter the maximum mass and stability conditions. Additionally, studying the behavior of QSs under more complex equations of state, particularly those including color superconducting phases, may provide deeper insights into the internal structure of these compact objects. Finally, the gravitational wave signatures from anisotropic QSs in $ R^2 $ gravity could be explored further, with the aim of distinguishing between different compact object models through future gravitational wave observations.

In this Appendix, we outline the derivation of the field Eqs. (14)−(17) starting from the Einstein frame action

By varying the action with respect to the metric $ \tilde{g}_{\mu\nu} $, one obtains the Einstein equations. For the static, spherically symmetric metric given in Eq. (13), the explicit forms of the components of the Einstein tensor (namely, $ \tilde{G}_{tt} $, $ \tilde{G}_{rr} $, and $ \tilde{G}_{\theta\theta} $) are derived accordingly:

The total energy-momentum tensor is constructed as the sum of the scalar field contribution,

and that of the anisotropic fluid, which transforms according to

with the Einstein frame matter variables related as

Furthermore, by varying the action with respect to φ, one obtains the scalar field equation

where the coupling constant is defined as

Finally, using the Bianchi identities in conjunction with the static metric (13), the conservation law yields the modified TOV equation. In this derivation, the anisotropic term $\dfrac{2}{r}(P_{\perp}-P_r)$ naturally emerges. For further details, see [30, 31].

The authors express gratitude for the editor's and reviewers' insightful feedback, which has substantially enhanced this manuscript's quality and clarity. A. Pradhan expresses gratitude to the IUCCA in Pune, India, for offering facilities under associateship programs. İ. S. expresses gratitude to EMU, TÜBİTAK, ANKOS, and SCOAP3 for their academic and/or financial support. İ. S. and T. T. also acknowledge COST Actions CA22113, CA21106, and CA23130 for their contributions to networking.