Radial pulsations, moment of inertia and tidal deformability of dark energy stars

We construct dark energy stars with Chaplygin-type equation of state (EoS) in the presence of anisotropic pressure within the framework of Einstein gravity. From the classification established by Iyer et al. (Class Quantum Grav 2:219, 1985), we discuss the possible existence of isotropic dark energy stars as compact objects. However, there is the possibility of constructing ultra-compact stars for sufficiently large anisotropies. We investigate the stellar stability against radial oscillations, and we also determine the moment of inertia and tidal deformability of these stars. We find that the usual static criterion for radial stability dM/dρc>0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$dM/d\rho _c >0$$\end{document} still holds for dark energy stars since the squared frequency of the fundamental pulsation mode vanishes at the critical central density corresponding to the maximum-mass configuration. The dependence of the tidal Love number on the anisotropy parameter α\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha $$\end{document} is also examined. We show that the surface gravitational redshift, moment of inertia and dimensionless tidal deformability undergo significant changes due to anisotropic pressure, primarily in the high-mass region. Furthermore, in light of the detection of gravitational waves GW190814, we explore the possibility of describing the secondary component of such event as a stable dark energy star in the presence of anisotropy.


I. INTRODUCTION
Different types of observations (such as Type Ia supernovae, structure formation and CMB anisotropies) indicate that our Universe is not only expanding, it is accelerating.Within the standard ΛCDM model (which is based on cold dark matter and cosmological constant in Einstein gravity), this cosmic acceleration is due to a smooth component with large negative pressure and repulsive gravity, the so-called dark energy.Such a model gives a good agreement with the recent observational data [1], but suffers from the well-known coincidence problem and the fine-tuning problem [2,3].The exact physical nature of dark energy is still a mystery and, consequently, the possibility that dark matter and dark energy could be different manifestations of a single substance has been considered [4][5][6][7].In that regard, it was shown that the inhomogeneous Chaplygin gas offers a simple unified model of dark matter and dark energy [8].It was also argued that if the Universe is dominated by the Chaplygin gas a cosmological constant would be ruled out with high confidence [9].
Using the Planck 2015 CMB anisotropy, type-Ia supernovae and observed Hubble parameter data sets, the full parameter space of the modified Chaplygin gas was measured by Li et al. [10].Based on recent observations of high-redshift quasars, Zheng and colleagues [11] investigated a series of Chaplygin gas models as candidates for dark matter-energy unification.The application of the Hamilton-Jacobi formalism for generalized Chaplygin gas models was carried out in Ref. [12].Additionally, it is worth mentioning that Odintsov et al. [13] considered two different equations of state for dark energy (i.e., powerlaw and logarithmic effective corrections to the pressure).* juanzarate@cbpf.brThey showed that the power-law model only yielded some modest results, achieved under negative values of bulk viscosity, while the logarithmic scenario provide good fits in comparison to the ΛCDM model.
Another way to give rise to an accelerated expansion of the Universe is by modifying the geometry itself [14,15], namely, considering higher curvature corrections to the standard Einstein-Hilbert action.Under this outlook, the cosmic acceleration can be modeled in the scope of a scalar-tensor gravity theory [16,17].Moreover, within the context of the so-called f (R) theories [18,19], the quadratic term in the Ricci scalar R leads to an inflationary solution in the early Universe [20], although such a model does not provide a late-time accelerated expansion.Nevertheless, the late-time acceleration era can be realized by terms containing inverse powers of R [21], though it was shown that this is not compatible with the solar system experiments [22].For a comprehensive study on the evolution of the early and present Universe in f (R) modified gravity, we refer the reader to the review articles [23][24][25] and references contained therein.On the other hand, the astrophysical implications due to the f (R) modified gravitational Lagrangian on compact stars have been intensively investigated in the past few years [26][27][28][29][30][31][32][33].
According to the aforementioned works, different dark energy models have been proposed in order to explain the mechanisms that lead to the cosmic acceleration.Only about 4% of the Universe is made of familiar atomic matter, 20% dark matter, and it turns out that roughly 76% of the Universe is dark energy [34].Within the context of General Relativity, dark energy is an exotic negative pressure contribution that can lead to the observed accelerated expansion.In the absence of consensus regarding a theoretical description for the current accelerated expansion of the Universe, theorists have proposed using the Chaplygin gas as a useful phenomenological description [4].If dark energy is distributed anywhere permeating ordinary matter, then it could be present in the interior of a compact star.Therefore, the purpose of this manuscript is to investigate the possible existence of compact stars with dark energy by assuming a Chaplygin-type EoS.For such stars to exist in nature, they need to be stable under small radial perturbations.
Adopting a description of dark energy by means of a phantom (ghost) scalar field, Yazadjiev [35] constructed a general class of exact interior solutions describing mixed relativistic stars containing both ordinary matter and dark energy.The energy conditions and gravitational wave echoes of such stars were recently analyzed in Ref. [36].Furthermore, the effect of the dynamical scalar field quintessence dark energy on neutron stars was investigated in [37].Panotopoulos and collaborators [38] studied slowly rotating dark energy stars made of isotropic matter using the Chaplygin EoS.Bhar [39] proposed a model for a dark energy star made of dark and ordinary matter in the Tolman-Kuchowicz spacetime geometry.For further stellar models with dark energy we also refer the reader to Refs.[40][41][42][43][44][45][46][47][48].
In addition, anisotropy in compact stars may arise due to strong magnetic fields, pion condensation, phase transitions, mixture of two fluids, bosonic composition, rotation, etc.Thus, regardless of the specific source of the anisotropy, it is more natural to think of anisotropic fluids when studying compact stars at densities above nuclear saturation density.In that regard, the literature offers some physically motivated functional relations for the anisotropy, see for example Refs.[49][50][51][52][53][54][55].However, we must point out that these anisotropic models are based on general assumptions (or ansatzes) that do not directly relate to exotic modifications of matter or gravity.Indeed, it has been argued that the deformation near the maximum neutron-star mass comes from the anisotropic pressure within these stars, which is caused by the distortion of Fermi surface predicted by the equation of state of the models [56].Becerra-Vergara et al. [57] showed that the contribution of the fourth order corrections parameter (a 4 ) of the QCD perturbation on the radial and tangential pressure generate significant effects on the mass-radius relation and the stability of quark stars.It has also been shown that the stellar structure equations in Eddington-inspired Born-Infeld theory with isotropic matter can be recast into GR with a modified (apparent) anisotropic matter [58].
Motivated by the several works already mentioned, we aim to discuss the impact of anisotropy on the macroscopic properties of dark energy stars with Chaplygin-like EoS.We will address the following questions: Do these stars belong to families of compact or ultra-compact stars?How does anisotropy affect the compactness and radial stability of dark energy stars satisfying the causality condition?In particular, by adopting the phenomenological ansatz proposed by Horvat et al. [51], we determine the radius, mass, gravitational redshift, frequency of the fundamental oscillation mode, moment of inertia and the dimensionless tidal deformability of anisotropic dark energy stars.The isotropic solutions are recovered when the anisotropy parameter vanishes, i.e. when α = 0.
The organization of this paper is as follows: In Sec.II we start with a brief overview of relativistic stellar structure, describing the basic equations for radial pulsations, moment of inertia and tidal deformability.We then introduce the Chaplygin-like EoS and discuss its relation to the cosmological context in Sec.III, as well as we present the anisotropy profile.Section IV provides a discussion of the numerical results for the different physical properties of dark energy stars.Finally, our conclusions are summarized in Sec.V.

II. STELLAR STRUCTURE EQUATIONS
In order to study the basic features of compact stars with dark energy, in this section we briefly summarize the stellar structure equations in Einstein gravity.In particular, we focus on hydrostatic equilibrium structure, radial pulsations, moment of inertia, and tidal deformability.
The theory of gravity to be used in this work is general relativity, where the Einstein field equations are given by with G µν being the Einstein tensor, R µν the Ricci tensor, R denotes the scalar curvature, and T µν is the energymomentum tensor.Since we are interested in isolated compact stars, we consider that the spacetime can be described by the spherically symmetric four-dimensional line element In addition, we model the compact-star matter by an anisotropic perfect fluid, whose energy-momentum tensor is given by where ρ is the energy density, σ ≡ p t − p r the anisotropy factor, p r the radial pressure, p t the tangential pressure, u µ the four-velocity of the fluid, and k µ is a unit fourvector.These four-vectors must satisfy u µ u µ = −1, k µ k µ = 1 and u µ k µ = 0. Notice that the stellar fluid becomes isotropic when σ = 0.

A. TOV equations
When the stellar fluid remains in hydrostatic equilibrium, neither metric nor thermodynamic quantities depend on the time coordinate.This allows us to write u µ = e −ψ δ µ 0 and k µ = e −λ δ µ 1 .Accordingly, the hydrostatic equilibrium of an anisotropic compact star is gov-erned by the TOV equations: which are obtained from Eqs. ( 1)-( 3) together with the conservation law ∇ µ T µ 1 = 0.The metric function λ(r) is determined from the relation e −2λ = 1 − 2m/r, where m(r) is the gravitational mass within a sphere of radius r.
By supplying an EoS for the radial pressure in the form p r = p r (ρ) and a defined anisotropy relation for σ, the system of differential equations ( 4)-( 6) is then numerically integrated from the center at r = 0 to the surface of the star r = R which correspond to a vanishing pressure.Therefore, the above equations will be solved under the requirement of the following boundary conditions where ρ c is the central energy density, and M ≡ m(R) is the total mass of the star calculated at its surface.The numerical solution of the TOV equations describes the equilibrium background and allow us to obtain the metric components and fluid variables.

B. Radial oscillations
A rigorous analysis of the radial stability of compact stars requires the calculation of the frequencies of normal vibration modes.Such frequencies can be found by considering small deviations from the hydrostatic equilibrium state but maintaining the spherical symmetry of the star.In the linear treatment, where all quadratic (or higher-order) or mixed terms in the perturbations are discarded, one assumes that all perturbations in physical quantities are arbitrarily small.The fluid element located at r in the unperturbed configuration is displaced to radial coordinate r + ξ(t, r) in the perturbed configuration, where ξ is the Lagrangian displacement.All perturbations have a harmonic time dependence of the form ∼ e iνt , where ν is the oscillation frequency to be determined.Consequently, defining ζ ≡ ξ/r, the adiabatic 1 radial pulsations of anisotropic compact stars are 1 In the adiabatic theory, it is assumed that the fluid elements of the star neither gain nor lose heat during the oscillation.
governed by the following differential equations [55] where ∆p r is the Lagrangian perturbation of the radial pressure and γ = (1 + ρ/p r )dp r /dρ is the adiabatic index at constant specific entropy.The above first-order time-independent equations ( 8) and (9) require boundary conditions set at the center and surface of the star, similar to a vibrating string fixed at its ends.Since Eq. ( 8) has a singularity at the origin, the following condition must be required while the Lagrangian perturbation of the radial pressure at the surface must satisfy C. Moment of inertia Suppose a particle is dropped from rest at a great distance from a rotating star, then it would experience an ever increasing drag in the direction of rotation as it approaches the star.Based on this description, we introduce the angular velocity acquired by an observer falling freely from infinity, denoted by ω(r, θ).Here we will calculate the moment of inertia of an anisotropic dark energy star under the slowly rotating approximation [59].This means that when we consider rotational corrections only to first order in the angular velocity of the star Ω, the line element (2) is replaced by its slowly rotating counterpart, namely and following Ref.[59], it is pertinent to define the difference ≡ Ω − ω as the coordinate angular velocity of the fluid element at (r, θ) seen by the freely falling observer.Keep in mind that Ω is the angular velocity of the stellar fluid as seen by an observer at rest at some spacetime point (t, r, θ, φ), and hence the four-velocity up to linear terms in Ω can be written as u µ = (e −ψ , 0, 0, Ωe −ψ ).To this order, the spherical symmetry is still preserved and it is possible to extend the validity of the TOV equations ( 4)- (6).Nonetheless, the 03-component of the field equations contributes an additional differential equation for angular velocity.By retaining only first-order terms in Ω, such component becomes As in the case of isotropic fluids, we follow the same treatment carried out by Hartle [59,60] and we assume that can be written as where P l are Legendre polynomials.Taking this expansion into account, Eq. ( 13) becomes At a distance far away from the star, where e −(ψ+λ) becomes unity, the asymptotic solution of Eq. ( 15) takes the form l (r) → a 1 r −l−2 + a 2 r l−1 .If spacetime is to be flat at large r, then ω → 2J/r 3 (or equivalently, → Ω − 2J/r 3 ) for r → ∞, where J is the total angular momentum of the star [59,61].Therefore, comparing this with the asymptotic behavior of l (r), we find that l = 1.As a result, is a function only of the radial coordinate, and Eq. ( 15) reduces to which can be integrated to give In view of Eq. ( 17), we can obtain the angular momentum J and hence the moment of inertia I = J/Ω of a slowly rotating anisotropic star: which reduces to the expression given in Ref. [61] for isotropic compact stars when σ = 0.For an arbitrary choice of the central value (0), the appropriate boundary conditions for the differential equation ( 16) come from the requirements of regularity at the center of the star and asymptotic flatness at infinity, namely Once the solution for (r) is found, we can then determine the moment of inertia through the integral (18).It is remarkable that the above expression for I is referred to as the "slowly rotating" approximation because it was obtained to lowest order in the angular velocity Ω [61].This means that the stellar structure equations are still given by the TOV equations ( 4)-( 6).

D. Tidal deformability
It is well known that the tidal properties of neutron stars are measurable in gravitational waves emitted from the inspiral of a binary neutron-star coalescence [62,63].In that regard, here we also study the dimensionless tidal deformability of individual dark energy stars.To do so, we follow the procedure carried out by Hinderer et al. [64] (see also Refs.[65][66][67][68][69][70] for additional results).The basic idea is as follows: In a binary system, the deformation of a compact star due to the tidal effect created by the companion star is characterized by the tidal deformability parameter λ = −Q ij /E ij , where Q ij is the induced quadrupole moment tensor and E ij is the tidal field tensor [68].Namely, the latter describes the tidal field from the spacetime curvature sourced by the distant companion.
The tidal parameter is related to the tidal Love number k 2 through the relation but it is common in the literature to define the dimensionless tidal deformability Λ = λ/M 5 , so in our results we will focus on Λ.The calculation of λ requires considering linear quadrupolar perturbations (due to the external tidal field) to the equilibrium configuration.Thus, the spacetime metric is given by g µν = g 0 µν + h µν , where g 0 µν describes the equilibrium configuration and h µν is a linearized metric perturbation.For static and even-parity perturbations in the Regge-Wheeler gauge [71], the perturbed metric can be written as [64] where H 0 , H 2 and K are functions of the radial coordinate, and Y lm are the spherical harmonics for l = 2. Since the perturbed energy-momentum tensor is given by δT ν µ = diag(−δρ, δp r , δp t , δp t ), the linearized field equations imply that: In addition, from δG 0 0 − δG 1 1 = −8π(δρ + δp t ), we can obtain the following differential equation [72] H + PH + QH = 0, (22) or alternatively, where we have defined with A ≡ dp t /dp r and v sr being the radial speed of sound.
By matching the internal solution with the external solution of the perturbed variable H at the surface of the star r = R, we obtain the tidal Love number [72] where C ≡ M/R is the compactness of the star, and y R ≡ y(R) is obtained by integrating equation ( 23) from the origin up to the stellar surface.

III. EQUATION OF STATE AND ANISOTROPY MODEL
As it is well known, a possible alternative to the Phantom and Quintessence fields is the Chaplygin gas, where the EoS assumes the form p r = −B/ρ, with B being a positive constant (given in m −4 units).In fact, it was argued that such gas could provide a solution to unify the effects of dark matter in the early times and dark energy in late times [4,11].Although the literature provides a more generalized version for such EoS in the context of the Friedmann-Lemaître-Robertson-Walker Universe [5][6][7][73][74][75][76][77], here we will use the simplest form plus a linear term corresponding to a barotropic fluid, namely where A is a positive dimensionless constant.Our model is characterized by two free parameters A and B. Nevertheless, we must emphasize here that Li et al. [10] considered an equation of state with three degrees of freedom, specifically p = Aρ − B/ρ α , where α is an extra parameter.They carried out a statistical treatment of astronomical data in order to constrain the parameter space.
In the light of the Markov chain Monte Carlo method, they found that at 2σ level, α = −0.0156+0.0982+0.2346−0.1380−0.2180and A = 0.0009 +0.0018+0.0030−0.0017−0.0030from CMB+JLA+CC data sets.In other words, the constants α and A are very close to zero and hence the nature of unified dark matterenergy model is very similar to the cosmological standard ΛCDM model.
On the other hand, at astrophysics level, compact stars obeying the EoS (28) have been investigated by several authors, see for example Refs.[38,41,[43][44][45].In this work we will adopt values of A and B for which appreciable changes in the mass-radius diagram can be visualized in order to compare our theoretical results with observational measurements of massive pulsars.
In order to describe physically realistic compact stars, the causality condition must be respected throughout the interior region of the star.In other words, the speed of sound (defined by v s ≡ dp/dρ) cannot be greater than the speed of light.Thus, in view of Eq. ( 28), we have and since the radial pressure vanishes at the surface of the star, then B = Aρ 2 .Thereby, the causality condition v 2 sr (R) = 2A < 1 implies that A < 0.5.Besides, it is more realistic to consider stellar models where there exists a tangential pressure as well as a radial one, since anisotropies arise at high densities, i.e. above the nuclear saturation density as considered in this work.Although the literature offers different functional relations to model anisotropic pressures at very high densities inside compact stars [49][50][51][52][53][54], here we adopt the simplest model, which was proposed by Horvat and collaborators [51] where α is a dimensionless parameter that controls the amount of anisotropy within the stellar fluid.This parameter can assume positive or negative values of the order of unity, see Refs.[26,32,51,52,55,[78][79][80][81][82].Notice that the isotropic solutions are recovered when the value of α vanishes.Specifically, the anisotropy ansatz (30) has two important characteristics: (i) the fluid becomes isotropic at the center generating regular solutions and (ii) the effect of anisotropy vanishes in the hydrostatic equilibrium equation in the Newtonian limit.Unlike this profile, the effect of anisotropy does not vanish in the hydrostatic equilibrium equation in the nonrelativistic regime for the Bowers-Liang model [49], which could be an unphysical trait as argued in Ref. [79].For a broader discussion on the different ways of generating static spherically symmetric anisotropic fluid solutions, we refer the reader to the recent review article [83].
Since the Eulerian perturbation for the metric potential λ can be written as δλ = −4πr(ρ + p r )e 2λ ξ [55], then δσ takes the form where it should be noted that the relation between the Eulerian and Lagrangian perturbations for radial pressure is given by ∆p r = δp r + rζp r .The above expression will be substituted in Eq. ( 9) when we discuss later the radial pulsations in the stellar interior for at least some values of α.

A. Equilibrium configurations
So far we do not know exactly whether the millisecond pulsars (observed in compact binaries from optical spectroscopic and photometric measurements) are hadronic, quark or hybrid stars.In fact, it has been theorized that cold quark matter might exist at the core of heavy neutron stars [84].Despite the precise measurements of masses [85][86][87] and radii [88][89][90], such constraints are still unable to distinguish the theoretical predictions coming from the different models for strange stars and (hybrid) neutron stars.This means that the dense matter EoS within compact stars still remains poorly understood.Furthermore, a realistic compact star possesses high magnetic fields and rotation properties, which significantly alter its internal structure.For comparison reasons, it is therefore common to use the observational mass-radius measurements (in view of the detection of gravitational waves and electromagnetic signals) on the mass-radius diagrams for any type of EoS even being of different microscopic compositions.In that perspective, our theoretical results will be compared with observational measurements.
We begin our discussion of dark energy stars by considering the isotropic case (i.e., when σ = 0 in the TOV equations).We numerically integrate Eqs. ( 4)-( 6) from the center up to the surface of the star through the boundary conditions (7).As usual, the radius R is determined when the pressure vanishes, and the total mass M is calculated at the surface.The felt panel of Fig. 1 exhibits the mass-radius relations of dark energy stars for different values of parameters A and B in the EoS (28).Remark that we have adopted values of A less than 0.5 in order to respect the causality condition.One can observe that small values of A (see black curve) do not provide compact stars that fit current observational data.However, higher values of maximum mass can be obtained for larger values of A, see for example red and green curves.
For a fixed value of A, the maximum mass decreases as the parameter B increases.We perceive that the secondary component resulting from the gravitational-wave signal GW190814 [91] can be consistently described as a compact star with Chaplygin EoS (28) for A = 0.4 and B ∈ [4,5]µ.Furthermore, the magenta curve fits very well with all observational data, but its maximum-mass value is above 3M .
Another interesting feature of these stars is their compactness, defined by C ≡ M/R.According to the classification adopted by Iyer et al. [92], the configurations shown in the mass-radius diagram correspond to compact stars, see the right plot of Fig. 1.Besides, we can appreciate that the compactness of dark energy stars is of the order of the compactness of hadronic-matter stars, as is the case of the SLy EoS [93], despite the fact that the maximum mass in the magenta configuration sequence can exceed 3M .Nonetheless, as we will see later, the introduction of anisotropy can turn such stars into ultracompact objects.Of course, this will depend on the amount of anisotropy in the stellar interior.
In order to include anisotropic pressures and investigate their effects on the internal structure of dark energy stars, we will adopt two specific models with the following parameters Similar to the isotropic case, we numerically solve the hydrostatic background equations ( 4)-( 6) with boundary conditions (7), but taking into account the anisotropy profile (30).For instance, for the model I and a central density ρ c = 2.0 × 10 18 kg/m 3 , Fig. 2 illustrates the mass density, pressure and squared speed of sound as functions of the radial coordinate for different values of the free parameter α.We can see that the internal structure of a dark energy star is affected by the presence of anisotropy.In effect, the radius of the star increases (decreases) for more positive (negative) values of α.In addition, we remark that the speed of sound, both radial and tangential, respect the causality condition.This has also been verified for other values of central density considered in the construction of Fig. 1.
Varying the central density, we obtain the mass-radius diagrams and mass-central density relations for models I and II, as shown in Fig. 3.We observe that the substantial changes introduced by anisotropy in dark energy stars occur in the high-mass branch (close to the maximum-mass point), while the effects are irrelevant at low central densities.The maximum-mass values increase as the parameter α increases (see also the data in Table I).Note that model I without anisotropic pressures is not capable of generating maximum masses above 2M .Here the constant B is given in µ = 10 −20 m −4 units.The gray horizontal stripe at 2.0M stands for the two massive NS pulsars J1614-2230 [85] and J0348+0432 [86].Yellow and blue regions represent the observational measurements of the masses of the highly massive NS pulsars J0740+6620 [87] and J2215+5135 [94], respectively.The filled pink band stands for the lower mass of the compact object detected by the GW190814 event [91], and the cyan area is the mass-radius constraint from the GW170817 event.Moreover, the NICER measurements for PSR J0030+0451 are displayed by black dots with their respective error bars [95,96].Right panel: Variation of the compactness with total gravitational mass, where the gray and orange stripes represent compact and ultra-compact objects, respectively, according to the classification given in Ref. [92].For comparison reasons, we have included the results corresponding to the SLy EoS [93] by blue curves in both plots.
Nevertheless, the inclusion of anisotropies (see the blue curve for α = 0.4) allows a significant increase in the maximum mass and hence a more favorable description of the compact objects observed in nature.On the other hand, model II with anisotropies (see orange curves) fits better with the observational measurements.In particular, in view of the lower mass of the compact object from the coalescence GW190814 [91], two curves are particularly outstanding.In other words, such object can be well described as an anisotropic dark energy star when α = 0.2 and α = 0.4.Moreover, model II with negative anisotropies (such as α = −0.4)favors the description of the massive pulsar J2215+5135 [94].
The left panel of Fig. 4 describes the behavior of compactness as a function of central density.Positive anisotropies lead to an increase in compactness, mainly in the high-central-density branch.Remarkably, for sufficiently large values of α (see purple curve), it is possible to obtain anisotropic dark energy stars as ultra-compact objects.
The gravitational redshift, conventionally defined as the fractional change between observed and emitted wavelengths compared to emitted wavelength, in the case of a Schwarzschild star is given by [61] In the right plot of Fig. 4, the surface gravitational red- shift is plotted as a function of the total mass for both models I and II.This plot indicates that the gravitational redshift of light emitted at the surface of a dark energy star is substantially affected by the anisotropy in the high-mass region, while the changes are negligible for sufficiently low masses.For a fixed value of central density, Table II shows that positive (negative) anisotropy increases (decreases) the value of the redshift.All plots correspond to model I and the black curves represent the isotropic solutions.Note that both the radial and tangential speed of sound obey the causality condition.Furthermore, one can observe that the increase in α leads to larger radii, and the anisotropy is more pronounced in the intermediate regions.

B. Oscillation spectrum
A necessary condition (the well-known M (ρ c ) method) for stellar stability is that stable stars must lie in the region where dM/dρ c > 0. According to the right plot of Fig. 3, the full blue and orange circles on each curve indicate the onset of instability for each family of equilibrium solutions.However, a sufficient condition is to calculate the frequencies of the radial vibration modes for each central density [61].Here we will analyze if both methods are compatible in the case of dark energy stars including anisotropic pressure.
Once the equilibrium equations ( 4)-( 6) are integrated from the center to the surface of the star, we then proceed to solve the radial pulsation equations ( 8) and ( 9) with the corresponding boundary conditions ( 10) and ( 11) using the shooting method.Namely, we integrate from the origin (where we consider the normalized eigenfunctions ζ(0) = 1) up to the stellar surface for a set of trial values ν 2 satisfying the condition (10).In this way, the appropriate eigenfrequencies correspond to the values for which the boundary condition ( 11) is fulfilled.For instance, for a central density ρ c = 1.5 × 10 18 kg/m 3 , α = 0.4 and parameters given by model I, Fig. 5 displays the radial behavior of the perturbation variables for the first five squared eigenfrequencies ν 2 n , where n indicates the number of nodes inside the star.This frequency spectrum forms an infinite discrete sequence, i.e.
, where the eigenvalue corresponding to n = 0 is the lowest one (or equivalently, the longest period of all the allowed vibration modes) and it is known as the fundamental mode.Such mode has no nodes,  I).The critical central density corresponding to the maximum point on the M (ρc) curve is modified by the presence of anisotropy for both models.The gray and light-green stripes represent compact and ultra-compact objects, respectively, according to the classification established by Iyer et al. [92].Positive anisotropy results in increased compactness for sufficiently high central densities, while the opposite occurs for negative anisotropy.Note also that dark energy stars would correspond to ultra-compact objects if α > 0.4 for model II, see for instance the purple curve for α = 0.7.Right panel: Surface gravitational redshift as a function of the total mass.In the high-redshift region it can be observed that positive (negative) anisotropy increases (decreases) the value of zsur.Meanwhile, the effect of anisotropy is irrelevant for sufficiently low redshifts.
whereas the first overtone (n = 1) has one node, the second overtone (n = 2) has two, and so on.Stable stars are described by their oscillatory behavior so that ν 2 n > 0 (i.e., ν n is purely real).On the other hand, if any of these is negative for a particular star, the frequency is purely imaginary and hence the star is unstable.
Since each higher-order mode has a squared eigenfrequency that is larger than in the case of the preceding mode, it is enough to calculate the frequency of the fun-damental pulsation mode for the equilibrium sequences presented in Fig. 3.With this in mind, in Fig. 6 we plot the squared frequency of the fundamental oscillation mode as a function of the central density (left panel) and gravitational mass (right panel).According to the left plot, the squared frequency of the fundamental mode is exactly zero at the critical-central-density value corresponding to the maximum-mass configuration as shown in the right plot of Fig. 3, see the full blue and orange cir-cles for both models.Furthermore, according to the right plot of Fig. 6, the maximum-mass values (that is, when dM/dρ c = 0) can be used as turning points from stability to dynamical instability.Therefore, we can conclude that the usual criterion to guarantee stability dM/dρ c > 0 is still valid for the case of anisotropic dark energy stars.In other words, the conventional M (ρ c ) method is compatible with the calculation of the eigenfrequencies of the normal vibration modes.
If the anisotropic dark energy star has a central density higher than one corresponding to the maximum-mass configuration (indicated by full blue and orange circles in Figs. 3 and 6), the star will become unstable against radial perturbations and collapse to form a black hole.For further details on the dissipative gravitational collapse of compact stellar objects we also refer the reader to Refs.[55,[97][98][99].Nonetheless, we must point out that there are EoS models that allow a compact star to migrate to another branch of stable solutions instead of forming a black hole when it is subjected to a perturbation.As a matter of fact, the first-order phase transition between nuclear and quark matter can generate multiple stable branches in the mass-radius diagram for hybrid stars [100].

C. Moment of inertia
To calculate the moment of inertia of anisotropic dark energy stars, we first need to solve the differential equation for the rotational drag (16) with boundary conditions (19).In particular, for model I and central density ρ c = 1.5 × 10 18 kg/m 3 , figure 7 illustrates the angular velocity everywhere for several values of α.As can be observed in the right plot, the dragging angular velocity outside the star has the behavior ω(r) ∼ r −3 , so that at infinity (where spacetime is flat) the distant local inertial frames do not rotate around the star, namely, ω(r) → 0 for r → ∞.Moreover, anisotropy significantly affects the angular velocity of the local inertial frames in the interior region of the star.More specifically, the dragging angular velocity increases (decreases) for positive (negative) values of the anisotropy parameter α.We can then determine the moment of inertia using the integral given in Eq. (18).For the above central density, we present the moment of inertia of some dark energy configurations for both models in Table II, where it can be noticed that I increases as the value of α increases.
We can now calculate the moment of inertia for a whole sequence of dark energy stars by varying the central density ρ c .The left panel of Fig. 8 displays the moment of inertia as a function of the gravitational mass for both models.Remarkably, model II provides larger values for the moment of inertia than model I. Indeed, the maximum value I max depends quite sensitively on the free parameters A and B in the EoS (28).In addition, the main effect of anisotropy on the moment of inertia for slow rotation occurs in the high-mass region, while its influence is irrelevant for sufficiently low masses.In order to better quantify the changes in the maximum values of the moment of inertia induced by the anisotropic pressure, we can define the following relative difference ∆I = I max,ani − I max,iso I max,iso , where I max,iso and I max,ani are the maximum values of the moment of inertia for isotropic and anisotropic configurations, respectively.In the right plot of Fig. 8 we present the dependence ∆I against the anisotropy parameter α.
The impact of anisotropy is getting stronger as |α| grows, reaching variations (with respect to the isotropic case) of up to ∼ 20% for α = 0.5.We can also note that such relative variations are almost independent of the model adopted.

D. Tidal properties
We will now investigate how the anisotropy parameter α affects the tidal properties of dark energy stars.Given a specific value of α, this requires solving the differential equation (23) for a range of central densities.The left panel of Fig. 9 is the result of calculating the tidal Love number (27) for a sequence of stellar configurations by considering different values of α, where the isotropic case corresponds to α = 0. Similar to the trends in strange quark stars, as reported in Ref. [70], the Love number of dark energy stars grows until it reaches a maximum value and then decreases as compactness increases.Note also that the maximum value of k 2 is sensitive to the value of α, indicating that the Love number decreases as the parameter α increases for both models.Although model II provides larger maximum masses (as well as redshift and moment of inertia) than model I, we see that the behavior is different for the maximum values in the tidal Love number.
Ultimately, in the right plot of Fig. 9, the dimensionless tidal deformability Λ = λ/M 5 is plotted as a function of mass, where it can be observed that smaller masses yield higher deformabilities.In each model, the presence of anisotropy has a negligible effect on Λ for small masses, while slightly more significant changes take place only in the high-mass region.

V. CONCLUSIONS AND OUTLOOK
In this work, we have focused on the equilibrium structure of dark energy stars by using a Chaplygin-like equation of state under the presence of both isotropic and anisotropic pressures within the context of standard GR.Our goal was to construct stable compact stars whose characteristics could be compared with the observational data on the mass-radius diagram.In this perspective, the global properties of a compact star such as radius, mass, redshift, moment of inertia, oscillation spectrum and tidal deformability have been calculated.To describe the anisotropic pressure within the dark energy fluid we have adopted the anisotropy profile proposed by Horvat et al. [51], where a free parameter α measures the degree of anisotropy.
We have discussed the possibility of observing stable dark energy stars made of a negative pressure fluid "−B/ρ" plus a barotropic component "Aρ".By way of comparison, the EoS parameters A and B have been cho-sen in such a way that they agree sufficiently with the observational data, e.g. the mass-radius constraint from the GW170817 event.For isotropic configurations, we have shown that various sets of values {A, B} can be chosen since they obey the causality condition and consistently describe compact stars observed in the Universe.Furthermore, we saw that the secondary component resulting from the gravitational-wave signal GW190814 [91] can be described as a dark energy star using A = 0.4 and  16) for a dark energy star described by model I and central ρc = 1.5 × 10 18 kg/m 3 in the presence of anisotropy for several values of the free parameter α.The solid and dashed lines represent the interior and exterior solutions, respectively.Right panel: Ratio of frame-dragging angular velocity to the angular velocity of the star, namely ω(r)/Ω = 1 − (r)/Ω.It can be observed that the outer solution behaves asymptotically at large distances from the surface of the star (this is, ω → 0 for r → ∞).Furthermore, appreciable changes in the angular velocity due to anisotropy can be noticeable, mainly in the interior region of the star.Based on these results, we have established two models with different values A and B in order to explore the effects of anisotropy in the interior region of a dark energy star.In particular, the maximum-mass values increase as the parameter α increases.We noticed that model I without anisotropic pressures is not capable of generating maximum masses above 2M .However, the inclusion of anisotropies (α = 0.4) allows a significant increase in the maximum mass and thus a more favorable description of the compact objects observed in nature.On the other hand, model II with anisotropies fits better with the observational measurements, although such a model can lead to the formation of ultra-compact objects for sufficiently large values of α.We also calculated the surface gravitational redshift for such stars, and our results indicated that z sur is substantially affected by the anisotropy in the high-mass branch, while the changes are irrelevant for sufficiently low masses.
A star exists in the Universe only if it is dynamically stable, so our second task was to investigate whether the dark energy stars are stable or unstable with respect to an adiabatic radial perturbation.Our results showed that the standard criterion for radial stability dM/dρ c > 0 still holds for dark energy stars since the squared frequency of the fundamental pulsation mode (ν 2 0 ) van- Note also that the Love number is substantially modified by the anisotropy parameter α for both models, while its greatest effect on tidal deformability Λ occurs only in the high-mass region.
ishes at the critical central density corresponding to the maximum-mass configuration.This has been examined in detail for both isotropic (α = 0) and anisotropic (α = 0) stellar configurations.
In the slowly rotating approximation, where only firstorder terms in the angular velocity are kept, we have also determined the moment of inertia of anisotropic dark energy stars.For this purpose, we first had to calculate the frame-dragging angular velocity for each central density.The presence of anisotropic pressure results in a substantial increase (decrease) of the angular velocity ω for more positive (negative) values of α.We found that the significant impact of the anisotropy on the moment of inertia occurs mainly in the high-mass branch for both models.Furthermore, the maximum value of the moment of inertia can undergo variations of up to ∼ 20% for α = 0.5 as compared with the isotropic case.
We have analyzed the effect of anisotropic pressure on the tidal properties of such stars.In particular, our outcomes revealed that the tidal Love number is sensitive to moderate variations of the parameter α, indicating that the maximum value of k 2 can increase as α decreases.In addition, the greatest effect of anisotropy on the dimensionless tidal deformability takes place only in the high-mass region.Based on the foregoing results, the present work thereby serves to develop a comprehensive perspective on the relativistic structure of dark energy stars in the presence of anisotropy.
Summarizing, we have explored the possible existence of stable dark energy stars whose masses and radii are not in disagreement with the current observational data.The Chaplygin-like EoS predicts maximum-mass values consistent with observational measurements of highly massive pulsars.Future research includes the adoption of widespread versions of Chaplygin gas that best fit key cosmological parameters.In future studies we will thereby take further steps in that direction, focusing on the different types of generalized Chaplygin gas models as discussed in Ref. [11].In addition, as carried out in the case of boson stars [101], it would be interesting to employ a Fisher matrix analysis in order to distinguish dark energy stars from black holes and neutron stars from tidal interactions in inspiraling binary systems.It is also worth mentioning that Romano [102] has recently discussed the effects of dark energy on the propagation of gravitational waves.In that regard, we expect that future electromagnetic observations of compact binaries and gravitationalwave astronomy will provide a better understanding of compact stars in the presence of dark energy, and even help us answer the most basic question: How did dark energy form in the Universe?Anyway, our results suggest that dark energy stars deserve further investigation by taking into account the cosmological aspects as well as the gravitational-wave signals from binary mergers.
Model I: A = 0.3, B = 6.0µ , Model II: A = 0.4, B = 5.2µ , which are models favored by observational measurements according to the left panel of Fig. 1.Moreover, model II precisely corresponds to the first model considered by Panotopoulos et al. [38].

SLyAFIG. 1 .
FIG.1.Left panel: Mass-radius diagrams for dark energy stars with Chaplygin-like EoS(28) and isotropic pressure (σ = 0) for several values of the positive parameters A and B.Here the constant B is given in µ = 10 −20 m −4 units.The gray horizontal stripe at 2.0M stands for the two massive NS pulsars J1614-2230[85] and J0348+0432[86].Yellow and blue regions represent the observational measurements of the masses of the highly massive NS pulsars J0740+6620[87] and J2215+5135[94], respectively.The filled pink band stands for the lower mass of the compact object detected by the GW190814 event[91], and the cyan area is the mass-radius constraint from the GW170817 event.Moreover, the NICER measurements for PSR J0030+0451 are displayed by black dots with their respective error bars[95,96].Right panel: Variation of the compactness with total gravitational mass, where the gray and orange stripes represent compact and ultra-compact objects, respectively, according to the classification given in Ref.[92].For comparison reasons, we have included the results corresponding to the SLy EoS[93] by blue curves in both plots.

2 FIG. 2 .
FIG.2.Radial behavior of the mass density (left panel), pressures (middle panel) and squared speed of sound (right panel) inside an anisotropic dark energy star with central density ρc = 2.0 × 1018 3 and several values of the parameter α.All plots correspond to model I and the black curves represent the isotropic solutions.Note that both the radial and tangential speed of sound obey the causality condition.Furthermore, one can observe that the increase in α leads to larger radii, and the anisotropy is more pronounced in the intermediate regions.

FIG. 3 .
FIG. 3. Mass-radius diagram (left panel) and mass-central density relation (right panel) for anisotropic dark energy stars as predicted by model I (blue curves) and II (orange curves) with anisotropy profile (30) for several values of α.The colored bands in the left plot represent the same as in Fig. 1.Moreover, the full blue and orange circles on the right plot indicate the maximum-mass points for model I and II, respectively.Note that the maximum-mass values for model II correspond to lower central densities than those for model I, however, model II allows larger masses (see also TableI).The critical central density corresponding to the maximum point on the M (ρc) curve is modified by the presence of anisotropy for both models.

FIG. 4 .
FIG.4.Left panel: Variation of the compactness with central density for several anisotropic dark energy star sequences.The gray and light-green stripes represent compact and ultra-compact objects, respectively, according to the classification established by Iyer et al.[92].Positive anisotropy results in increased compactness for sufficiently high central densities, while the opposite occurs for negative anisotropy.Note also that dark energy stars would correspond to ultra-compact objects if α > 0.4 for model II, see for instance the purple curve for α = 0.7.Right panel: Surface gravitational redshift as a function of the total mass.In the high-redshift region it can be observed that positive (negative) anisotropy increases (decreases) the value of zsur.Meanwhile, the effect of anisotropy is irrelevant for sufficiently low redshifts.

FIG. 5 . 2 ]FIG. 6 .
FIG.5.Numerical solution of the radial pulsation equations (8) and(9) in the case of an anisotropic dark energy star with central density ρc = 1.5 × 10 18 kg/m 3 , α = 0.4 and EoS parameters given by model I.The radius, mass and the fundamental mode frequency for such configuration are found in TableII.The lines with different colors and styles indicate different overtones so that the solution corresponding to the nth vibration mode contains n nodes in the internal structure of the star.Note that the eigenfunctions ζn(r) have been normalized assuming ζ = 1 at r = 0, and the Lagrangian perturbation of the radial pressure ∆pr,n(r) obeys the boundary condition(11) at the stellar surface.Since f0 is real, this configuration corresponds to a stable anisotropic dark energy star.

FIG. 7 .
FIG.7.Left panel: Numerical solution of the differential equation (16) for a dark energy star described by model I and central ρc = 1.5 × 10 18 kg/m 3 in the presence of anisotropy for several values of the free parameter α.The solid and dashed lines represent the interior and exterior solutions, respectively.Right panel: Ratio of frame-dragging angular velocity to the angular velocity of the star, namely ω(r)/Ω = 1 − (r)/Ω.It can be observed that the outer solution behaves asymptotically at large distances from the surface of the star (this is, ω → 0 for r → ∞).Furthermore, appreciable changes in the angular velocity due to anisotropy can be noticeable, mainly in the interior region of the star.

FIG. 8 .
FIG.8.Left panel: Moment of inertia versus mass for anisotropic dark energy stars, where a higher mass results in larger moment on inertia for both models.It is observed that the substantial impact of anisotropy on the moment of inertia occurs predominantly in the high-mass branch.Right panel: Relative deviation (33) as a function of the anisotropy parameter.The maximum value of the moment of inertia can undergo variations with respect to its isotropic counterpart of up to ∼ 20% for α = 0.5.

FIG. 9 .
FIG.9.Left panel: Tidal Love number plotted as a function of the compactness C ≡ M/R.Right panel: Dimensionless tidal deformability versus gravitational mass predicted by each model, where larger masses yield smaller deformabilities.Note also that the Love number is substantially modified by the anisotropy parameter α for both models, while its greatest effect on tidal deformability Λ occurs only in the high-mass region.

TABLE II .
Radius, mass, redshift, fundamental mode frequency (f0 = ν0/2π), moment of inertia and dimensionless tidal deformability of dark energy stars with central energy density ρc = 1.5 × 10 18 kg/m 3 as predicted by models I and II for several values of the anisotropy parameter α.Remarkably, with the exception of the fundamental mode frequency and tidal deformability, these properties undergo a significant increase as α increases.