Effects of anisotropic pressure on interacting quark star structure

Perturbative Quantum Chromodynamics (pQCD) corrections and color superconductivity predict that strongly interacting matter can reveal new physical phenomena under extreme conditions. Taking into account these interaction effects, we investigate the role of anisotropic pressure in quark stars composed of interacting quark matter. Adopting two physically well-motivated anisotropy profiles, we numerically solve the stellar structure equations in order to explore the consequences of anisotropic pressure on various macroscopic properties such as radius, gravitational mass, surface redshift, moment of inertia, tidal Love number and oscillation spectrum. Remarkably, for both anisotropy models, negative anisotropies increase the radial stability of interacting quark stars, while the opposite occurs for positive anisotropies. However, for the Bowers-Liang profile, the central density corresponding to the maximum-mass point does not coincide with the central density where the squared oscillation frequency vanishes, indicating that the existence of stable anisotropic interacting quark stars is possible beyond the maximum mass for negative anisotropies. Additionally, we compare our theoretical predictions with several observational mass-radius measurements and tidal deformability constraints, which suggest that both strong interaction effects and anisotropy effects play a crucial role in describing compact stars observed in the Universe.


I. INTRODUCTION
Einstein's theory of general relativity (GR) has stood like a pillar of modern theoretical physics regarding gravitational phenomena from the solar system to cosmological scales [1].The success of this theory comes through the first experimental test by Arthur Eddington in 1919 during a total solar eclipse, i.e., the deflection of the starlight next to the Sun.Since then GR is the simplest metric theory and it remains to be the most successful gravity theory for understanding the Universe.Among many astonishing predictions of GR, compact astrophysical objects such as black holes, neutron stars (NSs), white dwarfs have turned from purely mathematical objects to potentially real physical entities as a result of different observational astronomical measurements.
In 1967, the discovery of pulsars had a great impact on astronomers in general.This discovery proves the existence of NSs in the Universe and is crucial to understand the nature of ultra-dense compact objects.NSs are the incredibly dense remnants of high mass stars when they run out of fuel.Once the core of the massive stars has completely burned to iron, the energy production stops and the core rapidly collapses, squeezing electrons and protons together to form neutrons and neutrinos.Such stars are supported by neutrondegeneracy pressure, and thus, are the most compact stars in the Universe.A NS that has a mass between M ∼ 1 − 3M ⊙ where M ⊙ = 2 × 10 33 g with radius between 10 − 15 km [2,3].Consequently, their central densities are extremely high and easily exceed the nuclear saturation limit, i.e., ρ ≳ ρ nuc where ρ nuc = 2.8 × 10 14 g/cm 3 .Indeed, it is very hard to deal with dense matter at such extreme situation in a laboratory conducted on Earth, and hence no comprehensive picture has been achieved till date.Moreover, observed pulsars through electromagnetic (EM) signals have put a strong constraint on the equa-arXiv:2311.18770v2[gr-qc] 14 Dec 2023 tion of state (EoS) of dense matter in the interior of NSs.Meanwhile, the mass-radius measurements from spectroscopic observations of thermonuclear X-ray bursts, along with recent NICER (Neutron Star Interior Composition Explorer) data have significantly placed even tighter constraints on the EoS [4].Consequently, physicists have conjectured the existence of more exotic states such as strange quark matter (SQM) in the core of compact objects.As a matter of fact, was first speculated that compact stars could be partially or totally made of SQM [5][6][7].It has been suggested that SQM consists of almost equal numbers of u, d and s quarks, and a small number of electrons to attain the charge neutrality.Quarks are strongly interacting particles and may exist from a few fermis up to a large (kilometer-sized) ranging in size with the possibility of of forming consistent self-bound quark stars (QSs).The simple model proposed for SQM is the MIT bag model [8] in which the quarks are considered to be free inside a bag.Depending on this model, the internal structure of QSs has been explored by several authors in Einstein gravity [9][10][11][12][13][14] and modified gravity theories [15][16][17][18].
Most of the studies concerning the internal structure of NSs/QSs is fully determined based on the assumption of isotropic matter, where the stellar fluid is described by an isotropic perfect fluid.However, the discovery of more massive compact stars has turned out the situation more complex inside a compact stellar object, i.e., taking into consideration a transverse pressure in addition to a radial one [19].Under this anisotropic picture, Bowers & Liang [20] examined the possible effect of locally anisotropic EoS for relativistic spheres.They showed that anisotropy affects the maximum equilibrium mass and surface redshift of NSs.Indeed, it has been realized that local anisotropy in matter fields may arise for several reasons such as superfluid cores, strong magnetic fields, crystallization of the core or phase transitions, which play an important role in determining the microscopic composition of compact stars.But, the isotropy should be eventually restored if one considers an anisotropic matter source for compact stars.
With account of the anisotropy matter field, various solutions have been obtained in the literature [21][22][23][24][25][26][27][28][29].It was argued that the quasi-local variables can construct an EoS to describe the anisotropic fluid in spherical symmetry [30].The effects of anisotropy on the global stellar properties of NSs such as their mass-radius relation, the moment of inertia and tidal deformability were studied in Refs.[31][32][33][34].Within an anisotropic context, the radial stability was rigorously addressed for quark stars [10], neutron stars [31] and dark energy stars [35] through radial and adiabatic perturbations.Based on the Hartle-Thorne formalism, the mass correction and deformation of slowly rotating anisotropic NSs was discussed by Pattersons and Sulaksono [36].See also Refs.[37,38] for slowly rotating anisotropic NSs in modified gravity.In this work we will undertake how the anisotropic pressure affects the internal composition of compact stars with a more general and realistic equation of state for quark matter, namely, an interacting quark EoS.
The rest of the paper is organized as follows.In Sec.II, we preferred to describe the quark matter EoS which includes pQCD corrections and color superconductivity.In Sec.III, we present the gravitational field equations for an anisotropic fluid distribution assuming a static, spherically symmetric geometry.The Tolman-Oppenheimer-Volkoff (TOV) equations are also presented with proper boundary conditions.In the same section, we study various macroscopic properties of anisotropic interacting quark stars (IQSs) such as moment of inertia, tidal Love number and radial vibration frequencies.In Sec.IV, we employ two different anisotropy models for QSs in strong gravitational fields.In next, we continue our discussion for numerical findings specially focusing on the mass-radius relation and the radial stability with varying degrees of anisotropy in Section V. We further compare our theoretical predictions with observational data in Sec.VI.Finally, a summary is provided in Sec.VII.

II. INTERACTING QUARK MATTER EOS
Let us start by reviewing the main results contained in Ref. [39], where authors have established a general parametrization of the quark matter EoS which is specified by the perturbative QCD (pQCD) correction and color superconductivity [39].One interesting feature is that depending only on a single parameter one can convert the EoS into a dimensionless form.Moreover, the rescaling process reduces the number of degrees of freedom of the given EoS.In the current study, we are searching for the presence of anisotropic matter, and how it affects the macroscopic properties of SQSs.
To model the interacting quark matter, the relation between the radial pressure (P r ) to the energy density (ρ) is given by [39,40]: where the bag constant as an effective parameter is denoted by B eff that accounts for the nonperturbative con-tribution from the QCD vacuum and One can immediately check that when λ → 0, the Eq. ( 1) become the conventional noninteracting quark matter EoS P r = (ρ − 4B eff )/3.We can therefore say that the parameter λ is a relative measure of the strong interaction effects.In the following, ∆ denotes the gap parameter and m s represents the strange quark mass.The pQCD correction parameter a 4 can vary from small values to a 4 = 1, which gives a larger deviation from an ideal relativistic Fermi gas EoS.Here, sgn(λ) represents the sign of λ, and leads to a positive value when CFL phase characterizing the possible phases of color superconductivity.We can further introduce the dimensionless rescaling Consequently, removing the B eff parameter, Eq. (1) assumes the following dimensionless form It can be seen that for extremely large positive values of λ, Eq. ( 5) approaches the special form which is equivalent to P r = ρ − 2B eff .In Fig. 1 we show this linear behavior for a very large λ by a blue curve.We further note that the parameter λ gives a stiffer EoS for positive increasing values of λ, where the well known case λ = 0 describes ordinary noninteracting quark matter.With these increasing values one can reach sufficiently high masses [39,40] for IQSs.This can be clearly observed in the mass-radius diagram displayed in Fig. 2 into a dimensionless form.The advantage of using a dimensionless description for the different variables is that it becomes unnecessary to specify  (5) for several values of λ, where the linear behavior for λ → ∞ is described by the blue curve, see Eq. ( 6).Moreover, we recover the standard noninteracting quark matter EoS when λ = 0, which is represented by the black curve.
the value of the effective bag constant, unless we want to obtain the global characteristics of quark stars in physical units for comparison reasons.Furthermore, it was shown by Zhang [40] that it is possible to produce gravitational wave echoes for IQSs when λ ≳ 10 at large central pressure.In the present work, we will consider a value of λ for which the strong interaction effects are relatively appreciable, e.g.see the red curve in Fig. 1.However, in section VI we will adopt smaller values in order to compare our results with observational data.Of course, the qualitative behavior of the results must hold as λ varies.

A. TOV equations
In the present work we assume an anisotropic matter distribution where the radial pressure P r differs from the transverse pressure P ⊥ .The covariant form of the energy-momentum tensor of such source may be written as where u µ is the 4-velocity for a comoving fluid element, χ µ is the unit spacelike vector in the radial direction and g µν is the metric tensor.Moreover, ρ is the energy density while σ ≡ P ⊥ − P r is the anisotropic factor in the source, with σ = 0 corresponding to isotropic matter.The four-vectors u µ and χ µ must satisfy the following relations Considering static and spherically symmetric compact stars, we use the usual line element given by ds 2 = −e 2Φ(r) dt 2 + e 2Ψ(r) dr 2 + r 2 (dθ 2 + sin 2 θdϕ 2 ), (9) where Φ(r) and Ψ(r) are unknown metric potentials to be determined.The gravitational mass enclosed within a sphere of radius r, denoted by m(r), is related to the metric function Ψ(r) through e −2Ψ = 1 − 2m/r.Consequently, the Einstein field equation, R µν − g µν R/2 = 8πT µν (we choose units with G = c = 1), provides the following stellar structure equations: where the prime denotes a derivative with respect to r.
This set of equations ( 10)- (12) with the EoS (1) represents the internal structure of QSs.To solve them numerically we have to choose a specific value for ρ c (central energy density) with the following natural boundary conditions: and then integrate towards the surface of the star where the pressure drops to zero.This is identified as the stellar radius at P r (R) = 0. Additionally, boundary conditions are required to match this interior solution to an external metric, which is the Schwarzschild exterior solution.This is defined by with M = m(r = R) being the total mass of the star.Consequently, the surface gravitational redshift can be written in terms of radius and mass as follows

B. Moment of inertia
Under the slowly rotating approximation (i.e., we consider only first-order terms in the angular velocity of the star Ω) [41], the moment of inertia of an anisotropic compact star is given by [35] where ϖ = Ω − ω, with ω(r, θ) being the angular velocity acquired by an observer falling freely from infinity.In fact, it is possible to show that ϖ is a function only of the radial coordinate and is determined after solving the following differential equation (see Ref. [35] for further details) Suitable boundary conditions for the above differential equation come from the requirements of regularity at the origin and asymptotic flatness at very far distances from the star.Note that, at infinity (where spacetime is flat) the dragging angular velocity ω(r) goes to zero.Therefore, for an arbitrary choice of ϖ(0), Eq. ( 17) has to satisfy

C. Tidal Love number
Neutron stars are tidally deformed under the influence of an external tidal field (created by a companion star) [42].This deformation, quantified through a parameter termed the tidal deformability Λ, can be measured in gravitational waves emitted from the inspiral of a binary neutron-star coalescence.In that regard, here we calculate the tidal deformability of individual anisotropic QSs.The dimensionless tidal deformability is determined from the compactness C = M/R and the tidal Love number k 2 via [14,43,44] with k 2 being given by the following expression where α = y(R) − 4πR 3 ρ s /M.Here, ρ s stands for the energy density at the surface of the interacting QS.The correction term "−4πR 3 ρ s /M" is included because the energy density at the surface is finite and non-null [14,44].Actually, this becomes evident from the EoS (1).The function y(r) is obtained by solving the differential equation: where we have defined [14] with A = dP ⊥ /dP r and v sr = dP r /dρ being the radial speed of sound.Besides, the boundary condition for Eq. ( 21) at the stellar origin is y(0) = 2 [45].Therefore, y(R) is calculated after integrating equation ( 21) from the origin up to the surface of the anisotropic compact star.With that we determine α, and consequently the tidal Love number (20).

D. Radial stability
The TOV equations provide equilibrium solutions, however, such an equilibrium state can be stable or unstable [46].To examine the stellar stability of the equilibrium anisotropic configurations with interacting quark matter, it becomes necessary to determine the oscillation frequencies when a star is radially and adiabatically perturbed.The radial vibrations in the interior of an anisotropic compact star are governed by [31] where the perturbation ζ is related to the Lagrangian displacement ξ via ζ = ξ/r.This system of first-order differential equations has been obtained assuming that the perturbations have a harmonic time dependence of the form ∼ e iνt , where ν is the oscillation frequency to be determined.Moreover, ∆P r is the Lagrangian perturbation of the radial pressure, δσ is the Eulerian perturbation for the anisotropy factor, and γ = (1 + ρ/P r )dP r /dρ is the adiabatic index.Notice that, for any metric or fluid variable f , the relation between the Eulerian and Lagrangian perturbations is given by The expression for δσ will depend on the specific form of the anisotropic profile.At the surface of the star, where the pressure is zero, the Lagrangian perturbation of the radial pressure must satisfy ∆P r = 0.Meanwhile, the following boundary condition must be imposed at the stellar center

IV. OVERVIEW OF DIFFERENT ANISOTROPY MODELS
To generate numerical solutions in the present work, we will adopt two different anisotropy models, which have been widely used in the literature to model anisotropic matter at high densities and pressures: 1.The well-known Quasi-Local (QL) model as proposed by Horvat et al. [30], and define by where the free parameter β H measures the degree of anisotropy and the function µ = 2m(r)/r define the compactness.We focus in the range of −2 ≤ β H ≤ 2, and has been deeply discussed in [31,37,38,[47][48][49][50][51].
The Eulerian perturbation of the anisotropy factor takes the form [31] where 2. We are also interested in the anisotropy factor suggested by Bowers and Liang [20] where the free parameter β BL is an important quantity in describing the measure of anisotropy.Note that when β BL → 0, the anisotropy pressure vanishes at the center and this guarantees the regularity of Eq. ( 11).The meaningful value of β BL lies in the range of −2 ≤ β BL ≤ 2, see Refs.[26,37].
For this model, the Eulerian perturbation of the anisotropy factor is given by [31] where To describe realistic astrophysical objects, certain acceptability conditions have to be established or obeyed at any spacetime point of the star [52]: (a) Energy density, radial pressure and tangential pressure should be positive, (b) gradients for energy density and radial pressure must be negative, (c) the radial and tangential speed of sound should be less than the speed of light, (d) the interior solution must satisfy the energy conditions, (e) isotropy at the stellar center, among others.As we will see in our results, the models adopted in this work satisfy such physical acceptability conditions.
Using the Bowers-Liang model (30), Biswas and Bose showed that the anisotropic pressure in NSs can reduce their tidal deformability [26].Additionally, they claimed that certain equations of state that are ruled out by GW170817 observations (using isotropic matter) can become viable if the stars had a significant enough anisotropic pressure component.Rahmansyah and Sulaksono [51] have argued that the above two models are consistent with the recent NS multimessenger observations.Therefore, based on these interesting findings, there is a need to explore the astrophysical effect of each model on the internal structure of IQSs.

V. NUMERICAL RESULTS AND DISCUSSION
Given an anisotropy profile for σ and an EoS, we begin our analysis by solving the set of stellar structure equations ( 10)- (12).For our numerical calculations we will use B eff = 60 MeV/fm 3 and λ = 10 in the equation of state (1).Figure 3 illustrates the numerical solutions of the TOV equations for the anisotropy profile proposed by Horvat et al. [30] by considering a central density ρ c = 0.5 × 10 18 kg/m 3 and anisotropy parameter varying in the range β H ∈ [−0.6, 0.6].One observes that positive (negative) values of β H increase (decrease) the mass and radius of the interacting QS.Furthermore, the anisotropy factor σ vanishes at the origin (because of µ = 0 at r = 0) and at the surface (since P r = 0 at r = R), while the largest effect takes place in the intermediate zones.As is evident, the isotropic solutions are recovered when β H = 0.In Fig. 4 we present our numerical solutions for the anisotropy profile (30) with β BL ∈ [−0.3, 0.3].The radial behavior for the energy density and mass function is similar to the first model, however, the anisotropy factor σ is no longer zero at the surface.This last characteristic is due to the fact that the energy density does not disappear at the surface of a QS even if the radial pressure is zero, see EoS (1).Explicit values for the mass and radius of some stellar configurations represented in Figs. 3 and 4 are shown in Table I.It is worth mentioning that positive (negative) values of β H and β BL increase (decrease) the energy density and, according to Eq. (10), this should lead to increasing (decreasing) mass at any radial coordinate r within the stellar fluid.Consequently, all global properties of an interacting quark star will be altered due to the change in energy density and mass function.
The variation of the central density ρ c leads us to obtain a family of QSs in hydrostatic equilibrium, which can be presented in the so-called mass-radius diagram.In Fig. 5, we show such a diagram for a set of values of β H (see left plot) and β BL (see right plot).For both anisotropy models, we see that the main consequence of a positive anisotropy is a relevant increase in the maximum-mass values, while the opposite occurs for a negative anisotropy.The mass versus central density relations are displayed in Fig. 6, where it can be observed that the critical central density (corresponding to the maximum mass) is found to be at a lower value as the anisotropy increases.This means that, according to the standard criterion for stability dM/dρ c > 0, a positive (negative) anisotropy decreases (increases) the stability of an interacting QS.Notwithstanding, this classical criterion is a necessary but not sufficient condition for stellar stability.A sufficient condition involves determining the oscillation frequencies when a compact star is perturbed.We will return to this later.Once the mass and radius are known, we can then determine the gravitational redshift via expression (15).According to Fig. 7, z sur is affected by anisotropy mainly in the highmass branch, while the changes are irrelevant at small masses.This behavior was already expected due to the mass-radius relations shown in Fig. 5.
Before determining the moment of inertia of the anisotropic configurations presented in Fig. 5, we will analyze the radial behavior of the dragging angular velocity of a particle falling freely from infinity towards a slowly rotating star.Given a value of central density ρ c = 0.5 × 10 18 kg/m 3 , Fig. 8 shows the numerical solution of Eq. ( 17) inside and outside the star for both anisotropy models, where a positive anisotropy (with positive values of β H and β BL ) leads to increase ω, while the opposite occurs for negative anisotropies.The black curves correspond to isotropic solutions, and the dragging effects on the particle are smaller as the radial distance is sufficiently large, i.e. ω → 0 at infinity.Once we know the solution ϖ(r), we can find the moment of inertia through (16).According to Table I, given a central density value, I increases (decreases) as the anisotropies become more positive (negative).By repeating this procedure for a set of central densities, we can plot the moment of inertia as a function of mass, as presented in Fig. 9.One more time, we note that the effects of anisotropy become stronger only in the high-mass re-gion.
Our next step is to investigate the dependence of the tidal Love number on the anisotropic pressure.Given an anisotropy profile and a central density value, this involves solving Eq. ( 21), determining α and therefore k 2 via Eq.(20).Given a value of ρ c , the numerical results shown in Table I reveal that k 2 decreases with the increase of β H and β BL .In Fig. 10 we present the tidal Love number as a function of compactness for both anisotropy models.For the quasi-local ansatz, the impact of anisotropy on k 2 is negligible.Specifically, in the small-compactness region, k 2 decreases for more positive values of β H , but this behavior is reversed for large compactnesses (see inset on the left plot).Nonetheless, the effect of anisotropic pressure on k 2 is larger for the Bowers-Liang model, and the tidal Love number always increases with increasing β BL .Furthermore, the dimensionless tidal deformability (19) as a function of gravitational mass is plotted in Fig. 11.Once again we can observe that the most relevant effects on Λ due to anisotropy occur in the high-mass region.
Finally, to know whether the anisotropic configurations constructed in this work are dynamically stable, it is necessary to analyze the radial vibration modes.We begin our analysis by solving the differential equations (24) and (25) for a value of central density ρ c = 0.5 × 10 18 kg/m 3 , see Fig. 12.The upper (lower) plot of this figure corresponds to the first (second) anisotropy model for β H = 0.4 (β BL = 0.2), where we have shown the eigenfunctions ∆p r,n for the lowest six normal oscillation modes.For this stellar configuration, each nth vibration mode is associated with an eigenvalue or squared frequency ν 2 n so that we obtain a frequency spectrum of the form ν 2 0 < ν 2 1 < ν 2 2 < • • • , where n = 0 indicates the fundamental mode and it has no nodes, n = 1 is the first overtone (has one node), and so forth.Purely real frequencies describe dynamically stable stars (namely, when ν 2 n > 0), while imaginary frequencies describe unstable stars.For the sequence of equilibrium configurations presented in Fig. 5 with five values of β H , we have calculated the squared frequencies of the fundamental pulsation mode ν 2 0 for the quasi-local profile (27) and are shown in Fig. 13 as a function of the central density (left panel) and total mass (right panel).According to the right plot (see inset), the maximum mass corresponds precisely to ν 2 0 = 0, which means that the maximum-mass point can be used as an indicator of the transition from stability to instability of interacting quark stars with anisotropies described by the quasi-local profile.However, the right plot of Fig. 14 shows that the maximum-mass points do not coincide with ν 2 0 = 0 for β BL ̸ = 0. Therefore, the usual sta- The increase in either of the parameters β H or β BL implies that the critical configuration (corresponding to the maximum-mass point) is found at a smaller central-density value.For the horizontal axis in both plots, we have used the common logarithm (i.e., the logarithm with base 10).The same applies to Figs. 13 and 14.
bility criterion dM/dρ c > 0 is no longer valid for the Bowers-Liang model.This result is similar to that obtained for anisotropic neutron stars [31].Indeed, negative anisotropies increase the radial stability of interacting quark stars in the sense that it is possible to obtain stable stars slightly beyond the maximum mass for β BL < 0. Note further that, for both anisotropy models, the squared frequencies of the fundamental mode decrease with increasing central density.

VI. THEORETICAL PREDICTION AND OBSERVATIONAL DATA
Astrophysical measurements, based on observations of NSs in the X-ray band and on gravitational wave observations, play a crucial role in adopting or discarding an EoS.In that perspective, it is essential to compare our theoretical results with the recent observational measurements/constraints.In Fig. 15, we compare our mass-radius results with the NICER measurements ob-  17) for a central density ρ c = 0.5 × 10 18 kg/m 3 and EoS (1) for B eff = 60 MeV/fm 3 and λ = 10.We have also considered five values for β H and β BL .The solid lines represent the interior solutions, while dashed lines indicate the exterior region.A positive anisotropy leads to significantly increasing the angular velocity mainly in the innermost regions of the interacting quark star.As expected, the drag experienced (in the direction of rotation) by a particle falling toward the star is increasingly greater as it approaches the surface, so that ω → 0 for r → ∞.
tained from the pulsars PSR J0030+0451 [4,59] and PSR J0740+6620 [60,61].We have included the constraint from the GW170817 event [58], as well as the mass of the secondary component of the gravitational-wave signal GW190814 [56].We also incorporate the central compact object within the supernova remnant HESS J1731−347, whose mass and radius estimates have been reported in Ref. [57].All these measurements are being represented by different colored bands or contours in both plots of Fig. 15.
We begin our comparative analysis with the noninteracting case (i.e., when λ = 0) for B eff = 60 MeV/fm 3 , see black curves in the left panel of Fig. 15.As one can observe, isotropic solutions (β H = 0) do not provide astrophysical masses above 2 M ⊙ , which does not favor the description of highly massive pulsars.However, for positive anisotropies it is possible to exceed two solar masses, see the exemplary cases when β H = 0.3 and 0.6.Remarkably, when strong interaction effects take place (i.e., when λ > 0), we obtain results that satisfy most observational M − R measurements.In particular, for λ ∈ [0.1, 2.0], our outcomes suggest that the supernova remnant HESS J1731−347 can be described as an anisotropic IQS, regardless of the value of β H .This is allowed because the effect of anisotropy is irrelevant for small masses.Nonetheless, the anisotropic pressure plays an important role in the high-mass region by satisfying the constraints obtained from highly massive Pulsars observations, see specifically the green and orange curves.We further note that λ = 0.5 and β H ∈ [−0.4,0.6] favor a consistent description for the secondary component resulting from the GW190814 event.Besides, in the right plot of Fig. 15, we display the M − R outcomes corresponding to B eff = 90 MeV/fm 3 , where we can perceive that an increase in the effective bag constant leads to a decrease in the maximum mass, but it is still possible to obtain results compatible with the observational data, see for example the red curves for λ = 2.0.
For the anisotropic IQSs presented in the left plot of Fig. 15, the dimensionless tidal deformability as a function of the total mass is shown in Fig. 16, where we have included the estimate reported by LIGO-Virgo Collaboration for the tidal effects of the GW170817 event [58], namely Λ 1.4M ⊙ = 190 +390 −120 .It can be seen that, for λ ≲ 0.018, any value of the anisotropy parameter β H satisfies such an estimate.Nevertheless, we must point out that, this constraint is based on the assumption that the coalescing bodies were NSs.Actually, for quark stars in the present study, we must use the EoSindependent constraint Λ 1.4M ⊙ ≤ 800 (see the dark yel-  (27).According to the left plot, the central density value corresponding to ν 2 0 = 0 is found to be an increasing value as β H decreases.This means that negative anisotropies increase the radial stability of IQSs.Furthermore, according to the right plot, the maximum mass is reached exactly when ν 2 0 = 0, regardless of the value of β H . Notice that filled circles in the inner box mark the maximum-mass FIG.14. Squared oscillation frequencies as in Fig. 13, but now for the anisotropy model proposed by Bowers and Liang (30).As in the first model, ρ c corresponding to ν 2 0 = 0 increases as β BL becomes more negative.However, according to the right plot, the maximum mass is no longer found at ν 2 0 = 0 for β BL ̸ = 0.The existence of stable stars beyond the maximum mass for negative anisotropies is possible.low dashed vertical line) from the older LIGO/Virgo paper [62], which is satisfied for λ ≲ 0.1 if we adopt B eff = 60 MeV/fm 3 .As one can appreciate, this last restriction can further open the free parameter space for IQSs.Noticeably, when we assume B eff = 90 MeV/fm 3 , all our findings are consistent with the tidal deformability constraints for λ ∈ [0, 2.0], as shown in Fig. 17.Thus, quark stars made of interacting quark matter with large effective bag constant and proper λ choices can meet the tidal deformability constraint of the GW170817  [53] and J0348+0432 [54].The blue horizontal band indicates the observational mass measurement of the massive NS pulsar J2215+5135 [55], while the filled green stripe stands for the lower mass of the compact object detected by the GW 190814 event [56].The purple dot with its respective error bars represents the supernova remnant HESS J1731−347 [57] and the cyan region is the mass-radius constraint from the GW170817 event [58].Moreover, NICER constraints obtained from the pulsars PSR J0030+0451 [4,59] and PSR J0740+6620 [60,61] are indicated by different color contours.Note that, as a consequence of increasing B eff , the radius and mass of the maximum-mass configuration decrease.
event while satisfying the maximum mass restrictions, see also right panel of Fig. 15.

VII. CONCLUSIONS
Within the Standard Model of particle physics, Quantum chromodynamics (QCD) is the accepted theory of strong interaction force between quarks mediated by gluons.At the same time, perturbative expansions in QCD and color superconductivity allow us to describe interacting quark matter that includes the interquark effects induced by strong interaction.This may reveal new physical phenomena under extreme conditions such as inside compact stars.In this work, we have investigated the role of strong interaction on anisotropic QSs composed of interacting quark matter.Our study involves spherically symmetric anisotropic QSs adopting two physically well-motivated anisotropy profiles.To achieve this end, we have numerically solved the stellar structure equations in Einstein gravity using suitable boundary conditions and we discussed various global properties such as mass-radius relation, surface redshift, moment of inertia, tidal properties and oscillation spectrum.
The local anisotropy has been included with the help of two free parameter β H and β BL for the adopted models.By choosing particular values of these parameters, we have obtained the M − R relations for anisotropic interacting quark stars, as illustrated in Fig. 5.We noticed that the mass and radius simultaneously increase for positive increasing values of anisotropy factor, and this holds for both models, see Table I.According to the M(ρ c ) method in the M − ρ c diagram in Fig. 6, where the maximum-mass point indicates that stability ceases, the critical central density (corresponding to the maximum mass) decreases with the increase of the parameters β H and β BL .In other words, a positive (negative) anisotropy decreases (increases) the stability of IQSs.Depending on the obtained values of mass and radius, we also determined the gravitational redshift z sur , which is strongly affected by the presence of anisotropy −120 reported by LIGO-Virgo Collaboration by analyzing the tidal effects of the GW170817 event [58].The dark yellow dashed line stands for the EoS-independent constraint Λ 1.4M ⊙ ≤ 800 from Ref. [62].Here, all our results are consistent with the EoS-independent constraint Λ 1.4M ⊙ ≤ 800 [62].
in the high-mass region.
Our next step was to analyze the behavior of the frame-dragging angular velocity and moment of inertia in the slowly rotating approximation.In particular, given a central density value, we have shown that a positive anisotropy leads to significantly increasing the angular velocity mainly in the innermost regions of the IQSs, while the opposite occurs for negative anisotropies.The anisotropic pressure consequently increases the moment of inertia of the stars analyzed as β H and β BL increase from their negative values, showing its greatest impact on the most massive stars.
In our analysis of tidal properties we found that the tidal Love number k 2 is weakly altered by anisotropy for the model proposed by Horvat et al., while the changes are more pronounced for the Bowers-Liang model.Specifically, positive (negative) values of the parameter β BL increase (decrease) the tidal Love number of IQSs.Moreover, the main effect of anisotropy on the dimensionless tidal deformability Λ takes place in the high-mass branch.
The already mentioned M(ρ c ) method provides a simple and necessary condition for stellar stability, however, a sufficient condition is to determine the frequencies of normal pulsation modes when a star is perturbed radially and adiabatically.So we have also performed a rigorous examination of the anisotropic stellar configurations constructed here following a radial oscillations approach for each anisotropy profile.Noticeably, for both anisotropy models, the central density value corresponding to zero squared frequency of the fundamental vibration mode (namely, ν 2 0 = 0) was found to be an increasing value as β H and β BL decrease.In other words, negative anisotropies increase the radial stability of IQSs, as predicted by the M(ρ c ) method.For the quasi-local model, the M(ρ c ) method is compatible with the calculation of radial oscillation frequencies, so the maximum mass can be used as an indicator of the onset of instability for any value of β H . Notwithstanding, for the Bowers-Liang profile, the central density corresponding to the maximum-mass point does not coincide with the central density where the squared oscillation frequency vanishes, indicating that the existence of stable anisotropic IQSs is possible beyond the maximum mass for negative anisotropies.Meanwhile, stars stop being stable before reaching the maximum mass for positive values of β BL .This is an important outcome that substantially differentiates one anisotropy model from the other.
Finally, we have carried out a comparative analysis between our theoretical results and observational data.It is important to highlight that, when strong interac-tion effects take place (i.e., for λ > 0), we obtained results that satisfy most observational mass-radius measurements.Specifically, when λ ∈ [0.1, 2.0], our calculations suggest that the central compact object within the supernova remnant HESS J1731−347 can be described as an anisotropic interacting quark star, regardless of the value of β H . Furthermore, according to the outcomes shown in Figs.5-11, we observe that the anisotropy profiles yield identical quark-star characteristics, at least for the range of values of β H and β BL considered in this study.Thus, with the exception of the radial stability analysis, we must remark that the effect of anisotropy generated by both models on IQSs is similar.However, as shown in Refs.[31,51] from the mass-radius diagrams of anisotropic NSs, there is a notable difference between the findings of the two anisotropy profiles not only for high masses but also in the low-mass region.Here, our work has shown that this difference is more subtle for interacting quark matter, where the M − R relation is almost unalterable due to anisotropy in the smaller mass branch.Remarkably, our findings are comparable and consistent with the tidal deformability constraint from the GW170817 event [62] while satisfying the different observational mass-radius measurements.

FIG. 1 .
FIG.1.Equation of state involving interacting quark matter in its dimensionless form(5) for several values of λ, where the linear behavior for λ → ∞ is described by the blue curve, see Eq. (6).Moreover, we recover the standard noninteracting quark matter EoS when λ = 0, which is represented by the black curve.

6 FIG. 3 .FIG. 4 .
FIG. 3. Radial profiles of the energy density (left panel), mass function (middle panel) and anisotropy factor (right panel) inside an anisotropic interacting quark star with central density ρ c = 0.5 × 10 18 kg/m 3 and EoS (1) for B eff = 60 MeV/fm 3 and λ = 10.Here we have used the quasi-local model (27) for values of β H in the range β H ∈ [−0.6, 0.6], see the color bar on the right side of the third plot.The main consequence of positive (negative) anisotropy is a relevant increase (decrease) in the radius and mass of the star.Remark that the anisotropy is more pronounced in the intermediate regions, and the particular case β H = 0 corresponds to zero anisotropic pressure.

3 FIG. 7 .FIG. 8 .
FIG.7.Surface gravitational redshift versus total mass for both anisotropy profiles.In both plots, the impact of anisotropy is substantial only in the high-mass region.Larger redshifts are achieved for positive anisotropies.

3 FIG. 9 . 3 FIG. 10 .
FIG.9.Moment of inertia as a function of the mass for the anisotropic interacting quark stars presented in Fig.5, where an increase in the values of β H and β BL implies larger moments of inertia.Additionally note that, regardless of the anisotropy model, the changes in the moment of inertia due to anisotropic pressure are negligible in the small-mass region.

3 FIG. 11 .FIG. 12 . 5 FIG. 13 .
FIG.11.Dimensionless tidal deformability(19) as a function of the gravitational mass for both anisotropy profiles.One can see that the greatest impact of anisotropy on tidal deformability takes place in the high-mass branch.

FIG. 15 .
FIG.15.Comparison of mass-radius theoretical results for anisotropic interacting quark stars with observational data using the quasi-local model(27) and four values of λ = {0, 0.1, 0.5, 2.0} for B eff = 60 MeV/fm 3 (left panel) and B eff = 90 MeV/fm 3 (right panel).For each value of λ we have considered five values of β H , where the isotropic case corresponds to β H = 0.The gray line at 2.0 M ⊙ represents the two massive NS pulsars J1614−2230[53] and J0348+0432[54].The blue horizontal band indicates the observational mass measurement of the massive NS pulsar J2215+5135[55], while the filled green stripe stands for the lower mass of the compact object detected by the GW 190814 event[56].The purple dot with its respective error bars represents the supernova remnant HESS J1731−347[57] and the cyan region is the mass-radius constraint from the GW170817 event[58].Moreover, NICER constraints obtained from the pulsars PSR J0030+0451[4,59] and PSR J0740+6620[60,61] are indicated by different color contours.Note that, as a consequence of increasing B eff , the radius and mass of the maximum-mass configuration decrease.

TABLE I .
Global properties of interacting QSs with central density ρ c = 0.5 × 10 18 kg/m 3 and EoS (1) for B eff = 60 MeV/fm 3 and λ = 10.The radial behavior of the energy density, mass and anisotropy factor of these stars is shown in Figs.3 and 4for the quasi-local ansatz and Bowers-Liang profile, respectively.Note that the frequency of the fundamental mode and first overtone are given by f 0 = ν 0 /2π and f 1 = ν 1 /2π, respectively.
16G.16.Dimensionless tidal deformability versus total mass relation for the stellar configurations presented in the left panel of Fig.15, where we have additionally included the results for λ = 0.018 by gray curves.The blue vertical line denotes Λ 1.4M ⊙ = 190 +390