Equilibrium, radial stability and non-adiabatic gravitational collapse of anisotropic neutron stars

In this work we construct families of anisotropic neutron stars for an equation of state compatible with the constraints of the gravitational-wave event GW170817 and for four anisotropy ansatze. Such stars are subjected to a radial perturbation in order to study their stability against radial oscillations and we develop a dynamical model to describe the non-adiabatic gravitational collapse of the unstable anisotropic configurations whose ultimate fate is the formation of a black hole. We find that the standard criterion for radial stability $dM/d\rho_c>0$ is not always compatible with the calculation of the oscillation frequencies for some anisotropy ansatze, and each anisotropy parameter is constrained taking into account the recent restriction of maximum mass of neutron stars. We further generalize the TOV equations within a non-adiabatic context and we investigate the dynamical behaviour of the equation of state, heat flux, anisotropy factor and mass function as an unstable anisotropic star collapses. After obtaining the evolution equations we recover, as a static limit, the background equations.


I. INTRODUCTION
The most common matter-energy distribution for modeling the internal structure of compact stars is an isotropic perfect fluid. Nevertheless, there are strong arguments suggesting that nuclear matter at very high densities and pressures could naturally be described by an anisotropic fluid, that is, when the radial and tangential components of the pressure are not equal. As a matter of fact, the anisotropy could be generated by the presence of strong magnetic fields [1][2][3][4], solid cores [5,6], superfluidity [7,8], pion condensed phase configurations in neutron stars [9,10], etc. Furthermore, it is possible to obtain an anisotropic perfect fluid by combining the energy-momentum tensors of two isotropic perfect fluids [11][12][13].
As widely reported in the literature [14][15][16][17][18][19][20][21][22][23][24][25], the presence of anisotropy affects a number of important physical properties of compact stars such as the mass-radius relation, compactness, surface redshift, moment of inertia as well as the scalarization in scalar-tensor theories [26]. Indeed, the equation of state (EoS) plays a fundamental role in determining the internal structure of such stars and, consequently, in imposing stability limits. Therefore, it is important to carry out a stability analysis of anisotropic neutron stars taking into account the LIGO-Virgo constraints on the EoS for nuclear matter as a result of observation of the event GW170817 -the first direct detection of gravitational waves from the coalescence of a neutron star binary system [27,28]. In that respect, we are interested in considering a realistic EoS, which is compatible with the restriction obtained from this merger, and exploring the effects it can have on the physical characteristics of stable and unstable anisotropic * juanzarate@if.ufrj.br stars.
It is a well known fact that the solutions of the Tolman-Oppenheimer-Volkov (TOV) equations describe stellar configurations in hydrostatic equilibrium. Nonetheless, such equilibrium can be stable or unstable with respect to a compression or decompression caused by radial perturbations. Indeed, a conventional technique widely used to indicate the onset of instability is the M (ρ c ) method [17], also known as a necessary condition for stability analysis of compact stars [29,30]. This boundary between the stable and unstable stars describes the maximum amount of mass that can exist in a configuration before it must undergo a gravitational collapse. On the other hand, a sufficient condition for stability is to calculate the frequencies of the normal radial modes of relativistic vibrations [29][30][31], where Einstein field equations have to be linearized around the equilibrium configuration. If any of these squared frequencies is positive we have stable radial oscillations, whereas negative squared frequencies imply increasing or decreasing perturbations with time, i.e. the star is unstable.
In general relativity (GR), the stability analysis for isotropic compact stars with respect to radial perturbations has been widely discussed in the literature [32][33][34][35][36][37][38][39][40][41][42][43][44][45][46], whereas within the context of anisotropic configurations the normal radial modes technique has been used only in some specific cases, see e.g., [17,20,[47][48][49][50]. In particular, for anisotropic strange stars described by the MIT bag model EoS, it was shown that the M (ρ c ) method is not compatible with the calculation of frequencies to predict the onset of instability [20] for the anisotropy profile proposed by Bowers and Liang [14], this is, the maximummass stellar configurations do not correspond to the zero squared frequencies of the fundamental mode. Therefore, it is required to calculate the frequency of the oscillation modes in order to have absolute certainty about the radial stability of an anisotropic neutron star. But what could be the origin of the stellar oscillations? Just as the oscillations inside the Earth are excited by earthquakes and used in seismology to study the structure of the Earth, the fluid pulsations in neutron stars can be excited by cracks in the crust (a starquake) [51], as well as by tidal effects in binary inspirals [52,53].
In the present paper, in addition to generating hydrostatically stable anisotropic configurations, we are also interested in studying the process of gravitational collapse of the unstable configurations. In this regard, the pioneering work about gravitational collapse for a spherically symmetric distribution of matter in the form of dust cloud was carried out by Oppenheimer and Snyder [54]. Such idealized treatment has been improved by introducing a pressure gradient [55] and by replacing the exterior Schwarzschild metric by the Vaidya one [56]. In this way, the relativistic equations for the adiabatic, spherically symmetric gravitational collapse as given in Ref. [55] were modified by Misner [57] in order to allow an extremely simplified heat-transfer process (where the internal energy is converted into an outward flux of neutrinos). Moreover, the equations that govern the gravitational collapse of a ball of charged perfect fluid were derived by Bekenstein [58].
Over the years, it has been possible to construct physically viable gravitational collapse scenarios for isotropic fluids that include dissipative fluxes such as heat flow [59][60][61][62][63][64] and shear and bulk viscosities [65]. In turn, gravitational collapse models have also been developed for anisotropic fluids with dissipative processes [66][67][68][69][70][71], in the presence of electromagnetic field [72][73][74][75][76] and with cosmological constant [77,78]. Particularly the gravitational collapse of neutron stars considered as initial configurations has been investigated by the authors in Refs. [66,70,79,80], however, they did not perform an a priori analysis on the stellar stability of the initial static Schwarzschild configurations against radial pulsations. A complete analysis on stability and gravitational collapse for isotropic fluids was conducted by Ghezzi [81] (where the charged neutron stars in the unstable branch collapse directly to form black holes), and more recently in Ref. [45] for neutron stars with realistic EoSs.
We construct families of anisotropic neutron stars based on an EoS compatible with the recent observations, and we obtain the oscillation spectrum for each family by means of radial perturbations. In addition, we investigate the dynamical evolution of unstable anisotropic stars whose final fate is the formation of a black hole as a consequence of a non-adiabatic gravitational collapse. The spherical surface of the collapsing star divides spacetime into two different four-dimensional manifolds; an interior region with heat flux -described by a shear-free line element -and an exterior region which is described by the Vaidya metric for pure outgoing radiation. The paper is organized as follows: In Sec. II we present the basic formalism to describe equilibrium configurations for four anisotropy ansatze and deal with the normal radial modes. In Sec. III we develop a gravitational collapse model by introducing a time dependency on the metric functions that allows us to recover the static case under a certain limit. In Sec. IV we present a discussion of the numerical results. The paper ends with our conclusions in Sec. V. We adopt the signature (−, +, +, +) and physical units will be used throughout this work.

II. BASIC EQUATIONS FOR STELLAR STRUCTURE
The line element describing the interior spacetime of a spherically symmetric star, is written in the well-known form where x µ = (ct, r, θ, φ), and dΩ 2 = dθ 2 + sin 2 θdφ 2 is the line element on the unit 2-sphere. The metric functions ψ and λ, in principle, depend on x 0 and r.
With regard to the matter-energy distribution, we assume that the system is composed of an anisotropic perfect fluid, where the components of the pressure are not equal to each other, namely [14,18,82] with u µ being the four-velocity of the fluid, = c 2 ρ the energy density (where ρ denotes mass density), σ ≡ p t − p r the anisotropy factor, p r the radial pressure, p t the tangential pressure, and k µ is a unit spacelike four-vector. The four-vectors u µ and k µ must satisfy the following properties Within the context of GR, the spacetime geometry and the matter-energy distribution are related by the Einstein field equations where κ ≡ 8πG/c 4 , G µν is the Einstein tensor, R µν the Ricci tensor, and R denotes the scalar curvature. Here G is the gravitational constant and c is the speed of light in physical units.

A. Background and TOV equations
In the case of hydrostatic equilibrium none of the metric and thermodynamic quantities depends on the time coordinate x 0 , which entails that k µ = (0, e −λ , 0, 0) and hence the energy-momentum tensor contains only non-zero diagonal components T ν µ = diag(− , p r , p t , p t ). Consequently, from Eqs. (1)-(4) together with the conservation law of energy and momentum, the relativistic structure of an anisotropic star in the state of hydrostatic equilibrium is governed by the TOV equations where m(r) is the mass enclosed in the sphere of radius r. The metric function λ(r) is obtained by means of the relation As usual, we define the radius of the star when the radial pressure vanishes, i.e., the surface of the anisotropic star is reached when p r (r = R) = 0, and the total gravitational mass of the star is given by M ≡ m(R). Given a barotropic EoS of the form p r = p r (ρ) and a defined anisotropy relation for σ, Eqs. (5) and (6) can be integrated for a given central density and by guaranteeing regularity at the center of the star. Besides that, since the equilibrium system is a spherically symmetric star, the exterior spacetime of the anisotropic fluid must be described by the Schwarzschild metric so that the continuity of the metric on the surface imposes another boundary condition for the differential equation (7). Thus, the system of Eqs. (5)-(7) is solved under the requirement of the following boundary conditions B. Stability criterion through adiabatic radial oscillations In order to study the (in)stability of anisotropic neutron stars, it is necessary to calculate the frequencies of normal vibration modes. In fact, this involves an examination of radial perturbations in anisotropic fluid configurations. For such analysis we consider adiabatic vibrations, that is, we shall neglect the heat transfer between neighboring fluid elements.
Oscillation frequencies about the equilibrium state can be found by considering small deviations with respect to the state of hydrostatic equilibrium. In other words, the equilibrium configuration governed by the TOV equations (5)-(7) is radially perturbed in such a way that its spherical symmetry is maintained. Such a perturbation will cause motions in the radial directions so that a fluid element located at radial coordinate r in the unperturbed configuration is displaced to radial coordinate r + ξ(x 0 , r) in the perturbed configuration, where ξ is the Lagrangian displacement. This involves solving the perturbed Einstein equations δG µν = κδT µν for small radial oscillations from equilibrium. Therefore, in order to do a tractable analysis of the pulsations, all equations are linearized in the Eulerian perturbation functions δh, where the quantity h represents any metric or fluid variable and is decomposed as h(x 0 , r) = h 0 (r) + δh(x 0 , r). The quantities denoted by a subscript zero stand for the solutions in the equilibrium configuration.
It is important to note here that the relation between Eulerian and Lagrangian perturbations is given as follows where ∆h is the Lagrangian perturbation, that is, the change measured by an observer who moves with the fluid. If in the perturbed state we define v ≡ ∂r/∂x 0 = ∂ξ/∂x 0 , to first order in ξ the non-zero components of the energy-momentum tensor (2) take the form T 00 = e 2ψ , T 0r = T r0 = −( 0 + p r0 )e 2λ0 v, T rr = e 2λ p r , T θθ = r 2 p t , T φφ = r 2 p t sin 2 θ, (11) where the four-velocity and the unit four-vector are given by u µ = (e −ψ , ve −ψ0 , 0, 0) and k µ = (ve λ−2ψ , e −λ , 0, 0), respectively. Afterward, by retaining the terms only of first order in ξ and δh, the linearized field equations are In addition, the µ = r component for the fourdivergence of the energy-momentum tensor ∇ ν T ν µ = 0, provides the following relation (15) or alternatively, Through Einstein equations it was possible to express some perturbations in terms of the Lagrangian displacement ξ and the unperturbed variables. Nonetheless, we still need to have an expression for the perturbation δp r , and to obtain it, an additional condition is necessary (the conservation of the baryon number). If n is the number of baryons per unit volume, its conservation in GR is given by ∇ µ J µ = 0, where J µ ≡ nu µ is the baryon number current. At the same time, if we consider that the EoS has the general structure n = n( , p r ), we have where γ ≡ n pr dpr dn is a dimensionless quantity that determines the changes of radial pressure associated with variations in the particle number density. If the entropy is conserved, we obtain the adiabatic index given by We now assume that the Lagrangian displacement has a harmonic time dependence as ξ(x 0 , r) = ξ(r)e iαx 0 where cα ≡ ω is a characteristic frequency to be determined. The same procedure is applied for the metric functions and thermodynamic quantities. Thus, the substitution of Eq. (13) into (16), leads to By means of Eqs. (7), (14), (17) and the θθ-component of the field equations, the last expression can be written in terms of the unperturbed variables Because all terms are now the amplitudes of the perturbations and quantities of the static background, we can delete all reference to subscripts zero. Furthermore, let us define a new variable ζ ≡ ξ/r, so that Eqs. (17) and (20) can be rewritten as two first-order time-independent where the specific form of the perturbation δσ depends on the anisotropy ansatz that we are going to use afterwards. These equations govern the adiabatic radial oscillations in the interior of an anisotropic spherical star, and we highlight that such equations differ considerably from those obtained by the authors in Ref. [20]. We have a Sturm-Liouville type problem for determining the radial oscillation modes with eigenvalues ω 2 0 < ω 2 1 · · · < ω 2 n < · · · , where n stands for the number of nodes inside the anisotropic stellar fluid. It is evident that when σ = 0, the above system of equations assumes the Gondek's form for isotropic fluids [37,39,45].
To numerically solve Eqs. (21) and (22), it becomes necessary to specify some physically meaningful boundary conditions. Analogous to a vibrating string that is fixed at its endpoints, the radial pulsations in the inner region of a star occur between its center and the surface. Indeed, since Eq. (21) has a singularity at the origin, it is required that as r → 0 the coefficient of 1/r term must vanish, namely Meanwhile, at the stellar surface where p r (R) = 0, the appropriate boundary condition is that Lagrangian perturbation of the radial pressure vanishes, this is, Notice that the Lagrangian displacement must vanish at the center due to spherical symmetry, that is ξ(0) = 0. Nevertheless, our equations are in terms of ζ, so a particularly simple approach that is often adopted is to normalize the eigenfunctions so that ζ(0) = 1 at the center. In addition, we point out that in the treatment carried out by Misner et al. [83] for isotropic configurations, the radial oscillations are described by a second-order ordinary differential equation in the "renormalized displacement function" given by η ≡ r 2 ξe −ψ . In this regard, for anisotropic fluids we obtain the following differential equation governing the adiabatic radial pulsations where Eq. (25) leads to a self-adjoint eigenvalue problem for the oscillation frequencies. Indeed, when the anisotropy vanishes (i.e., p t = p r ), such equation is reduced to the form given in Refs. [38,83].

C. Equation of state and anisotropy ansatz
The conventional way to tackle the problem of anisotropic configurations is by specifying a barotropic EoS for radial pressure, i.e. p r = p r (ρ), and additionally an anisotropy function σ ≡ p t − p r must be assigned. Here, we are going to use a well known EoS and four different functions for σ.
Based on the SLy effective nucleon-nucleon interaction, Douchin and Haensel [84] calculated an unified equation of state (the so-called SLy EoS) that covers three main regions of neutron-star interior: outer crust, inner crust and liquid core (consisting of neutrons, protons, electrons and negative muons). Such an EoS is compatible with the constraints of the gravitational-wave event GW170817 (observed by the LIGO-Virgo detectors [27]), and the analytical parameterization of pressure as a function of density for non-rotating stars is as follows [85] A(y) = a 1 + a 2 y + a 3 y 3 1 + a 4 y K 0 (a 5 (y − a 6 )) +(a 7 + a 8 y)K 0 (a 9 (a 10 − y)) +(a 11 + a 12 y)K 0 (a 13 (a 14 − y)) +(a 15 + a 16 y)K 0 (a 17 (a 18 − y)), where A ≡ log(p r /dyn cm −2 ), y ≡ log(ρ/g cm −3 ), and K 0 (x) ≡ 1/(e x + 1). The fitting parameters of this expression a i can be found in Ref. [85]. Furthermore, below we specify the anisotropy functions provided in the literature to model anisotropic matter at high densities:

Quasi-local ansatz
Horvat et al. [17] suggested an anisotropy ansatz as being a bilinear function in the radial pressure and compactness, namely where β H is a dimensionless parameter that measures the degree of anisotropy within the star, and µ(r) ≡ 2Gm/c 2 r is known as compactness. The advantage of this ansatz is that the fluid becomes isotropic at the stellar center since µ ∼ r 2 when r → 0. On the other hand, (30) is important only for relativistic configurations which is in agreement with the assumption that the anisotropy may arise at high densities. A similar ansatz for the anisotropy measure was used to describe gravastar models in Ref. [86]. Moreover, in the calculations it is common to assume −2 ≤ β H ≤ 2 [4,17,26,82,87]. According to Eq. (8), the Eulerian perturbation for the anisotropy function (30), can be written as

Bowers-Liang ansatz
Another relation for anisotropic models was proposed by Bowers and Liang [14], given by where the anisotropy factor depends nonlinearly on the radial pressure and energy density. The anisotropy vanishes at the origin in order to yield regular solutions, and it is (in part) gravitationally induced since 1 − µ = e −2λ .
Here the literature offers a similar range for the free parameter β BL as in the first model mentioned above (see Refs. [26,87,88]). In this case, δσ assumes the form

Herrera-Barreto ansatz
An ansatz that was initially studied in Ref. [89] in order to find a family of non-isotropic configurations from any isotropic model, and later summarized by Herrera and Barreto [18], is the following where h is a constant throughout the sphere, and for h = 1 we recover the isotropic case. Here we are going to define β HB ≡ (h − 1)/2h so that σ = β HB rp r . The possibility (34), like the other ansatze, guarantees that the anisotropy must vanish at the center of symmetry of the fluid. In this case β HB cannot be positive in order for the tangential pressure to be always positive throughout the stellar interior. We must point out that physically relevant solutions correspond to p r , p t ≥ 0 for r ≤ R. The Eulerian perturbation for the ansatz (34) is given by

Covariant ansatz
Finally we will consider an additional ansatz that has recently been introduced by Raposo et al. [90] in order to study the dynamical properties of anisotropic self-gravitating fluids in a covariant framework, that is, σ ≡ −Cf ( )k µ ∇ µ p r = −Cf ( )e −λ p r , where C is a free constant that measures the deviation from isotropy. For simplicity, we are going to consider f ( ) = , and since we are using physical units, we have where now the free parameter is β R ≤ 0 and, unlike previous models, it has units of cubic meters. At the stellar origin the fluid becomes isotropic since the radial pressure is maximum at this point, and at the surface both pressures vanish. For this anisotropic model, δσ takes the following form

III. NON-ADIABATIC GRAVITATIONAL COLLAPSE
The goal in this section is to study the dynamical evolution of unstable anisotropic neutron stars whose ultimate fate is the formation of an event horizon, and to describe the formation of black holes it is necessary to use models that involve non-ideal fluids [91]. The gravitational collapse is a highly dissipative phenomenon in which massless particles (photons and neutrinos) carry thermal energy for exterior spacetime [92,93]. In this respect, we are going to deal with the problem following the standard procedure, namely, the spherical hypersurface Σ of the collapsing star divides the spacetime into two different regions where each one is described by a particular matter-energy distribution.

A. Interior spacetime and generalized TOV equations
We model the collapsing configuration by means of a locally anisotropic fluid, bounded by Σ, and that undergoes dissipation in the form of heat flow. Accordingly, in the diffusion approximation, the energy-momentum tensor is given by [92] where, as in the static background, represents the energy density, p r the radial pressure and p t the tangential pressure. The four-velocity u µ , the heat flux q µ and the unit four-vector along the radial direction k µ , must satisfy the following relations In order to introduce a time dependence on the metric functions such that under a certain limit we can recover the static case, we assume that the geometry of the interior spacetime is described by the following spherically symmetric, shear-free line element [45] where x µ − = (ct, r, θ, φ) are the coordinates in the interior manifold and df /dx 0 − < 0 in the case of collapsing configurations. From relations (39) and by using comoving coordinates, we get being q = q(x 0 − , r) the rate of energy flow per unit area along the radial coordinate. In the static limit f (x 0 − ) → 1 we recover Eq. (1) which describes the initially static anisotropic star.
The explicit form of the Einstein field equations for the energy-momentum tensor (38) and metric (40), is given by where the dot and the prime denote differentiation with respect to x 0 − and r, respectively. The mass entrapped within the radius r and at time t, is given by the following expression [92,94] As a consequence, the generalized TOV equations within a non-adiabatic context are derived from the fourdivergence of the energy-momentum tensor (38) along with Eq. (43), that is, where it should be noted that now the energy density, radial pressure, heat flux, anisotropy ansatz, and mass function depend on both t and r, whereas ψ = ψ(r), λ = λ(r) and f = f (t). It becomes apparent that at the static limit such equations are reduced to those already known in (6) and (7), respectively. Furthermore, we remark that Eq. (46) is the dissipative analogue of the relation (8).

B. Junction conditions on the stellar surface
Since the star is radiating energy, the exterior spacetime is described by the Vaidya metric [56,95], given as follows where x µ + = (cυ, χ, θ, φ), m(υ) is the mass function that depends on the retarded time υ, and e = ±1 describes the incoming (outgoing) flux of radiation around the source of gravitational field. In our collapse model the radiation is expelled into outer region so that dm/dυ ≤ 0.
Since the collapsing star is described by two spacetime regions with distinct geometric properties, it becomes necessary to invoke the junction conditions on Σ established by Israel [96,97]. Such conditions require continuity of the line element and extrinsic curvature through the hypersurface, namely where the intrinsic metric to Σ is given by The extrinsic curvature tensor is defined by where x µ ± are the coordinates of the exterior and interior spacetime, ς i = (cτ, θ, φ) are the coordinates that define the comoving timelike hypersurface, and n µ ± are the unit normal vectors to Σ which have already been calculated by Santos [59]. The non-vanishing extrinsic curvature components K ± ij are given explicitly in Appendix A. As a result, the junction conditions (50) imply that with z Σ being the boundary redshift of the radial radiation emitted by the non-adiabatic sphere. Eq. (53) is the equality of the proper radii as measured from the perimeter of Σ. The expression (54) is a measure of the total mass of the star as it collapses, and Eq. (56) indicates that the radial pressure at the surface of the star is different from zero unless the heat flow vanishes. It is evident that in the static limit m Σ corresponds to the total mass of the initial Schwarzschild configuration M , and the redshift is reduced to z Σ = e λ(R) − 1, namely, the gravitational redshift of light emitted at the surface of the initially static neutron star.

C. Evolution quantities
The fact that the radial pressure does not vanish at the surface leads to an additional differential equation that allows us to fix the time evolution of our model. By taking into account Eqs. (43) and (45) into (56), we obtain which can be integrated to generate the following equa- where the integration constant has been determined by applying the static limit. Then the solution of Eq. (58) is given by For an observer at rest at infinity, the redshift (55) diverges at the time of formation of an event horizon. This means that a black hole has been formed as outcome of the gravitational collapse of an unstable anisotropic star when For systems describing gravitational collapse we must have 0 < f ≤ 1, i.e. the time function f decrease monotonically from f = 1 to f = f bh . In other words, the time goes from t = −∞ (when the model is static) to t = t bh (when the star becomes a black hole), but a time displacement can be done without loss of generality. Consequently, one obtains from Eq. (54) the mass of the formed black hole, which reads: Notice that Eq. (59) provides t as a function of f , nonetheless, is more useful to obtain f (t) in order to analyze the dynamical quantities as a function of time as a stellar configuration collapses. Thus, it is convenient to numerically solve Eq. (57) as a final value problem by specifying a value of f (t) and df (t)/dt at time t = t bh , where the two final conditions are established through Eqs. (58)- (60). The relevant physical quantities during the collapse such as energy density, radial pressure, tangential pressure and heat flow, are given by where a ≡ GM/c 2 R 2 , and the anisotropy factor takes the following form Finally, a kinematic quantity that provides information about the rate of expansion of the fluid sphere is given by the four-divergence of the four-velocity [13,92] that is, the expansion scalar and whose action is to change the volume of the spherical star but it preserves the principal axes. In fact, when f → 1, the heat flow and the expansion scalar vanish and, therefore, the fluid becomes perfect. The other physical quantities are reduced to those already known in the static background.

D. Energy conditions
The dissipative anisotropic fluid must satisfy the energy conditions throughout the gravitational collapse in order for it to be physically acceptable. This means that the energy-momentum tensor (38) has to be diagonalized through equation |T − µν − Υg − µν | = 0, so that the eigenvalues Υ take the explicit form where we defined ∆ ≡ ( + p r ) 2 − 4q 2 andq ≡ q c e λ √ f . Thus, the following energy conditions must hold [98] Weak energy conditions (WEC) Such inequalities entail that c) Υ 0 + Υ i ≤ 0, for i = 1, 2, 3.
The first two inequalities have already been included in weak energy conditions. As regards the third inequality, we obtain Strong energy conditions (SEC) Specifically, the first inequality implies that whereas the second inequality has been considered in the other energy conditions. Therefore, to know whether our collapse model is physically acceptable, it is only necessary to verify that the energy conditions (71), (74), (75) and (76) are respected.

A. Equilibrium configurations and radial oscillations
Given different values of central density and an anisotropy ansatz with EoS (29), the background equations (5)-(7) can be numerically solved under the conditions (9) to produce a family of anisotropic neutron stars in hydrostatic equilibrium. Figure 1 illustrates the landscape of these configurations for each anisotropy ansatz. In the case of the ansatze proposed by Horvat et al. (30) and Raposo et al. (36), for low enough central densities, the configurations have a mass similar to the isotropic case. Nevertheless, above a certain central density value, the masses deviate considerably from those provided by isotropic solutions. It is evident that depending on the degree of anisotropy within the configurations, the curves in the M versus R plane can move significantly away from the isotropic case.
One of the most interesting approaches to determining the neutron-star matter EoS is through measurements of the masses and radii of these stars. The observations made in the last few years are allowing us to improve our understanding of the properties of cold dense matter and, therefore, to constrain the EoS. In that regard, we expect that the effects generated by anisotropic pressure are also within the constraints obtained by recent observations. Indeed, we consider one of both mass and radius measurements for the millisecond pulsar PSR J0030+0451 from Neutron Star Interior Composition Explorer (NICER) data, which was obtained by using a Bayesian inference approach to analyze its energy-dependent thermal X-ray waveform [99]. According to Fig. 1, it is possible to construct anisotropic neutron stars that lie within the region provided by the NICER data. Furthermore, from the observation of the GW event GW170817 there is a recent restriction on the maximum mass of neutron stars [100]. To ensure that our results are consistent with such a restriction, bounds on anisotropy parameters must be established. This hints that β H 0.40, β BL 0.23, β HB −0.042, and β R −4.6 × 10 11 m 3 .
In table I, we list the mass and radius of the maximummass configurations for each anisotropic model with different values of the free parameter. According to the M (ρ c ) method, the first maximum on the M (ρ c ) curve corresponds to a critical central density ρ c = ρ crit which delimits a family of stars that is stable against gravitational collapse. This means that the unstable branch in the sequence of stars is located after the critical density where dM/dρ c < 0. Due to its simplicity, this condition has been widely used in the literature. However, such a condition is just necessary but not sufficient to determine the limits of stability.
Once the equilibrium quantities are known by solving the TOV equations, our second task is to verify if the M (ρ c ) method is compatible with the calculation of frequencies. To that end, we have to solve the system of coupled first-order equations (21) and (22) with boundary conditions (23) and (24). The numerical solution of these equations is carried out using the shooting method, that is, we integrate the equations for a set of trial values of ω 2 satisfying the condition (23). Besides that, we consider that normalized eigenfunctions correspond to ζ(0) = 1 at the center, and we integrate to the stellar surface. The values of the squared frequency for which the boundary condition (24) is satisfied are the correct frequencies of the radial oscillations. In particular, for a central mass density ρ c = 2.0 × 10 18 kg/m 3 with anisotropy parameters β H = 0.50, β BL = 0.25, β HB = −0.05 and β R = −5.0×10 11 m 3 , we display in Fig. 2 the Lagrangian perturbation of the radial pressure for a set of test values ω 2 , where each minimum indicates the appropriate frequency. In other words, for a given stellar configuration there are different eigenvalues ω 2 n with their respective eigenfunctions ζ n (r) and ∆p r,n (r), where n represents the number of nodes inside the anisotropic star. In Fig. 2, the first (leftmost) minimum represents the fundamental oscillation mode and it has no nodes between the center and the surface, whereas the first overtone (n = 1) has a node, the second overtone (n = 2) has two, and so forth. These radial pulsations in the anisotropy ansatz can be visualized in more detail in Fig. 3 and whose fundamental mode oscillation frequencies are shown in table II.
Using the central density as a parameter, the family of anisotropic stars that are actually stable is shown in the left plot of Fig. 4. The squared frequency of the fundamental oscillation mode against the central density is shown in the right plot of the same figure. Unlike the case of strange stars with MIT bag model EoS (where the frequency of the fundamental mode always decreases with increasing central density), in neutron stars ω 2 0 increases to a maximum value and then decreases with ρ c regardless of the anisotropic model. In addition, for larger values of β H and β BL , the onset of instability is indicated at a smaller and smaller central density value. On the other hand, for larger values of β HB and β R , the onset of instability is found at a greater and greater central density value as we approach the isotropic case.
Taking into account the data recorded in table I, we see that only for the ansatze proposed by Horvat et al. (30) and Raposo et al. (36), the onset of instability indicated by the M (ρ c ) method is located exactly at the configuration that has vanishing frequency of the fundamental mode. In other words, for these two anisotropy profiles, the maximum-mass point M max and ω 2 0 = 0 are reached at the same central density value. Nevertheless, for the anisotropic models suggested by Bowers-Liang (32) and , the squared frequency of the fundamental mode does not pass through zero at the critical central density corresponding to the maximum-mass configuration. Therefore, it is evident that anisotropy affects the stellar stability and the critical central density does not always correspond to the onset of instability. Anisotropic neutron stars in hydrostatically stable equilibrium oscillate with a purely real (fundamental) frequency when are subjected to a radial perturbation, whereas the unstable stars (with imaginary frequency of the lowest oscillation mode) undergo a gravitational collapse from rest to form a black hole. We assume that the unstable configurations are initially in a state of hydrostatic equilibrium and then gradually begin to collapse until the formation of an event horizon. In the case of unstable anisotropic neutron stars with SLy EoS and anisotropic model proposed by Bowers and Liang (32), for an initial central mass density ρ c = 2.6 × 10 18 kg/m 3 and anisotropy parameter β BL = 0.2, we solve Eq. (57) with final conditions f = 0.150 and df /dt = −8.613 × 10 3 s −1 at time t = t bh = −1.723 × 10 −5 s. For this particular configuration, the tangential pressure dominates the radial pressure and ω 2 0 = −13.001 × 10 6 s −2 . Then we perform a time displacement so that this unstable star evolves from the initial instant t = 0 (when the interior structure is governed by the background equations and the external solution is Schwarzschild-type) until the moment of horizon formation t bh = 1.981 ms, that is, when the star has collapsed and the mass of the resulting black hole is m bh = 1.335 M .
The energy density (62) and radial pressure (63) as functions of the radial coordinate at different times are displayed in the upper and lower panels of Fig. 5, respectively. Both thermodynamic quantities present their maximum values at the stellar center and change significantly in the last moments of the collapse, while near the surface the changes are relatively small. Unlike the static case, the radial pressure at the surface is no longer zero during the dynamical evolution of the gravitational collapse because there exist a radial heat flux according to Eq. (56). As a result, we can investigate how the EoS behaves as the star collapses. The upper panel of Fig.  6 reveals that the EoS for the radial pressure undergoes sudden changes during the collapse of an unstable neutron star, where the maximum and minimum values in each curve correspond to the center and the surface of the star, respectively.
In the lower panel of Fig. 6 we plot the radial heat flux (65), which undergoes great alterations in the intermediate regions of the collapsing configuration and its value is not zero at the surface. Indeed, when the heat flux vanishes, the radial pressure also vanishes at the surface and the exterior solution is the Schwarzschild vacuum solution. According to the upper panel of The isotropic case is shown in all plots as a benchmark by a solid black line, and the cigar-shaped yellow region is one of both mass and radius measurements for PSR J0030+0451 from NICER data [99]. The horizontal narrow band in cyan color stands for the recent restriction of maximum mass of neutron stars as a result of observation of the GW event GW170817 [100].
In the lower panel of Fig. 7 we illustrate the radial behaviour of the expansion scalar (67) during the collapse process. For any instant of time, it can be seen that Θ < 0 and ∂Θ/∂r > 0, which means that the stellar system is collapsing. In addition, it is clear that Θ → 0 as t → 0.
At long last, the collapsing anisotropic neutron star with initial central mass density ρ c = 2.6 × 10 18 kg/m 3 , SLy EoS for radial pressure and anisotropy parameter β BL = 0.2, is physically reasonable because it obeys the energy conditions in the full extent of the star and throughout the collapse process. It is worth emphasizing that we have tested this procedure for the other anisotropy ansatze (30), (34) and (36), obtaining a similar behaviour during the dynamical evolution of the col-lapse.

V. CONCLUSIONS
In this paper, we have constructed families of anisotropic neutron stars for an EoS compatible with the recent observation of the event GW170817 and for four different anisotropy ansatze, aiming to settle bounds on the maximum masses reachable in this kind of anisotropic models. We have carried out an analysis of adiabatic radial pulsations for such stars in order to study their radial stability against gravitational collapse. We have also developed a dynamical model that describes the non-adiabatic gravitational collapse of the unstable anisotropic configurations, and TOV equations have been generalized within this context. To summarize, the main conclusions of this work are: For the SLy EoS in the radial pressure, we have constrained the anisotropy parameters provided in the literature in order to satisfy the recent restriction of the maximum mass of neutron stars based on gravitational wave observations. Our results thus suggest that β H 0.40, β BL 0.23, β HB −0.042, and β R −4.6 × 10 11 m 3 . However, we must point out that these bounds can be slightly altered for other EoSs.
Anisotropy affects the stellar stability and the critical central density (where dM/dρ c = 0) does not always correspond to the onset of instability. In particular, the maximum point on the M (ρ c ) curve does not indicate the onset of instability for the anisotropic models proposed by Bowers-Liang [14] and Herrera-Barreto [18]. Nevertheless, we remark that this criterion is compatible with the calculation of the oscillation frequencies in the case of the ansatze suggested by Horvat et al. [17] and Raposo et al. [90].
Given an initial value of central mass density and a specific anisotropy profile for an unstable anisotropic neutron star, we have investigated the evolution of the equation of state for radial pressure and anisotropy ansatz as the star undergoes a non-adiabatic gravitational collapse. The sudden changes in all relevant physical quantities occur near the formation of the event horizon as a consequence of a radial heat flow.