The effect of pressure anisotropy on quark stars structure in the Starobinsky model

The study of compact objects, such as QSs, has driven extensive research efforts due to their unique role as astrophysical laboratories for high-density matter and potential indicators of physics beyond general relativity (GR). Numerous modified gravity theories 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 not only addresses early-universe inflation but is also 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, thus making them attractive for modeling compact stars like QSs [15−23]. The gravitational wave signal from GW 190814 has initiated the possibility on the secondary object as a neutron star with a mass 2.6 $ M_{\odot} $. The modified gravity theories have been proposed on this phenomena [22, 23]

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 interactions among particles at ultra-high densities [24, 25]. 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 [16, 26]. 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 [27].

This paper 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 the modified TOV [28] 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.

The plan of the paper is 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 Lagrangian is replaced by a general function of the Ricci scalar, and the action is given by

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

respectively. For our purposes, we shall consider a specific form of Starobinsky's model, i.e. $ f(R)=R+aR^2 $, where a is taken to be $ a \geq 0 $ in agreement with the inequalities given in Eq. (2) [29]. 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 again.

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

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

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

Taking variation with respect to the metric $ \tilde{g}_{\mu \nu} $ and the 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{dV(\phi)}{d\phi} $. Moreover, the coupling constant $ \alpha (\varphi) $ is related with [32, 33]

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, which is described by

where ρ is the energy density, $ P_{\perp} $ and $ P_r $ are the pressures perpendicular and parallel to the spacelike vector $ \chi_\nu $. 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 [32, 33]

where the tilde indicates the Einstein frame.

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

where $ d \Omega^2 = d\theta^2 + \sin^2\theta 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. In order to provide a comprehensive description of the configuration that is being examined, 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 the exterior problem simultaneously. Thus, we need boundary conditions in the natural Einstein frame to maintain regularity at the origin,

where $ \rho_c $ and $ \phi_c $ are the central value of energy density and the 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 [32], the condition $ \frac{d\varphi}{dr}(0) = 0 $ guarantees the regularity of the scalar field ϕ, which consequently secures the regularity of Φ at the stellar interior $ r = 0 $. Moreover, $ \Lambda(0)=0 $ is necessary for the Einstein frame geometry at the center to be 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. On the other hand, we choose the boundary conditions in such a manner that they could 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 (20) guarantee asymptotic flatness in both the Einstein and Jordan frames.

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

Next, we will show the structure equation for a QCD-driven EoS and 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 [34]. This means that this method could be useful for studying stars that are very densely packed together.

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

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 $ 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 [36]. In the above expression, the term with $ m_s $ 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³ [38, 39].

Finally, the analytic expression related to energy density and radial pressure is [37, 40]

where the constant coefficients λ described by

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

As demonstrated in Ref. [37], the $ B_{\rm eff} $ parameter can be avoided using the subsequent 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 to be seen in Fig. 2 and Table 2, results in an increase in 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 still 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 are obtained within the scenarios of interacting quark matter (IQM) EoS (28) admixed with quasi-local EoS (30). For the given central density, shown in Table 2, we plot all the diagrams for the given parameter set: $ \beta = 0.35 $, $ B = 60 $ MeV/fm³, $ a = 50.0 $ km² in variation of $ \bar{\lambda} \in \{0.05, 0.10, 0.20 \} $. The presented observational constraints are same as of Fig. 1. The $ (M-R) $ diagram also includes mass measurements from the GW190425 event. [54].

$ \bar{\lambda} $	$ M_{max} $ $ (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.  The fundamental features of QSs in $ R^2 $ gravity are presented. Parameters of the considered QSs are $ B = 60 $ MeV/fm³, $ \beta = 0.35 $, $ a = 50.0 $ km² in variation of $ \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=\frac{1}{3}(\bar{\rho}-1) $. Moreover, sufficiently large value of $ \bar{\lambda} $ i.e., when $ \bar{\lambda}>>0 $, we get

or, equivalently, $ P_r={\rho}-2B_{\rm eff} $, using Eq. (26). In Refs. [37, 40, 41], 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 with the interacting quark matter (IQM) EoS. Starting with the quasi-local EoS [42], that has a broad range of applications in both modified gravity theory and GR. The explicit form of the EoS,

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, whereas negative values reduce the tangential pressure, potentially affecting stability in alternative ways. Our subsequent stability analysis confirms that this form of anisotropy maintains 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 [43−50]. Our investigations indicate that values with $ |\beta|>0.35 $ tend to push the solutions toward the limits of physical acceptability-positive values bring the configuration uncomfortably close to the Buchdahl limit while negative values may compromise the 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, i.e., $ \Delta \to 0 $, resulting in the recovery of the isotropic solution. In [42] it was shown that two pressure components vanish at the surface of the star, i.e., $ P_r \left( r \rightarrow R \right) = P_{\perp} \left( r \rightarrow R \right) = 0 $. Note that $ \beta = 0 $ maintains the regularity condition at the stellar interior also.

Finally, we are in a position to solve the TOV equations (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, we measure the mass of the star in solar mass $ (M_{\odot}) $, radius is measured in $ km $, the energy density and the pressures in $ MeV/fm^3 $, and the parameter a in $ \text{ km}^2 $, whereas the dimensions of the scalar field and the factor of compactness 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 $, respectively. As one can see, the $ (M-R) $ curves are almost indistinguishable from each other for different values of the parameter a ranging from $ a = 5 $ to $ a = 500 $ km². This situation is similarly to Ref. [51], which supports the existence of anisotropic QSs with an integrating perturbative QCD corrections. Table 1 also includes a list of the properties of QSs, including their maximum masses and their corresponding radii. From the table, it is seen 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² in our investigation. 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 [52] and PSR J0740+6620 [53]. Moreover, the solid blue contour regions stand for the GW170817 event [55], whereas the estimation for the mass and radius of the HESS J1731-347 [56] are shown in dark and light pink. Beside that, we also 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 are obtained within the scenarios of interacting quark matter (IQM) EoS (28) admixed with quasi-local EoS (30). For the given central density, shown in Table 1, we plot all the diagrams for the given parameter set: $ B = 60 $ MeV/fm³, $ \bar{\lambda} = 0.2 $, $ \beta = -0.35 $ in variation of $ 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 [52] and PSR J0740+6620 [53]. The solid blue contour regions stand for the GW170817 event [55], whereas the estimation for the mass and radius of the HESS J1731-347 [56] 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.  The fundamental features of QSs in $ R^2 $ gravity are presented. Parameters of the considered QSs are $ B = 60 $ MeV/fm³, $ \bar{\lambda} = 0.2 $, $ \beta = -0.35 $ in variation of $ 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 (e.g., tidal deformability) may be needed to distinguish between different values of a. Moreover, the table also indicates that the Buchdahl limit remains intact, i.e., $ M/R < 4/9 $ [57].

We further carry out the analysis in variation of $ \bar{\lambda} $ with other parameters for our numerical calculations: $ \beta = 0.35 $, $ B = 60 $ MeV/fm³ and $ a = 50.0 $ km², respectively. 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. It is evident that the tendencies of 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 their 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 found in Table 1. The table further demonstrates that the Buchdahl limit is preserved, i.e., $ M/R < 4/9 $ [57].

Finally, we have examined the physical effect of anisotropy parameter β on the properties of QSs. In the Fig. 3 we plot the $ (M-R) $ and $ (M-M/R) $ curves in variation of $ \beta \in [-1, 1] $ with other parameters for our numerical calculations: $ B = 60 $ MeV/fm³, $ \bar{\lambda} = 0.2 $ and $ a = 50.0 $ km², respectively. It is clearly seen from the figure that the $ (M-R) $ curves vary significantly when the anisotropic pressure is varied, i.e., when we vary the parameter β. In principle, the maximum mass and radius of QSs more than 2$ M_{\odot} $, and goes up to 2.55 $ M_\odot $ at $ \beta = 0.35 $ from our calculation. This illustrates that there the maximum mass of the QS increases for increasing values β. Simultaneously, we also measured the QS's maximum mass for the isotropic scenario, or $ \beta = 0 $, which decreases 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 comes out to be 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 trend of $ (M-M/R) $ curves follows the same pattern as $ (M-R) $ relations; i.e., the maximum value of compactness increases when β grows. This value lies within the range of $ 0.266 < M/R < 0.294 $. We summarize our results for various values of the parameter β in Table 3. Interestingly, the compactness does not exceed the Buchdahl bound, $ M/R < 4/9 $ [57].

Figure 3.  (color online) Mass-radius $ (M-R) $ and mass-compactness $ (M-M/R) $ relations are obtained within the scenarios of interacting quark matter (IQM) EoS (28) admixed with quasi-local EoS (30). For the given central density, shown in Table 3, we plot all the diagrams for the given parameter set: $ B = 60 $ MeV/fm³, $ \bar{\lambda} = 0.2 $, $ a = 50.0 $ km² in variation of $ \beta \in \{-0.35, 0.00, 0.35\} $. The presented observational constraints are same as of Fig. 1. The green horizontal bar represents the maximum mass constraints from 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.  The fundamental features of QSs in $ R^2 $ gravity are presented. Parameters of the considered QSs are $ B = 60 $ MeV/fm³, $ \bar{\lambda} = 0.2 $, $ a = 50.0 $ km² in variation of $ \beta \in \{-0.35, 0.00, 0.35\} $.

We addressed the astrophysical feasibility of our proposed model, we will assess the stability of anisotropic QSs utilizing the static stability criterion, adiabatic index, and sound speed. We shall examine each of them sequentially with a visual representation.

We will now examine the stability of the equilibrium configuration using the static stability criterion (SSC) [58, 59]. It is a necessary condition, but not a sufficient one. 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 all considered cases mentioned above. We utilize the same parameters for plotting as illustrated in Figs. 1 and 3. It is clear from those figures 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 $ dM/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 $ dM/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 the anisotropy plays a crucial role in stabilizing the QS configurations.

Figure 4.  (color online) The profiles illustrate the 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 [60]. For stability, the adiabatic index must satisfy the condition $ \gamma > \frac{4}{3} $, which is critical for preventing gravitational collapse [61]. 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: [62−66]

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 sound speeds 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 > \frac{4}{3} $. However, in anisotropic stars, the condition $ \gamma > \frac{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 the stability.

According to Chandrasekhar's approach, the condition for stability is that both adiabatic indices $ \Gamma_r $ and $ \Gamma_t $ must exceed $ \frac{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 Figure 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) The adiabatic indices in radial ($ \Gamma_r $) and tangential ($ \Gamma_{\perp} $) directions with the same parameter set in 1 to 3.

The sound speed $ 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:

In the case of 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 sound speeds for different parameter sets. Our results confirm that the sound speeds 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 sound speed is smooth, indicating that the stellar configurations are stable and consistent with the physical requirements for compact stars.

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

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

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

Therefore, the computed sound speeds 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 QS configurations considered.

Due to the study of gravitational waves, physicists can investigate binary neutron or QSs via tidal effects. We will show the relation between the tidal Love number ($ k_2 $) and the 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 EoS and plays a major role to the 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 is determined by

where

We also introduce the dimensionless tidal deformability

The lower of the Λ value becomes, the higher of the stars could hold against the tidal deformation.

We also plot the dimensionless tidal deformability Λ against the total mass $ M_{\odot} $ in Fig. 9 for the parameter set defined in Fig. 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 [67, 68]) due to increased compactness for positive β may lead to higher-frequency post-merger signals.

Figure 9.  (color online) The tidal deformation of the anisotropic QSs in Starobinsky model has been plotted for the parameter sets defined in Fig. 1, 2 and 3, respectively.

In this work, 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 quark matter equations of state (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 sound speeds 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 open questions remain. 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 $ \frac{2}{r}(P_{\perp}-P_r) $ naturally emerges. For further details, the reader is referred to [32, 33].

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.