Nonradial oscillation modes of a neutron star possess valuable information about its internal structure and nuclear physics. Starting from the quadrupolar order such modes under general relativity are known as quasinormal modes since they dissipate energy through gravitational radiation and their frequencies are complex. The stability of these modes is governed by the sign of the imaginary part of the frequency which determines whether the mode would decay or grow over time. In this paper we develop a fully consistent framework in general relativity to study quasinormal modes of neutron stars with anisotropic pressure whose motivation includes strong internal magnetic fields and nonvanishing shear or viscosity. We employ parametrized models for the anisotropy with a specific form under perturbation that is consistent with previous studies under the Cowling approximation and solve the perturbed Einstein field equations numerically. The assumed form of the perturbed anisotropy can be easily extended to more general cases. We find that unlike the case for isotropic neutron stars the imaginary parts of some of the pressure inlineformulammlmath displayinlinemmlmipmmlmimmlmathinlineformulamodes flip signs as the degree of anisotropy deviates from zero depicting a transition from stable modes to unstable modes. This finding indicates that some anisotropic neutron star models are unstable potentially restricting the form of sustained anisotropy.
['Astrophysics and astroparticle physics', 'General Relativity and Quantum Cosmology', 'Astrophysics - High Energy Astrophysical Phenomena', 'High Energy Physics - Phenomenology']
Nonradial instabilities in anisotropic neutron stars
{'Nonradial instabilities in anisotropic neutron stars': 'Shu Yan Lau, ∗ Siddarth Ajith, Victor Guedes, and Kent Yagi † Department of Physics, University of Virginia, Charlottesville, Virginia 22904, USA (Dated: September 17, 2024)  \nNon-radial oscillation modes of a neutron star possess valuable information about its internal structure and nuclear physics. Starting from the quadrupolar order, such modes under general relativity are known as quasi-normal modes since they dissipate energy through gravitational radiation and their frequencies are complex. The stability of these modes is governed by the sign of the imaginary part of the frequency, which determines whether the mode would decay or grow over time. In this paper, we develop a fully consistent framework in general relativity to study quasi-normal modes of neutron stars with anisotropic pressure, whose motivation includes strong internal magnetic fields and non-vanishing shear or viscosity. We employ parametrized models for the anisotropy, with a specific form under perturbation that is consistent with previous studies under the Cowling approximation, and solve the perturbed Einstein field equations numerically. The assumed form of the perturbed anisotropy can be easily extended to more general cases. We find that, unlike the case for isotropic neutron stars, the imaginary parts of some of the pressure ( p -)modes flip signs as the degree of anisotropy deviates from zero, depicting a transition from stable modes to unstable modes. This finding indicates that some anisotropic neutron star models are unstable, potentially restricting the form of sustained anisotropy.', 'I. INTRODUCTION': "Pulsation modes provide the basis to study the dynamical deformations of a neutron star (NS). The frequency of these modes characterizes the stellar internal composition and density profile. In addition, the stability of the modes also helps predict whether certain configurations of a star can exist in nature. With the advancement in multimessenger astronomy over the years, the detection of some of the pulsation modes is expected in the future [1-8].  \nWithin Newtonian theory, the adiabatic pulsation modes are normal modes described by real frequencies in the absence of dissipation, radiation reaction, and various types of instabilities (see, e.g., [9-14]). For even parity (polar) non-radial modes, the well-known ones include the fundamental ( f -)modes, the pressure ( p -)modes, and the gravity ( g -)modes. In the original classification by Cowling [15], the naming of these modes is based on the local propagation behavior of the waves that form the modes. For example, the p -modes, which are of interest in this study, can be conceived as the standing waves formed by propagating acoustic waves within the bulk via a local analysis [13, 16, 17].  \nUnder general relativity (GR), non-radial oscillation modes (except the dipole modes) can dissipate energy through the emission of gravitational radiation. These modes are described by oscillations with complex frequencies, denoted by ω hereafter, and are known as the quasinormal modes (QNMs). The complex values are the characteristics of an open system. Despite many of those being referred to as the 'fluid modes' [18], which are analogous to the Newtonian stellar pulsation modes, they are  \nfundamentally different due to the boundary conditions in the exterior vacuum region of the star.  \nThe stability of the QNMs is determined by the sign of the imaginary part [19]. For our sign convention of the time dependence of the QNMs expressed as e iωt , the mode will grow exponentially over time if the imaginary part is negative, thus being unstable.  \nTypical NSs are described as isotropic fluids, in which the local stress is invariant under rotational transformations. Stability analysis on the QNMs of non-rotating isotropic fluid stars has shown that the unstable modes must have pure imaginary frequencies, i.e., the mode must be non-propagating [20-22]. This stability criterion can be extended to certain regimes, including the two-fluid NS models [23], where the argument given in [20] still applies. A similar condition has also been shown to hold in the scalar-tensor theories for the radial scalarmodes [24, 25], an example of QNMs in non-GR theories.  \nMeanwhile, local anisotropy can exist within a NS in various scenarios involving elasticity [26-28], superfluidity [29, 30], pion condensation [31], strong magnetic field [32-36], and viscosity [37]. Certain exotic relativistic stars are also predicted to have anisotropic stress, like boson stars [38-40], strange stars [41, 42], dark matter stars [43], and gravastars [44, 45]. The solution to a static spherically symmetric anisotropic NS was first studied in [46]. After that, much work has been done on calculations of the NS structure with different anisotropic equations of state (EOSs) as well as assessing the radial stability [47-55]. One interesting finding is that depending on the anisotropy analysis, the maximum-mass configuration of the anisotropic NSs may or may not correspond to a change in radial stability [53-55], in contrast to the isotropic NSs [56].  \nThe non-radial deformations of anisotropic NSs are also studied in previous literature. In the static limit, such perturbations within the slow rotation and small  \ntidal deformation approximation allow one to calculate the moment of inertia, spin-induced quadrupole moment, and tidal deformability [57, 58] that are known to enjoy some universal relations [59-61]. Anisotropic NSs were used to also study how universal relations for multipole moments of NSs approach the black hole limit [62, 63]. For non-radial time-dependent perturbations, Refs. [6467] employed the relativistic Cowling approximation in calculating the pulsation modes of anisotropic NSs by ignoring the metric perturbations. Such an approximation can accurately predict the frequencies of non-radial oscillation modes for isotropic stars with an error of < 10% for p -modes, which decreases with the order of the mode [68], while its validity is currently unknown for anisotropic stars. Reference [48] has also investigated the non-radial instability of anisotropic stars, demonstrating that the non-radial pulsation modes are stable except for infinitely large anisotropy. They, however, considered only the Newtonian limit of incompressible stars for a specific anisotropy ansatz that turns out to be unphysical as it causes a divergence at the stellar center.  \nIn this work, we develop a consistent framework to treat stellar perturbations for anisotropic fluids in full GR. Although there are many works studying stellar perturbations for isotropic stars in full GR since the early 1970s [20-22], there was no such formulation available for anisotropic stars for more than 50 years. The new formulation developed here allows us to check the validity of the Cowling approximation and study the stability of anisotropic stars for the first time. Our numerical result shows that the anisotropy causes certain p -modes to acquire a negative imaginary part in the eigenfrequencies while having a non-zero real part. This can be interpreted as an oscillatory mode with exponentially growing amplitude, which is in stark contrast with the case for isotropic NSs. We further demonstrate how the mathematical proof of the non-existence of oscillatory unstable modes for isotropic NSs does not apply to the anisotropic case by considering a similar integral formula as Eq. (16) of [20]. We also found that the Cowling approximation is valid for anisotropic NSs with an error of < 30% for f -modes and < 10% for p -modes from the fully consistent analysis.  \nThe rest of this paper is organized as follows: In Sec. II, we describe the two anisotropy models we employ in this study. Next, we give the formalism for computing the non-radial pulsations of a spherically symmetric nonrotating anisotropic NS in Sec. III. We then show the results of the mode instability in Sec. IV. Lastly, we discuss the interpretations of the results and possible extensions of the work in Sec. V.  \nUnless otherwise specified, we employ the geometrized unit system with G = c = 1 in this paper.", 'II. ANISOTROPIC NEUTRON STARS': "The fluid stress within an anisotropic NS has directional dependence. The stress-energy tensor describing the anisotropic fluid with spherical symmetry takes the form [46]  \nT αβ = ρu α u β + p r h αβ -σ Ω αβ , (1)  \nwhere ρ is the energy density, σ = p r -p t is the degree of anisotropy, p r and p t represent the pressures in the radial and tangential direction respectively, h αβ = g αβ + u α u β is the transverse metric on a 3D space, Ω αβ = h αβ -k α k β is another transverse metric on a 2D sphere, u α is the four-velocity vector of a fluid, and k α is the unit normal vector in the radial direction orthogonal to u α .  \nIn this paper, we consider two anisotropy ansatzes described by quasi-local parameters. These anisotropy models contain not only the local thermodynamic parameters, like ρ , or p r , but also quasi-local parameters that depend on the geometry at the specified point of the spacetime. The first one is modified from the anisotropy model given in Horvat et al. [51]:  \nσ = βp r µ 2 , (2)  \nwhere µ = 2 m/r , m is the gravitational mass enclosed within a radius r , and β is a dimensionless parameter governing the degree of anisotropy. We denote this ansatz as the 'H-model'. Note that compared to [51], σ here contains an extra factor of µ . Such an additional factor is introduced to satisfy the regularity conditions of the perturbations at r = 0, which requires ∂σ/∂µ = O ( r 2 ) as shown in Eq. (31). The above model still satisfies the conditions used to construct the original model, namely (i) the anisotropy term to vanish from the hydrostatic equilibrium equation in the non-relativistic limit ( p r ≪ ρ ), and (ii) the anisotropy vanishes at the stellar center. Due to the additional factor of µ in the new model which is of the order ∼ 0 . 1, we use the range | β | ≤ 10 instead of | β | ≤ 2 in the original model. Note that in the relativistic Cowling approximation used in [64], the original anisotropy model satisfies the regularity conditions automatically since there are no metric perturbations.  \nThe second model is proposed by Bowers and Liang [46]:  \nσ = β ( ρ + p r )( ρ +3 p r ) 1 -µ r 2 , (3)  \nwhere | β | ≤ 0 . 5 as in [55]. This anisotropy model permits an analytic solution to the stellar structure for incompressible stars. We denote it as the 'BL-model'. Unlike the H-model, we do not need to modify the BL-model as it satisfies the regularity conditions at the center.  \nThe relation between ρ and p r is provided by an EOS. We choose the polytropic EOS named 'EOS II' in [64] (denoted by 'Poly' hereafter) and a more realistic EOS, SLy4 [69], to demonstrate our major numerical results.  \nWith these anisotropy ansatzes and EOSs, we numerically integrate the structural equations (see, e.g., [46]) to construct the anisotropic NS models at equilibrium. The stress-energy tensor of all models considered satisfies the dominant, weak, and strong energy conditions. Moreover, the speeds of sound in the radial and tangential directions within the stars (see, e.g., [66, 70]) both satisfy causality. We also check that both p r and p t are positive. Note that for some models with large positive β outside the range specified above, p t can become negative near the stellar surface, which is defined at the radius where p r = 0.", 'III. FORMULATION': "Here we provide the full GR formalism for computing the QNMs for the first time. We construct the background non-rotating solution of an anisotropic star and introduce non-radial time-dependent small perturbations to the metric and fluid variables. We numerically solve the equations to obtain the QNM frequencies and eigenfunctions. For the perturbations, we assume they follow a specific form (given in Eq. (27) described below), which reduces to the Cowling approximation [64] when we take the metric perturbations to zero. Note that this is not derived from the fundamental properties of anisotropic fluids. One important assumption of this model is that the perturbed tangential pressure is always the same for both directions along the 2-sphere. We suggest a further generalization in Appendix A but leave it for future work.  \nWe first define the static spherically symmetric background metric as  \nds 2 = -e ν dt 2 + e λ dr 2 + r 2 d Ω 2 , (4)  \nwhere ν and λ are functions of r while d Ω is the unit 2sphere line element. The Einstein field equations of the anisotropic fluid background give (e.g., [46])  \nν ' =2 m +4 πr 3 p r r 2 e λ , (5)  \np ' r = -( ρ + p r ) ν ' 2 -2 σ r , (6)  \nλ ' =2 -m +4 πr 3 ρ r 2 e λ , (7)  \ne λ = ( 1 -2 m r ) -1 , (8)  \nwhere the superscript prime denotes derivative with respect to r .  \nIn the pulsating star, the metric and the matter dynamics are governed by the perturbed Einstein field equa-  \ntions and the stress-energy conservation:  \nδG αβ = 8 πδT αβ , (9)  \nδ ( ∇ α T αβ ) = 0 , (10)  \nwhere ∇ α denotes the covariant derivative, G αβ is the Einstein tensor and T αβ is the stress-energy tensor. The prefix δ represents Eulerian perturbations and is related to the Lagrangian perturbations ∆ through the Lie derivative [71], ∆ = δ + L ξ , given the Lagrangian displacement vector ξ α of the fluid.  \nWe employ the same convention of the perturbed metric as in [72, 73], defined in the Regge-Wheeler gauge [74]:  \nδg αβ dx α dx β = -∑ ℓ,m [ e ν r ℓ H 0 dt 2 +2 iωr ℓ +1 H 1 dtdr + e λ r ℓ H 2 dr 2 + r ℓ +2 Kd Ω 2 ] e iωt Y ℓm . (11)  \nHere ( H 0 , H 1 , H 2 , K ) are functions of r , and Y ℓm are the spherical harmonics. The radial parts of the expansion depend on the index ℓ , which is suppressed for conciseness in the following expressions. This applies to all the radial parts of the expansions in spherical harmonics. Since we focus on non-rotating stars with spherical backgrounds, perturbations with different ℓ do not couple with each other.  \nThe stress-energy tensor takes the same form as [64]. One major assumption about the perturbed density and radial pressure is the adiabatic relation  \n∆ p r = γp r ρ + p r ∆ ρ, (12)  \nwhere γ is the adiabatic index, and is assumed to be  \nγ = ρ + p r p r ( ∂p r ∂ρ ) s = ρ + p r p r ( ∂p r ∂ρ ) eq . (13)  \nThe subscript ' s ' represents an adiabatic derivative while the subscript 'eq' denotes the derivative of the background p r against ρ . This essentially means the EOS remains a single variable function in the dynamical timescale. The Lagrangian displacement vector is defined as  \nξ α ∂ α = ∑ ℓ,m r ℓ -1 e iωt [ e -λ/ 2 WY ℓm ∂ r -V r ( ∂Y ℓm ∂θ ∂ θ +csc 2 θ ∂Y ℓm ∂ϕ ∂ ϕ ) ] , (14)  \nwhere V and W are functions of r .  \nUsing Eqs. (9)-(14), we arrive at the set of equations governing the stellar pulsations:  \nH ' 1 = [ 4 π ( ρ -p r ) e λ r -2 me λ r 2 -ℓ +1 r ] H 1 + e λ r K + e λ r H 0 -16 π ( ρ + p r ) r e λ (1 -¯ σ ) V, (15)  \nK ' = ℓ ( ℓ +1) 2 r H 1 + ( ν ' 2 -ℓ +1 r ) K -8 π ( ρ + p r ) e λ/ 2 r W + 1 r H 0 , (16)  \nW ' = re λ/ 2 (1 -¯ σ ) K + ( -ℓ +1 r + 2¯ σ r ) W + re ( λ -ν ) / 2 γp r X + re λ/ 2 2 H 0 -ℓ ( ℓ +1) r e λ/ 2 (1 -¯ σ ) V, (17)  \nX ' = ρ + p r 2 e ν/ 2 [ rω 2 e -ν + ℓ ( ℓ +1) 2 r (1 -2¯ σ ) ] H 1 + ρ + p r 2 e ν/ 2 [( 3 2 -2¯ σ ) ν ' -(1 -6¯ σ ) 1 r -4¯ σ 2 r ] K -ρ + p r r e ( λ + ν ) / 2 [ 4 π ( ρ + p r ) + ω 2 e -ν -F ] W -1 r ( ℓ -2 ρ + p r γp r ¯ σ ) X + ρ + p r 2 e ν/ 2 ( 1 r -ν ' 2 ) H 0 + ℓ ( ℓ +1) e ν/ 2 p ' r r 2 (1 -¯ σ ) V + 2 e ν/ 2 r S. (18)  \nHere, the overbarred anisotropy parameter is defined by  \n¯ σ = σ ρ + p r , (19)  \nwhile X is the radial part of the expansion  \n∑ ℓ,m r ℓ XY ℓm e iωt = -e ν/ 2 ∆ p r , (20)  \nand F in Eq. (18) is given by  \nF = e -λ/ 2 { r 2 2 ( e -λ/ 2 ν ' r 2 ) ' -e -λ/ 2 [( 6 r 2 -2 ν ' r ) ¯ σ -2 r 2 ( rσ ' ρ + p r ) -4 r 2 ¯ σ 2 ] } . (21)  \nThe Eulerian perturbation of the anisotropy parameter, S , is defined by  \nand is related to other perturbation variables by  \n∑ ℓ,m r ℓ SY ℓm e iωt = δσ, (22)  \nS = -[ ( ∂σ ∂p r ) eq +( ρ + p r ) ( A p ' r + 1 γp r )( ∂σ ∂ρ ) eq ] e -λ/ 2 p ' r r W -[ ( ∂σ ∂p r ) eq + ρ + p r γp r ( ∂σ ∂ρ ) eq ] e -ν/ 2 X -( ∂σ ∂µ ) eq e -λ H 0 . (23)  \nHere, A is the relativistic Schwarzschild discriminant, which we take to be zero due to our assumption on γ (i.e., assuming the EOS is a single variable function of p r ), namely  \nA ≡ ρ ' ρ + p r -p ' r γp r = 0 . (24)  \nWe further assume the partial derivatives on σ follow the equilibrium background [64], similar to the case of γ . In general, these derivatives should depend on the thermodynamic relations in the perturbation timescale [75]. The Einstein field equations also give H 2 = H 0 as in the isotropic case. The variables ( H 0 , V ) are related to the dependent variables through the algebraic relations  \n[ 3 m + ( ℓ -1)( ℓ +2) 2 r +4 πr 3 p r ] H 0 =8 πr 3 e -ν/ 2 X -[ ℓ ( ℓ +1) 2 ( m +4 πr 3 p r ) -ω 2 r 3 e -λ -ν ] H 1 + [ ( ℓ -1)( ℓ +2) 2 r -ω 2 r 3 e -ν -e λ r ( m +4 πr 3 p r )(3 m -r +4 πr 3 p r ) ] K -16 πre -λ/ 2 ( ρ + p r )¯ σW, (25)  \nX = ω 2 ( ρ + p r ) e -ν/ 2 (1 -¯ σ ) V + ρ + p r 2 e ν/ 2 H 0 -p ' r r e ( ν -λ ) / 2 W -e ν/ 2 S. (26)  \nLet us summarize how the above equations are derived. Equations (15) and (16) are the ( t, θ )-component and the ( t, r )-component of Eq. (9). Equation (17) can be obtained by combining the adiabatic relation Eq. (12) with the t -component of Eq. (10). It can also be derived with an alternative approach starting from the continuity equation as described in Appendix A. Equation (18) is obtained by eliminating H ' 0 , K ' from the r -component of Eq. (10) using the ( t, r ) and ( θ, r )-components of Eq. (9). The algebraic relation Eq. (25) is derived by eliminating the H ' 0 , K ' terms with ( t, r ), ( r, r ), and ( θ, r )-components  \nof Eq. (9), and Eq. (26) is found from the θ -component of Eq. (10). Equation (27) comes from the assumed relation  \nδσ = ∂σ ∂p r δp r + ∂σ ∂ρ δρ + ∂σ ∂µ δµ, (27)  \nwhich reduces to the Cowling limit if we take δµ = 0. To solve for the QNMs, we integrate Eqs. (15)-(18) numerically to obtain the perturbed solution of the NS interior. We begin the integration near the stellar center. Assuming the solutions are regular yields the following constraints:  \nH (0) 1 = 16 π ℓ ( ℓ +1) ( ρ (0) + p (0) r ) W (0) + 2 ℓ +1 K (0) , (28)  \nX (0) =( ρ (0) + p (0) r ) e ν (0) 2 { [ 4 π 3 ( ρ (0) +3 p (0) r ) + 2 σ (2) ρ (0) + p (0) r -ω 2 e -ν (0) ℓ ] W (0) + K (0) 2 } , (29) S (0) =0 , (30)  \nwhere the superscript ( n ) denotes the coefficient of the n th order expansion of the quantity in series of r about the stellar center. At the surface, the vanishing Lagrangian pressure perturbation gives X ( R ) = 0. The interior solution is then matched with the metric perturbation in vacuum at the stellar surface. We follow the Zerilli's method described in [72] (also check [76] for amendments of several typos) to solve for the complex frequency corresponding to purely outgoing solutions at far-field region.  \nNote that Eq. (30) combined with Eq. (27) effectively places an extra constraint on σ near the center:  \n( ∂σ ∂x ) eq = O ( r 2 ) , (31)  \nwhere x represents ( p r , ρ , µ ). This is similar to how the background structure equation (Eq. (6)) restricts σ to go as r 2 at the center (at least r to avoid irregularity). The H anisotropy model is modified from the original form in [51] to satisfy the above condition.", 'IV. UNSTABLE QUASINORMAL MODES': 'Here we show the results of the ℓ = 2 QNMs of anisotropic NS models where the p -mode frequencies can have a negative imaginary part for a certain range of anisotropic parameters. These modes are interpreted as unstable QNMs. These results are obtained by numerically solving the linearized perturbed Einstein equations and stress energy conservation as described in Sec. III. The equations have been checked to correctly reduce to those in the Cowling limit of anisotropic stars given in [64]. We performed our calculations using two independent numerical codes and found consistent results with a difference of ∼ 1%. The major code is publically available at [77].  \nIn Fig. 1, we present the imaginary parts of the QNM frequencies of the f -mode, p 1 -mode, and the p 2 -mode against β for various combinations of the ρ -p r EOS (Poly and SLy4) and the anisotropy models (H and BL). First, the imaginary part of all modes in the isotropic limit  \nPoly-H model  \nPoly-BL model  \nFIG. 1. The imaginary part of the first three ℓ = 2 QNM frequencies, f = ω/ (2 π ), against the anisotropy parameter β for different combinations of the EOSs (Poly vs SLy4) and the anisotropy models (H vs BL). The p -mode frequencies are scaled up by the constant factors specified in the legends to provide better visualization of the zero-crossings. The central densities of the Poly models and the SLy4 models are set at 7 . 455 × 10 14 g cm -3 and 9 . 88 × 10 14 g cm -3 respectively, such that the isotropic ( β = 0) NS model has a mass of 1.40 M ⊙ .  \n<!-- image -->  \n( β = 0) is also positive, which serves as a sanity check of the numerical results, as the governing equations readily reduce to the isotropic case given in [73, 76]. Next, we see that the f -modes have positive imaginary frequencies for any anisotropy parameter β . On the other hand, the imaginary parts of the p -modes cross zero as we vary β away from 0, and become negative in a certain range of β . For the p 2 -modes of the Poly-H model, it even crosses zero three times within the range of β we consider. We also found similar zero-crossing behaviors in p -modes at higher radial orders. This illustrates that there are multiple unstable p -modes within the anisotropic NSs. For a particular value of β , even the lowest order modes being stable does not necessarily imply a stable star. In principle, one needs to solve for an infinite number of p -modes to check for stability.  \nWe further illustrate the instability in contrast to the isotropic case, in which a smooth sequence of QNMs is forbidden to cross the real axis [20, 21]. We focus on the p 1 -mode frequencies and present in Fig. 2 a set of anisotropic NS models with various β in the complex frequency plane. We start increasing β from the corresponding lower bound value in the direction indicated by the arrow in each sequence of models. The frequency crosses the real axis multiple times depending on the EOSs and anisotropy models. In the case of isotropic stars, these crossings are forbidden by the equation of motion [20] as we describe in more detail below.  \nWe next perform a further analytic study by following [20] that was devised for variational principle. They showed, in the isotropic case, that the QNM frequency can be related to the eigenfunctions as an integral for-  \n<!-- image -->  \nFIG. 2. The ℓ = 2 p 1 -modes in the complex frequency plane of the models in Fig. 1 for Poly (left) and SLy4 (right) EOSs. The arrow indicates the direction of increasing β , in which the absolute values are chosen to be | β | ≤ 10 and | β | ≤ 0 . 5 for the H-model and the BL-model respectively. The value of β corresponding to the first point within the plotting range is shown for the SLy4-H model. Note that the curves intersect in the isotropic case (the intersection with Re( f ) = 6304 Hz for the right panel.  \n<!-- image -->  \nmula consisting of pairs of the perturbation variables, schematically written in the form  \nω 2 ∫ ∞ 0 I 1 dr = ∫ ∞ 0 I 2 dr + I 3 ( R ) + I 4 ( ∞ ) , (32)  \nwhere ω is the QNM frequency and I 1 to I 4 are sums of pairs of the perturbation variable eigenfunctions, which are in general complex (see Appendix B for the explicit expressions). The term I 3 is evaluated at the stellar radius R while I 4 is evaluated at infinity. In I 1 to I 3 , the perturbation variables are expressed in complex conjugate pairs such that they are guaranteed to be real. For I 4 , the outgoing wave condition for QNMs makes its imaginary part non-vanishing, except when ω is purely imaginary. As a result, ω cannot be real unless there is no outgoing wave at infinity, meaning that the QNM frequencies cannot cross the real axis in the complex frequency plane except at the origin.  \nFor the anisotropic case, we find the following extra terms on the right-hand side of Eq. (32):  \n-ω 2 ∫ R 0 J 1 dr + ∫ R 0 J 2 dr + J 3 ( R ) , (33)  \nwhere J 1 and J 3 are real, while J 2 is in general complex (see Appendix B for the explicit forms of these terms). Together with I 4 ( ∞ ) in Eq. (32), there is the possibility of ω becoming purely real. As a result, ω 2 is allowed to cross the real axis for β = 0.  \n̸', 'V. DISCUSSIONS AND CONCLUSION': 'In this work, we demonstrated that the p -modes can become unstable in some range of the anisotropy parameter β with two different anisotropy models and two NS EOSs. We analytically illustrated that such instabilities are allowed to exist based on the mathematical structure of the eigenvalue problem, in contrast to the isotropic case. In Appendix D, we compared our calculations in full GR with those within the Cowling limit and found that the latter is valid for anisotropic NSs with an error of < 10% for p -modes. Additional data are also provided for numerical comparisons.  \nThe implication of the instability still requires further analysis. We can expect the nonlinear effects will become important as the unstable mode amplitude grows exponentially with time, causing the linearized theories to eventually break down. A feedback mechanism is expected to explain the source of energy for the growth of QNMs while gravitational waves carry away the energy. This should cause some permanent changes in the stellar background and/or the anisotropy model while extracting the energy for the growth.  \nWe can draw some insights from the ChandrasekharFriedman-Schutz (CFS) instability of a rotating star [10, 12], in which a QNM is seen as a prograde mode in the inertial frame, but retrograde in the co-rotating frame, thus having opposite (canonical 1 ) angular momenta (see e.g., [79]). When it loses positive angular  \nmomentum through radiation, it is driven to even larger amplitudes since the mode carries a negative angular momentum. In this case, the energy source that drives the instability comes from the rotation of the background star. Such an effect cannot be captured without considering the feedback on the background rotation.  \nIn the setup of this study, there is no background rotation to source the growth in the mode amplitude. However, some permanent changes to the background star are expected for the mode to increase its amplitude spontaneously. One possibility is that the anisotropy will decay in a short time due to certain feedback mechanisms from the unstable mode, which can only be obtained through calculations beyond linearized perturbation theories. This may restrict the form of the anisotropy models that can be sustained in a NS.  \nThe effects of viscous damping within the NS are also neglected in this study. These effects are important in determining the suppression of the instability due to other dissipative processes. As a reference, the dominant (shear) viscosity damping rate of typical isotropic NSs is ≲ 10 -8 Hz for ℓ = 2 modes [80], which is much smaller than the growth rate of the unstable p -modes of the majority of the parameter range reported here. This means the known viscous damping effect does not suppress the instability.  \nThis work opens up several new avenues in the study of anisotropy-driven instabilities. First, to explore the physical origin of the instability, we need to identify the energy source for the growth of the mode amplitude. A variational approach may allow us to construct an energy balance law as in [20]. Second, the stability of other QNMs, like the f -modes and g -modes, is also an interesting problem to consider. Although our numerical results suggest that the former are stable within the range of β considered, it requires a thorough study with a larger range of NS parameters and anisotropy models to draw a meaningful conclusion. Finally, we can further explore the prospects of constraining the anisotropy models through the unstable QNMs by requiring stability.  \nLastly, we re-emphasize that the form of the anisotropic pressure and its perturbation in this paper are parametrized models of assumed forms instead of being derived from fundamental thermodynamic principles, which can be unrealistic and are considered effective models only. The recent series of preprints [81-83] address the issues of these existing anisotropic models that enforce an explicit r -dependence to fix the problem of singularity at the stellar center. They point out that these anisotropic EOSs do not inform how matter perturbs and can violate the principle of equivalence. For this reason, they developed the theory of a relativistic anisotropic fluid starting from the action using a consistent model inspired by the theory of liquid crystals. This theory can be applied to non-radially symmetric configurations and allows a self-consistent derivation of the pulsation equations. Similar analyses of anisotropy from elasticity in relativity can be found in [27, 50]. As  \nwe discuss in Appendix A, our formalism can be generalized by modifying the form of the perturbed anisotropy to describe anisotropy from different origins (e.g., elasticity, magnetic field, etc). It is worth exploring whether the instability reported in this paper exists for QNMs derived from such a theory, i.e., whether it emerges from some flawed assumptions in the anisotropy model or physically exists.  \nNote added.While we were working on this project, a preprint [84] appeared on arXiv that provided a derivation of perturbation equations for oscillating anisotropic NSs in full GR. However, we noticed that the authors have incorrectly assumed an isotropic deformation of the fluid when applying the thermodynamic relations, while the deformation along each direction should be different for a generic nonradial deformation (see Eqs. (A3) and (A4) of Appendix A). Hence, the perturbation equations presented in [84] do not correctly reduce to those in the Cowling limit [64] and are inconsistent with those in this paper that were derived through a different approach than [84]. Furthermore, Ref. [84] used the original anisotropy model by Horvat et al . [51] that does not satisfy the correct boundary conditions for the perturbation equations. The main focus of this paper is on the stability of QNMs, which is also different from [84].', 'ACKNOWLEDGMENTS': 'K.Y. acknowledges support from NSF Grant PHY2339969 and the Owens Family Foundation. The authors thank Phil Arras for useful discussions. The authors also thank the anonymous referee for providing important comments on the origin of the instability and pointing out interesting references about the CFS instability.', 'Appendix A: An alternative derivation with the continuity equation': "In this section, we provide an alternative way to derive the pulsation equations, Eq. (17) in specific. Instead of using the t -component of Eq. (10), we utilize the continuity equation of particle number density and the relation of thermodynamics. At the end, we explain how this method allows incorporating a more general form of the EOS in the perturbed configuration.  \nWe begin by reviewing the derivation of the pulsation equations for isotropic stars. It is common to introduce the particle number density, n , which is governed by the continuity equation (see, e.g., [71]):  \n∇ α ( nu α ) = 0 = ⇒ ∆ n n = -1 2 h αβ ∆ g αβ . (A1)  \nThe number density is an independent variable in the EOS, i.e., p = p ( n ), ρ = ρ ( n ), but does not enter the equation of motion directly once we have ρ as a function of p . The pressure p here is the isotropic pressure.  \nThe Lagrangian perturbation of n is related to that of ρ through the thermodynamic relation  \n∆ ρ = ( ρ + p ) ∆ n n . (A2)  \nCombining Eqs. (A1), (A2) and the adiabatic relation γ ∆ ρ/ ( ρ + p ) = ∆ p/p gives the t -component of the stress energy conservation (Eq. (10)) of an isotropic fluid.  \nWhile Eq. (A1) still holds in the anisotropic case, Eq. (A2) needs to be modified, since the work done on the fluid does not depend only on the volume change, but is directional dependent 2 . This requires us to explicitly write down the change in energy due to work done in each direction (see, e.g., [50, 85], for the generalized thermodynamic relations in an anisotropic medium):  \n∆ ρ = ( ρ + p r ) ∆ n n + σ (∆ U θ θ +∆ U ϕ ϕ ) , (A3)  \nwhere ∆ U αβ is the perturbative Lagrangian strain tensor [86] 3 , given by  \n∆ U αβ = ∆ h αβ = 1 2 h α µ h β ν ∆ g µν , (A4)  \nand ∆ U α β = g αµ ∆ U µβ . Equation (A3), substituted with the adiabatic relation, Eq. (12), reads  \n∆ p r p r = γ [ ∆ n n + ¯ σ (∆ U θ θ +∆ U ϕ ϕ ) ] . (A5)  \nOne can show that this is equivalent to Eq. (17) using Eqs. (A1) and (A4), which are explicitly written as  \n∆ n n = ∑ ℓ,m [ 1 2 H 0 + K -e -λ/ 2 r W '  \n-ℓ +1 r 2 e -λ/ 2 W -ℓ ( ℓ +1) r 2 V ] r ℓ Y ℓm e iωt , (A6)  \n∆ U θ θ +∆ U ϕ ϕ = ∑ ℓ,m [ -K +2 e -λ/ 2 W r 2 + ℓ ( ℓ +1) r 2 V ] r ℓ Y ℓm e iωt . (A7)  \nNotice that, in general, the anisotropic fluid can have different proportionality constants in each direction depending on the anisotropy EOS. That is,  \n∆ p j = -∑ γ ij ∆ U i i , (A8)  \np j i where γ ij denotes the proportionality constant in the stress-strain relation for the spatial indices i, j . This suggests a further generalization of the pulsation equations described in the formulation section. An even more general form can include non-zero contributions from offdiagonal elements of ∆ U i j like in the elastic materials. We leave such studies for future work.", 'Appendix B: Integral relation of the eigenvalue and eigenfunctions': "The eigenfrequencies and the eigenfunctions of the pulsation equations are related through the integral equations Eqs. (32) and (33), where J 1 to J 3 in Eq. (33) vanish in the isotropic limit and are non-zero only inside the star. In this section, we provide the explicit forms of the terms I 1 to I 4 and J 1 to J 3 .  \nThe terms I 1 to I 4 are obtained from Eq. (16) of [20], and are explicitly written as:  \nI 1 = e λ -ν 2 r 2 { ( ρ + p r ) [ 1 r 4 ∣ ∣ ∣ ˜ W ∣ ∣ ∣ 2 + ℓ ( ℓ +1) r 2 ∣ ∣ ∣ ˜ V ∣ ∣ ∣ 2 ] -ℓ ( ℓ +1) 16 πe λ r 2 ∣ ∣ ∣ ˜ H 1 ∣ ∣ ∣ 2 -1 16 π [ ∣ ∣ ∣ ˜ K ∣ ∣ ∣ 2 +2Re ( ˜ K ˜ H ∗ 0 ) ]} , (B1) I 2 = e λ + ν 2 r 2 { -( ρ + p r ) ν ' A 2 e λ r 4 ∣ ∣ ∣ ˜ W ∣ ∣ ∣ 2 + 1 γp r ∣ ∣ ∣ ˜ P ∣ ∣ ∣ 2 + [ ℓ ( ℓ +1) 16 πr 2 -3( ρ + p r ) 4 ] ∣ ∣ ∣ ˜ H 0 ∣ ∣ ∣ 2 -1 16 πe λ ∣ ∣ ∣ ˜ K ' ∣ ∣ ∣ 2 + ρ + p r γp r Re ( ˜ P ˜ H ∗ 0 ) -( ρ + p r ) A e λ/ 2 r 2 Re( ˜ W ˜ H ∗ 0 ) + 1 8 πe λ Re [( ν ' ˜ H 0 + ˜ H ' 0 ) ˜ K ' ∗ ] } , (B2) I 3 = e ν ρ 2 [ 2 e -ν/ 2 Re ( ˜ W ˜ H ∗ 0 ) + ν ' r 2 ∣ ∣ ∣ ˜ W ∣ ∣ ∣ 2 ] , (B3) (B4)  \nI 4 = e ν r 2 16 π [ ˜ K ∗ ( ˜ H ' 0 + ν ' ˜ H 0 -˜ K ' ) -˜ K ' ˜ H ∗ 0 ] .  \nHere, ˜ P is defined through δp r = ∑ ℓ,m ˜ PY ℓm e iωt , and the symbol Re denotes taking the real part. The perturbation variables are redefined in the following ways to simplify the expressions:  \nr ℓ K = ˜ K, r ℓ H 0 = ˜ H 0 , r ℓ V = -˜ V , r ℓ S = ˜ S, r ℓ +1 H 1 = ˜ H 1 , r ℓ +1 W = ˜ W. (B5)  \n2  \nWe see that I 1 to I 3 are real while I 4 is in general complex if we impose the purely outgoing boundary condition at infinity [20, 87].  \nThe terms J 1 to J 3 are derived using the procedures described in Appendix A of [20], with the Einstein equations and stress-energy conservation for the anisotropic case:  \nJ 1 = -ℓ ( ℓ +1) e λ -ν 2 ( ρ + p r )¯ σ ∣ ∣ ∣ ˜ V ∣ ∣ ∣ , (B6) J 2 = e ν 2 [ ( ρ + p r )¯ σ ( -˜ K ' + 2 e -λ 2 A r 3 ˜ W ) -2¯ σ r ( 1 + ρ + p r γp r ) ˜ P + 2 r ˜ S ] ˜ W ∗ + e ν + λ 2 ( ¯ σ ˜ P + ρ + p r 2 ¯ σ ˜ H 0 -˜ S ) [ r 2 ˜ K ∗ + ℓ ( ℓ +1) ˜ V ∗ ] , (B7) J 3 = 2 e ν ρ 3 ¯ σ ∣ ∣ ∣ ˜ W ∣ ∣ ∣ 2 . (B8)  \nr  \nThe perturbation variables of J 1 and J 3 are paired in the same way as I 1 , while J 2 is in general complex as the complex perturbation variables cannot be paired in the same ways as I 1 to I 3 . Although it is not obvious that the complex pulsation variables within J 2 cannot be paired up to form real terms, we verified this by numerically computing the integral using the solutions of the QNMs.  \nAnother perspective is to focus on ˜ S , which depends on the anisotropy ansatz itself. The choice of the functional form for σ affects the perturbation variable ˜ S , which in general causes J 2 to be complex. In [48], it is shown that the non-radial modes are stable as long as the anisotropy is finite in the Newtonian incompressible star with ˜ S = β ˜ P . We point out here that it holds specifically for this anisotropy ansatz and therefore does not contradict our findings in this paper. In addition, the instability comes from the imaginary part of the eigenfrequency, which is caused by the GR effect.", 'Appendix C: Eigenfrequencies in Cowling limit': "Here, we show that the instability reported in this paper does not occur under Cowling approximation by illustrating that the perturbation variables in J 2 (see Eq. (B7)) all pair up to be real when we drop the metric perturbations.  \nThe full GR formalism can be reduced to the relativistic Cowling limit given in [64]. Notice that the approximation procedure is not simply taking all metric perturbations to zero. The steps involved are the same as the isotropic case described in, e.g., [68, 88], where we set H 0 , K , H 1 to zero in Eqs. (17), (18), and (26), and further drop the term -4 π ( ρ + p r ) 2 e ( λ + ν ) / 2 W/r in Eq. (18) 4 . The resulting formulas are equivalent to Eqs. (27), (29) and (30) of [64] with a barotropic EOS, and can further be reduced to the isotropic limit [90].  \nWriting the Cowling approximation pulsation equations in terms of the variables with tilde introduced in Eq. (B5), we have  \n˜ W ' = ( 2¯ σ r -p ' r γp r ) ˜ W -re λ/ 2 γp r ˜ P -ℓ ( ℓ +1) e λ/ 2 (1 -¯ σ ) ˜ V , (C1)  \n˜ P ' = ρ + p r r 2 e λ/ 2 [ ω 2 e -ν + ν ' 2 e -λ A ] ˜ W -( 1 + ρ + p r γp r ) ν ' 2 ˜ P -2 r ˜ S, (C2) ˜ P = -( ρ + p r )(1 -¯ σ ) ω 2 e -ν ˜ V + ˜ S. (C3)  \nSubstituting Eqs. (C1)-(C3) into Eq. (B7) and setting all metric perturbations to zero, we can obtain  \nJ 2 = -( e ν 2 ˜ P ˜ W ∗ ) ' + ω 2 e λ -ν 2 ( ρ + p r ) [ 1 r 2 ∣ ∣ ∣ ˜ W ∣ ∣ ∣ 2 -(1 -¯ σ ) ℓ ( ℓ +1) ∣ ∣ ∣ ˜ V ∣ ∣ ∣ 2 ] + e ν -λ 2 ( ρ + p r ) ν ' 2 A r 2 ∣ ∣ ∣ ˜ W ∣ ∣ ∣ 2 , (C4)  \nway as H 0 and K . An argument about H 1 being larger than the other metric perturbations in the Newtonian limit is given in [89]  \nwhere the first term becomes a real boundary term after integration. Hence, combining I 1 to I 3 and J 1 to J 3 under Cowling approximation, ω 2 must be real.", 'Appendix D: Additional numerical data': "In this section, we provide additional numerical data of the QNMs computed for the EOSs and anisotropy models shown in the main text. We demonstrate that the anisotropy has a small effect on the real part of the frequency, even though it changes the stability of the mode. We also quantify the validity of the Cowling approximation for computing the oscillation modes of anisotropic stars for the first time.  \nIn Fig. 3, we show the real and imaginary parts of the QNM frequencies for the first six ℓ = 2 modes of the Poly-H model. The real part of the frequencies of the isotropic and the anisotropic case is of similar values  \n- [1] D. Tsang, J. S. Read, T. Hinderer, A. L. Piro, and R. Bondarescu, Resonant shattering of neutron star crusts, Phys. Rev. Lett. 108 , 011102 (2012).\n- [2] C. Chirenti, R. Gold, and M. C. Miller, Gravitational waves from f-modes excited by the inspiral of highly eccentric neutron star binaries, Astrophys. J. 837 , 67 (2017).\n- [3] H. Yu and N. N. Weinberg, Dynamical tides in coalescing superfluid neutron star binaries with hyperon cores and their detectability with third-generation gravitationalwave detectors, Mon. Notices Royal Astron. Soc. 470 , 350-360 (2017).\n- [4] G. Pratten, P. Schmidt, and T. Hinderer, Gravitationalwave asteroseismology with fundamental modes from compact binary inspirals, Nat. Commun. 11 , 10.1038/s41467-020-15984-5 (2020).\n- [5] Z. Pan, Z. Lyu, B. Bonga, N. Ortiz, and H. Yang, Probing crust meltdown in inspiraling binary neutron stars, Phys. Rev. Lett. 125 , 201102 (2020).\n- [6] H. H.-Y. Ng, P. C.-K. Cheong, L.-M. Lin, and T. G. F. Li, Gravitational-wave asteroseismology with f-modes from neutron star binaries at the merger phase, Astrophys. J. 915 , 108 (2021).\n- [7] S. Y. Lau and K. Yagi, Probing hybrid stars with gravitational waves via interfacial modes, Phys. Rev. D 103 , 063015 (2021).\n- [8] A. Passamonti, N. Andersson, and P. Pnigouras, Dynamical tides in neutron stars: the impact of the crust, Mon. Notices Royal Astron. Soc. 504 , 1273-1293 (2021).\n- [9] J. P. Cox, Non-adiabatic stellar pulsation. i., Astrophys. J. 127 , 194 (1958).\n- [10] S. Chandrasekhar, Solutions of two problems in the theory of gravitational radiation, Phys. Rev. Lett. 24 , 611 (1970).\n- [11] P. R. Wood, Pulsation of high luminosity helium stars, Mon. Not. R. Astron. Soc. 174 , 531-539 (1976).\n- [12] J. L. Friedman and B. F. Schutz, Secular instability of rotating newtonian stars, Astrophys. J. 222 , 281 (1978).  \nfor all the modes. For the isotropic case, the imaginary part of the modes is all above zero and decreases with the mode order. Meanwhile, that of the anisotropic case does not show a monotonic behavior and can become negative.  \nIn Fig. 3, we also show the frequencies obtained by the Cowling approximation. The f -mode frequencies show 20-30% deviations from the full GR real frequencies. The percentage deviation decreases if we go to higher radial order p -modes. In the case of the p 5 -modes, the deviation reduces to a 1% level.  \nMore numerical data are provided in Table. I, including the real and imaginary frequencies of the first three ℓ = 2 QNMs, as a reference for readers interested in a numerical comparison. The β = 0 case of the SLy4 EOS has been checked against numerical data available (e.g., [91]), and attains good agreement. We see that with the same ρ -p r relation, the real part of the frequencies is not very sensitive to the change in β , while the imaginary part can change by orders of magnitude and even change the sign.  \n- [13] J. P. Cox, Theory of stellar pulsation (Princeton University Press, Princeton, NJ, 1980).\n- [14] A. Gautschy, W. Glatzel, A. Gautschy, and W. Glatzel, On highly non-adiabatic stellar pulsations and the origin of strange modes, Mon. Not. R. Astron. Soc. 245 , 597-597 (1990).\n- [15] T. G. Cowling, The non-radial oscillations of polytropic stars, Mon. Not. R. Astron. Soc. 101 , 367-375 (1941).\n- [16] M. Gabriel and R. Scuflaire, Properties of non-radial stellar oscillations., Acta Astron. 29 , 135 (1979).\n- [17] W. Unno, Y. Osaki, H. Ando, and H. Shibahashi, Nonradial Oscillation of Stars (University of Tokyo Press, Tokyo, Japan, 1979).\n- [18] K. D. Kokkotas and B. G. Schmidt, Quasi-normal modes of stars and black holes, Living Rev. Relativ. 2 , 10.12942/lrr-1999-2 (1999).\n- [19] K. S. Thorne and A. Campolattaro, Non-radial pulsation of general-relativistic stellar models. i. analytic analysis for l ≥ 2, Astrophys. J. 149 , 591 (1967).\n- [20] S. L. Detweiler and J. R. Ipser, A variational principle and a stability criterion for the non-radial modes of pulsation of stellar models in general relativity, Astrophys. J. 185 , 685 (1973).\n- [21] J. R. Ipser, On the implications of the nonexistence of unstable outgoing normal modes for the stability of spherical stellar models to nonradial perturbations in general relativity, Astrophys. J. 199 , 220 (1975).\n- [22] N. Lebovitz, W. Reid, and P. Vandervoort, Theoretical principles in astrophysics and relativity (University of Chicago Press, Chicago, IL, 1981).\n- [23] K. S. Wong, L. M. Lin, and P. T. Leung, Universality in oscillation modes of superfluid neutron stars?, Astrophys. J. 699 , 1809-1821 (2009).\n- [24] T. Harada, Stability analysis of spherically symmetric star in scalar-tensor theories of gravity, Prog. Theor. Phys. 98 , 359-379 (1997).\n- [25] R. F. P. Mendes and N. Ortiz, New class of quasinormal modes of neutron stars in scalar-tensor gravity, Phys.  \nFIG. 3. The real (top panel) and imaginary (bottom panel) frequencies of the first six ℓ = 2 QNMs of the Poly-H model with central density 7.455 × 10 14 g cm -3 . The left panel has β = 0 and the right panel has β = 10. We also present the normal mode frequencies computed using the relativistic Cowling approximation (Cowling) described in [64].  \n<!-- image -->  \n𝑓  \n𝑝  \n1  \nPoly-H (  \n𝛽  \n= 10)  \n𝑝  \n2  \n𝑝  \n3  \n𝑝  \n4  \n𝑝  \n5  \nTABLE I. Numerical data of the first three ℓ = 2 QNM frequencies of the anisotropic NS models. The central densities of the Poly models and the SLy4 models are set at 7 . 455 × 10 14 g cm -3 and 9 . 88 × 10 14 g cm -3 respectively.  \nRev. Lett. 120 , 201104 (2018).  \n[26] H. Heintzmann and W. Hillebrandt, Neutron stars with an anisotropic equation of state: mass, redshift and stability., Astron. Astrophys. 38 , 51 (1975).  \n[27] M. Karlovini and L. Samuelsson, Elastic stars in general  \nrelativity: I. foundations and equilibrium models, Class. Quantum Gravity 20 , 3613-3648 (2003).  \n- [28] S. Nelmes and B. M. A. G. Piette, Phase transition and anisotropic deformations of neutron star matter, Phys. Rev. D 85 , 123004 (2012).  \n)  \nz  \nH  \n(  \n)  \n𝑓  \n(  \ne  \nR  \n)  \nz  \nH  \n(  \n)  \n𝑓  \n(  \nm  \nI  \n1×10  \n5×10  \n1×10  \n4  \n3  \n3  \n0.5  \n0  \n-0.5  \n- [29] P. S. Letelier, Anisotropic fluids with two-perfect-fluid components, Phys. Rev. D 22 , 807 (1980).\n- [30] L. Herrera, Local anisotropy in self-gravitating systems, Phys. Rep. 286 , 53-130 (1997).\n- [31] R. F. Sawyer and D. J. Scalapino, Pion condensation in superdense nuclear matter, Phys. Rev. D 7 , 953 (1973).\n- [32] C. Y. Cardall, M. Prakash, and J. M. Lattimer, Effects of strong magnetic fields on neutron star structure, Astrophys. J. 554 , 322-339 (2001).\n- [33] X.-G. Huang, M. Huang, D. H. Rischke, and A. Sedrakian, Anisotropic hydrodynamics, bulk viscosities, and r -modes of strange quark stars with strong magnetic fields, Phys. Rev. D 81 , 045015 (2010).\n- [34] E. J. Ferrer, V. de la Incera, J. P. Keith, I. Portillo, and P. L. Springsteen, Equation of state of a dense and magnetized fermion system, Phys. Rev. C 82 , 065802 (2010).\n- [35] S. S. Yazadjiev, Relativistic models of magnetars: Nonperturbative analytical approach, Phys. Rev. D 85 , 044030 (2012).\n- [36] M. Sinha, B. Mukhopadhyay, and A. Sedrakian, Hypernuclear matter in strong magnetic field, Nucl. Phys. A 898 , 43-58 (2013).\n- [37] W. Barreto and S. Rojas, An equation of state for radiating dissipative spheres in general relativity, Astrophys. Space Sci. 193 , 201-215 (1992).\n- [38] D. J. Kaup, Klein-gordon geon, Phys. Rev. 172 , 1331 (1968).\n- [39] F. E. Schunck and E. W. Mielke, General relativistic boson stars, Class. Quantum Gravity 20 , R301-R356 (2003).\n- [40] M. Alcubierre, J. Barranco, A. Bernal, J. C. Degollado, A. Diez-Tejedor, V. Jaramillo, M. Megevand, D. N'u˜nez, and O. Sarbach, Extreme ℓ -boson stars, Class. Quantum Gravity 39 , 094001 (2022).\n- [41] D. Deb, S. R. Chowdhury, S. Ray, F. Rahaman, and B. Guha, Relativistic model for anisotropic strange stars, Ann. Phys. (N. Y.) 387 , 239-252 (2017).\n- [42] A. Errehymy, M. Daoud, and E. H. Sayouty, A spherically symmetric model of anisotropic fluid for strange quark spheres, Eur. Phys. J. C 79 , 10.1140/epjc/s10052-0196862-9 (2019).\n- [43] P. H. R. S. Moraes, G. Panotopoulos, and I. Lopes, Anisotropic dark matter stars, Phys. Rev. D 103 , 084023 (2021).\n- [44] C. Cattoen, T. Faber, and M. Visser, Gravastars must have anisotropic pressures, Class. Quantum Gravity 22 , 4189-4202 (2005).\n- [45] C. B. M. H. Chirenti and L. Rezzolla, How to tell a gravastar from a black hole, Class. Quantum Gravity 24 , 4191-4206 (2007).\n- [46] R. L. Bowers and E. P. T. Liang, Anisotropic spheres in general relativity, Astrophys. J. 188 , 657 (1974).\n- [47] S. m. c. i. m. c. Bayin, Anisotropic fluid spheres in general relativity, Phys. Rev. D 26 , 1262 (1982).\n- [48] W. Hillebrandt and K. O. Steinmetz, Anisotropic neutron star models: stability against radial and nonradial pulsations., Astron. Astrophys. 53 , 283 (1976).\n- [49] K. Dev and M. Gleiser, Anisotropic stars II: Stability, Gen. Relativ. Gravit. 35 , 1435 (2003).\n- [50] M. Karlovini, L. Samuelsson, and M. Zarroug, Elastic stars in general relativity: Ii. radial perturbations, Class. Quantum Gravity 21 , 1559-1581 (2004).\n- [51] D. Horvat, S. Iliji'c, and A. Marunovi'c, Radial pulsations and stability of anisotropic stars with a quasi-local equa-  \n- on of state, Class. Quantum Gravity 28 , 025009 (2010).\n- [52] D. Horvat, S. Iliji'c, and A. Marunovi'c, Radial stability analysis of the continuous pressure gravastar, Class. Quantum Gravity 28 , 195008 (2011).\n- [53] J. D. Arba˜nil and M. Malheiro, Radial stability of anisotropic strange quark stars, J. Cosmol. Astropart. Phys. 2016 (11), 012-012.\n- [54] A. A. Isayev, General relativistic polytropes in anisotropic stars, Phys. Rev. D 96 , 083007 (2017).\n- [55] J. M. Z. Pretel, Equilibrium, radial stability and nonadiabatic gravitational collapse of anisotropic neutron stars, Eur. Phys. J. C 80 , 10.1140/epjc/s10052-020-83013 (2020).\n- [56] J. M. Bardeen, K. S. Thorne, and D. W. Meltzer, A catalogue of methods for studying the normal modes of radial pulsation of general-relativistic stellar models, Astrophys. J. 145 , 505 (1966).\n- [57] H. O. Silva, C. F. B. Macedo, E. Berti, and L. C. B. Crispino, Slowly rotating anisotropic neutron stars in general relativity and scalar-tensor theory, Class. and Quantum Gravity 32 , 145008 (2015).\n- [58] K. Yagi and N. Yunes, I-love-q anisotropically: Universal relations for compact stars with scalar pressure anisotropy, Phys. Rev. D 91 , 123008 (2015).\n- [59] K. Yagi and N. Yunes, I-Love-Q, Science 341 , 365 (2013), arXiv:1302.4499 [gr-qc].\n- [60] K. Yagi and N. Yunes, I-Love-Q Relations in Neutron Stars and their Applications to Astrophysics, Gravitational Waves and Fundamental Physics, Phys. Rev. D 88 , 023009 (2013), arXiv:1303.1528 [gr-qc].\n- [61] K. Yagi and N. Yunes, Approximate Universal Relations for Neutron Stars and Quark Stars, Phys. Rept. 681 , 1 (2017), arXiv:1608.02582 [gr-qc].\n- [62] K. Yagi and N. Yunes, Relating follicly-challenged compact stars to bald black holes: A link between two no-hair properties, Phys. Rev. D 91 , 103003 (2015), arXiv:1502.04131 [gr-qc].\n- [63] K. Yagi and N. Yunes, I-love-q relations: from compact stars to black holes, Class. Quantum Gravity 33 , 095005 (2016).\n- [64] D. D. Doneva and S. S. Yazadjiev, Nonradial oscillations of anisotropic neutron stars in the cowling approximation, Phys. Rev. D 85 , 124023 (2012).\n- [65] E. J. A. Curi, L. B. Castro, C. V. Flores, and C. H. Lenzi, Non-radial oscillations and global stellar properties of anisotropic compact stars using realistic equations of state, Eur. Phys. J. C 82 , 527 (2022), arXiv:2206.09260 [gr-qc].\n- [66] S. R. Mohanty, S. Ghosh, P. Routaray, H. C. Das, and B. Kumar, The impact of anisotropy on neutron star properties: Insights from i-f-c universal relations (2023).\n- [67] O. P. Jyothilakshmi, L. J. Naik, and V. Sreekanth, Nonradial oscillations in anisotropic dark energy stars (2024).\n- [68] L. Lindblom and R. J. Splinter, The Accuracy of the Relativistic Cowling Approximation, Astrophys. J. 348 , 198 (1990).\n- [69] F. Douchin and P. Haensel, A unified equation of state of dense matter and neutron star structure, Astron. Astrophys. 380 , 151-167 (2001).\n- [70] L. M. Becerra, E. A. Becerra-Vergara, and F. D. LoraClavijo, Realistic anisotropic neutron stars: Pressure effects (2024).\n- [71] J. L. Friedman and B. F. Schutz, On the stability of relativistic systems, Astrophys. J. 200 , 204 (1975).\n- [72] L. Lindblom and S. L. Detweiler, The quadrupole oscillations of neutron stars, Astrophys. J., Suppl. Ser. 53 , 73 (1983).\n- [73] S. Detweiler and L. Lindblom, On the nonradial pulsations of general relativistic stellar models, Astrophys. J. 292 , 12 (1985).\n- [74] T. Regge and J. A. Wheeler, Stability of a schwarzschild singularity, Phys. Rev. 108 , 1063 (1957).\n- [75] K. S. Thorne, The General Relativistic Theory of Stellar Structure and Dynamics, in High Energy Astrophysics, Volume 3 , Vol. 3, edited by C. Dewitt, E. Schatzman, and P. V'eron (1967) p. V.\n- [76] J.-L. Lu and W.-M. Suen, Determining the long living quasi-normal modes of relativistic stars, Chin. Phys. B 20 , 040401 (2011).\n- [77] The code for computing the non-radial modes of anisotropic neutron stars in Fortran is publically available on GitHub at https://github.com/SYLau/code-ansqnm-share.git.\n- [78] Y. Levin and G. Ushomirsky, Non-linear r-modes in a spherical shell: issues of principle, Mon. Notices Royal Astron. Soc. 322 , 515-522 (2001).\n- [79] M. Maggiore, Gravitational Waves. Vol. 2: Astrophysics and Cosmology (Oxford University Press, London, England, 2018).\n- [80] C. Cutler, L. Lindblom, and R. J. Splinter, Damping times for neutron star oscillations, Astrophys. J. 363 , 603 (1990).\n- [81] T. Cadogan and E. Poisson, Self-gravitating anisotropic  \n- fluids. i: Context and overview (2024).\n- [82] T. Cadogan and E. Poisson, Self-gravitating anisotropic fluid. ii: Newtonian theory (2024).\n- [83] T. Cadogan and E. Poisson, Self-gravitating anisotropic fluid. iii: Relativistic theory (2024).\n- [84] S. Mondal and M. Bagchi, Non-radial oscillation of anisotropic neutron stars in full general relativity (2023).\n- [85] L. Landau and E. Lifshits, Theory of Elasticity (Pergamon, 1959).\n- [86] B. Carter and H. Quintana, Proc. R. Soc. A 331 , 57-83 (1972).\n- [87] R. Price and K. S. Thorne, Non-radial pulsation of general-relativistic stellar models. ii. properties of the gravitational waves, Astrophys. J. 155 , 163 (1969).\n- [88] T. Zhao, C. Constantinou, P. Jaikumar, and M. Prakash, Quasinormal g modes of neutron stars with quarks, Phys. Rev. D 105 , 103025 (2022).\n- [89] L. S. Finn, Relativistic stellar pulsations in the Cowling approximation, Mon. Not. R. Astron. Soc. 232 , 259 (1988), https://academic.oup.com/mnras/articlepdf/232/2/259/3316314/mnras232-0259.pdf.\n- [90] P. N. McDermott, H. M. van Horn, and J. F. Scholl, Nonradial g-mode oscillations of warm neutron stars, Astrophys. J. 268 , 837 (1983).\n- [91] A. Kunjipurayil, T. Zhao, B. Kumar, B. K. Agrawal, and M. Prakash, Impact of the equation of state on f - and p - mode oscillations of neutron stars, Phys. Rev. D 106 , 063005 (2022)."}