General Relativistic Stability and Gravitational Wave Content of Rotating Triaxial Neutron Stars

Triaxial neutron stars can be sources of continuous gravitational radiation detectable by ground-based interferometers. The amplitude of the emitted gravitational wave can be greatly affected by the state of the hydrodynamical fluid flow inside the neutron star. In this work we examine the most triaxial models along two sequences of constant rest mass, confirming their dynamical stability. We also study the response of a triaxial figure of quasiequilibrium under a variety of perturbations that lead to different fluid flows. Starting from the general relativistic compressible analog of the Newtonian Jacobi ellipsoid, we perform simulations of Dedekind-type flows. We find that in some cases the triaxial neutron star resembles a Riemann-S-type ellipsoid with minor rotation and gravitational wave emission as it evolves towards axisymmetry. The present results highlight the importance of understanding the fluid flow in the interior of a neutron star in terms of its gravitational wave content.


Introduction
One of the most profound predictions of general relativity is that a system which possesses time varying multipole moments higher than a quadrupole, generates gravitational waves.The most common systems that satisfy such criterion are the ones that not symmetric about their rotation axis, with prime examples being those of binary compact objects.Therefore it is not a surprise that the first direct detection of a gravitational wave came from a binary black hole [1,2].In the first three observational periods (O1-O3), the LIGO/Virgo [3,4] collaboration has discovered gravitational waves from almost 100 binaries [5,6,7] including two binary neutron stars [8,9,10,11] and two black hole-neutron star mergers [12].Another exciting possibility is to detect gravitational waves from a single neutron star that exhibits some kind of asymmetry [13,14,15,16].Although such gravitational waves are much weaker than the ones emerging from a binary system 1 arXiv:2312.16728v1[gr-qc] 27 Dec 2023 (and this is one of the reasons that they have not been detected yet), they have the potential of providing important information regarding the nature of a neutron star such as regarding various fluid instabilities or its elastic, thermal and magnetic characteristics.
A hydrodynamical instability is one well-known mechanism that can produce nonaxisymmetric neutron stars which emit gravitational waves [17].One important parameter that characterizes unstable rotating neutron stars is β := T /|W |, where T is the rotational kinetic energy and |W | the gravitational binding energy [17,18].As the rotation of the star increases, there are two critical points (nonaxisymmetric instabilities) that are associated with two different physical mechanisms.In the presence of some dissipative mechanism such as viscosity or gravitational radiation, at β = β s the star becomes secularly unstable to a bar-mode deformation.The timescale of this instability is set by the dissipation and is much longer than the dynamical (free-fall) timescale.At even higher rotation rates, when β = β d > β s the star becomes dynamically unstable to a bar-mode deformation.This instability emerges regarless of any possible dissipation and its growth is set by the dynamical timescale.For incompressible stars in Newtonian gravity β Newt s = 0.1375 and β Newt d = 0.2738 [19].Although the values of β at these critical points can change in general relativity, with compressible equations of state and differential rotation, the overall idea (existence of distinct secular and dynamical instability points) remains 1 .
The broadbrush picture above can be further refined by the fact that there are two categories of secular instabilities: i) The viscosity-driven instability which as the name suggests manifests itself in the presence of viscous dissipation [23], and ii) the Chandrasekhar-Friedman-Schutz (CFS) instability which is driven by gravitational radiation reaction [24,25,26].For Newtonian incompressible fluids, an axisymmetric rotating body is described by a Maclaurin spheroid [19], an oblate spheroid having R x = R y ̸ = R z .At the point of secular instability, when β = 0.1375, two families of triaxial (R x ̸ = R y ̸ = R z ) solutions emarge: a) The Jacobi ellipsoids, which are uniform rotating ellipsoidal figures of equilibrium in the inertial frame, and thus emit gravitational waves, and b) the Dedekind ellipsoids, which are ellipsoidal figures of equilibrium stationary in the inertial frame, and therefore do not emit gravitational waves 2 .The Dedekind ellipsoids have constant vorticity and nonzero internal fluid circulation.Equilibrium solutions (a) and (b) are related to the processes (i) and (ii) respectively as follows [27,28,29].Viscosity dissipates mechanical energy but conserves angular momentum and a Jacobi ellipsoid has less mechanical energy, T + W , than a Maclaurin spheroid of the same rest mass and angular momentum.Thus the viscous-driven evolution (i) of an unstable Maclaurin spheroid would proceed towards a Jacobi ellipsoid (a).On the other hand, gravitational radiation preserves circulation along any closed path on a plane parallel to the equator, but not angular momentum.A Dedekind ellipsoid has less mechanical energy, than a Maclaurin spheroid of the same rest mass and circulation.Thus, in the absence of viscosity, the CFS-driven (ii) evolution of an unstable Maclaurin spheroid would proceed towards a Dedekind ellipsoid (b).The presence of both viscosity and gravitational radiation tends to stabilize the star against these competing mechanisms [30].In the limit where the gravitational wave timescale equals the viscous timescale, the Maclaurin spheroid is secularly stable all the way to the dynamical instability point.
One important difference between the viscosity-driven instability and the CFS instability is that the latter is generic while the former is absent in sufficiently slowly rotating stars.In addition the viscosity-driven instability emerges only for sufficiently stiff equations of state in which the bifurcation point exists before the mass-shedding (Keplerian) limit.In Newtonian gravity with a polytropic equation of state, p = kρ Γ , the triaxial sequence exists only if Γ ⪆ 2.24 [31].In general relativity the critical adiabatic index does not change significantly but slightly increases ∼ 2.4 [32,33,34,35].At the same time the critical parameter β s also increases relative to the Newtonian value (0.1375) by a factor that depends on the compactness of the neutron star [36].On the other hand, the CFS instability becomes stronger in general relativity and sets in at β < 0.1375 [37,38] so that the two instabilities no longer occur at the same value of β s .
Sequences of triaxial solutions in general relativity have been investigated in [39,40,41] using a select set of stiff equations of state.It was found that the triaxial sequence becomes shorter (a smaller deformation is allowed) as the compactness increases, while supramassive [42] triaxial equilibria are possible, depending on the equation of state.In [43] the first full general relativistic simulations of triaxial uniformly rotating neutron stars have been performed, and the dynamical stability of certain normal and supramassive solutions was established.It was found that all triaxial models evolve toward axisymmetry, maintaining their uniform rotation, while losing their triaxiality through gravitational wave emission.Similar results were reported in [44] where triaxial quark stars (having finite surface density) were evolved in general relativity.
In this work we investigate the fate and stability of triaxial models against a variety of perturbations.First we establish the dynamical stability of the most triaxial figure of quasiequilibrium along two constant rest-mass sequences, one that corresponds to compactness 0.1, and another one that corresponds to compactness 0.19.Second, by replacing the Jacobi-like velocity flow with a Dedekind-like one we explore the fate of the resulting ellipsoidal neutron star.We find that in some cases this procedure leads to a Riemann-S-type ellipsoidal figure of quasiequilibrium that barely rotates while largely preserving its nonaxisymmetric shape.This object emits gravitational waves whose amplitude is approximately 20% the one coming from the original triaxial neutron star as it evolves towards axisymmetry, and highlights the importance of the fluid flow in accurate gravitational wave analysis.

Numerical methods and model parameters
For the construction of the initial models we use the cocal code as described in [39,40,41] while for the evolution we use the Einstein Toolkit [45,46,47,48,49].Below we summarize the most import features of our initial data models and their evolutions.

Initial data
We construct uniformly rotating triaxial neutron stars having angular velocity Ω and velocity with respect to the inertial frame v i = Ωϕ i = Ω(−y, x, 0).The fluid's 4-velocity can be written as where u t is a scalar.The spacetime of the rotating star possesses a helical Killing vector, k α , where with the fluid variables being Lie-dragged along k α , Here ρ, h, s are the rest-mass density, enthalpy, and the entropy per unit rest-mass.We have ρh = ϵ + p, where ϵ is the total energy density and p is the pressure.
In order to ensure the existence of triaxial uniformly rotating models we use a stiff polytropic equation of state with Γ = 4.For the polytropic constant we chose k = 1 in G = c = M ⊙ = 1.Similar to [41,43] the value of Γ used is simply to prove a point of principle, rather than to address physical EOS parameters.
The models have been computed with the cocal code, a second-order finite-difference code whose methods are explained, for example, in [40,50].For computational convenience, we employ the Isenberg-Wilson-Mathews (IWM) formulation [51,52,53,54,50].Therefore the 3-metric is γ ij = ψ 4 f ij , where ψ is the conformal factor and f ij the flat metric.The unknown gravitational variables in the 3+1 formulation are the lapse α, the shift β i , and the conformal factor ψ. In the cocal code the full system of equations (wavless formulation) is also used but the differences from the conformally flat IWM scheme are small [40].A number of diagnostics are used to describe the initial solutions and explicit formulae are given in the appendix of [40].The most important diagnostics are: 1) The angular momentum of the triaxial star J (where J is the Arnowitt-Deser-Misner (ADM) angular momentum), which is computed via a surface integral at infinity or a volume integral over the spacelike hypersurface.2) The kinetic energy, which is defined as T := 1  2 JΩ (although we are not in axisymmetry we still use this formula because it is gauge-invariant), and 3) the gravitational potential energy, which is defined as W := M − M p − T .Here M is the (ADM) mass and M p is the proper mass (rest-mass plus internal energy) of the star (see e.g.[55]).These expressions are used then to compute the rotation parameter β = T /|W |.We also define the moment of inertial as I := J/Ω.
To measure of accuracy of the initial data we use two diagnostics: The first one is the difference between the Komar and ADM mass [40], For stationary and asymptotically flat spacetimes M K = M [56].The second diagnostic is the relativistic virial equation [57].In our calculations both diagnostics are ∼ O(10 −4 ).   2 , Ω, M , M 0 , (M/R) s , T /|W |, I, are the central rest-mass density, the coordinate radii, the proper eccentricity with respect to the z-axis, the angular velocity, the ADM mass, the rest mass, the corresponding (with the same rest mass) spherical compactness, the ratio of kinetic over gravitational potential energy, and the moment of inertia, respectively.To convert to cgs units, use the fact that 1 = 1.477 km = 4.927 µs = 1.989 × 10 33 g.We start our calculations by computing two sequences of axisymmetric rotating neutron stars having constant rest mass that correspond to spherical compactions (M/R) s = 0.1 and 0.19.These sequences which are shown with orange ((M/R) s = 0.1) and green ((M/R) s = 0.19) color in Fig. 1, are the analogues of the Maclaurin spheroids in Neutonian gravity [19].In the left panel of Fig. 1, the mass is plotted against the central density, while in the right panel we plot the the rotation parameter β versus the eccentricity e := 1 − ( Rz / Rx ) 2 with respect to the proper radii.Note that the mass and the density can by rescaled to any number using the polytropic constant k3 , hence the axes in the left panel are normalized accordingly.In the left panel of Fig. 1 we also show the sequence of spherically symmetric Tolman-Oppenheimer-Volkov (TOV) solutions (black curve), and the sequence of maximally uniformly rotating (at the mass-shedding limit, also called the Kepler limit) solutions (red curve).With black (red) circle we denote the solution of maximum nonrotating (uniformly rotating) mass.The shaded area corresponds to densities where the speed of sound c s is larger than the speed of light.All the solutions used in this work are causal.
For sufficiently high rotation rate (β) a second branch of solutions appears.These are triaxial solutions (R x ̸ = R y ̸ = R z ) that correspond to the Newtonian Jacobi ellipsoids [19].The two sequences that correspond to (M/R) s = 0.1 and 0.19 are shown with blue ((M/R) s = 0.1) and magenta ((M/R) s = 0.19) colors.In the left panel of Fig. 1 all triaxial solutions that correspond to each sequence have masses and densities very close to each other so they appear as a single triangle point (magenta or blue).In the right panel though, the triaxial sequences are clearly seen.For a fixed eccentricity a triaxial model has less T /|W | than the corresponding axisymmetric model.In particular for a fixed eccentricity the triaxial solution has less gravitational (ADM) mass M (thus more negative gravitational potential energy W ), angular momentum J , angular velocity Ω, and moment of inertia I than the corresponding axisymmetric model.Therefore it has less kinetic energy too.On the other hand it has larger proper mass, hence less total energy T + W = M − M p 4 and thus it is the preferred figure of equilibrium.
From the left panel of Fig. 1 we notice that the bifurcation point happens at larger β or e as the compactness increases.The triaxial sequence also shrinks the larger the compactness, which intuitively means that it is harder to construct a triaxial neutron star of large compactness.For  incompressible fluids [36] in general relativity it was found that where β Newt s = 0.1375 at eccentricity e Newt s = 0.8127.Equation ( 5) predicts that β s = 0.15 at (M/R) s = 0.1 while for (M/R) s = 0.19, it is β s = 0.166, which are in broad agreement with the right panel of Fig. 1.Notice also that the IWM formulation sligthly overestimates β s as well as e s at the bifurcation point with respect to a full solution to the Einstein equations [40].
The models used in this study are shown in the right panel of Fig. 1 as blue and magenta dots.They constitute the most triaxial solutions along the corresponding sequences of constant rest mass.In Table 1, these two solutions are dubbed as C010s17 (magenta corresponds to compactness 0.1) and C019s08 (blue corresponds to compactness 0.19).
We have employed the single star module of the cocal code to compute the quasi-equilibrium solutions of this work.This module uses the KEH method [58,59] on a single spherical patch (r, θ, ϕ) with r ∈ [r a , r b ], θ ∈ [0, π], and ϕ ∈ [0, 2π], where r a = 0, r b = O(10 6 M ) (no compactification used), to achieve convergence through a Green's function iteration.The grid structure in the angular dimensions is equidistant but not in the radial direction.The definitions of the grid parameters can be seen in Table 2, along with the specific values used here.

Evolutions
For the evolution we use the baikal [60] code, which solves the Einstein field equations in the BSSN formalism and the Illinois GRMHD [46,47] to evolve fluid quantities.The code is built on the cactus infrastructure and uses carpet [48] for mesh refinement, which allows us to focus numerical resolution on the strong-gravity regions, while also placing outer boundaries at large distances well into the wave zone for accurate GW extraction and stable boundary conditions.The evolved geometric variables are the conformal metric γij , the conformal factor ϕ, (γ ij = e 4ϕ γij ), the conformally-rescaled, tracefree part of the extrinsic curvature, Ãij , the trace of the extrinsic curvature, K, and three auxiliary variables Γi = −∂ j γij , a total of 17 functions.For the kinematical variables we adopt the puncture gauge conditions [61,62], which are part of the family of gauge conditions using an advective "1 + log" slicing for the lapse, and a 2 nd order "Gamma-driver" for the shift [63].
The equations of hydrodynamics are solved in conservation-law form adopting high-resolution shock-capturing methods [64] 3: Grid parameters used for the evolution of each model.Parameter N corresponds to the number of points used to cover the largest radius of the star.Parameter dx is the step interval in the coarser level.density, ρ, the pressure p and the coordinate three velocity v i = u i /u 0 .The enthalpy is written as h = 1 + e + p/ρ, and therefore the stress energy tensor is T αβ = ρhu α u β + pg αβ .Here, e is the specific internal energy 5 .
The grid structure used in these evolutions is summarized in Table 3. Typically we use five refinement levels with the innermost level half-side length being approximately ∼ 1.5 times larger than the radius of the star in the initial data (R x ).We use 240 × 240 × 240 cells for the innermost refinement level, which means that we have approximately 160 points across the neutron star largest diameter.(For the initial data construction we used 256 points across the largest neutron star diameter.)For the innermost refinement level this implies a ∆x ∼ 5.53 × 10 −3 (C010s17) and ∆x ∼ 5.60 × 10 −3 (C019s08).This number of points was necessary in order to have accurate evolutions of such stiff equations of state (Γ = 4) which present a challenge for any evolution code.

Results
We perform full general relativistic simulations of the two most triaxial models C010s17 and C019s08 under a variety of perturbations in order to probe their stability and more importantly their fate especially with respect to their nonaxisymmetric shape.As a first test for dynamical stability we evolve these triaxial models by applying a 5% pressure depletion in their interior.Note that the dynamical timescale is a fraction of the period P of rotation of any system In Fig. 2 the maximum (central) density oscillations versus time are shown for five rotation periods.Overall, both models in Table 1 show the same oscillatory behavior when we pressure-deplete them, thus they are stable against quasiradial perturbations on dynamical timescales.Since these are the most triaxial models along a sequence of constant rest mass, the presented triaxial sequence (quasiequilibria along magenta or blue lines in right panel of Fig. 1) would also be stable.
Having established the dynamical stability against quasiradial perturbations, we now focus on the velocity flow of the triaxial figures of quasiequilibrium and investigate how it affects their global hydrodynamical stability.Let us recall that in Newtonian gravity the velocity of a Riemann-S ellipsoid in the inertial frame is where Ω is the angular velocity of the ellipsoid and Λ the angular frequency of the internal fluid circulation [19,27].When there is no internal fluid circulation (Λ = 0) the fluid velocity describes the velocity field of a Jacobi ellipsoid with vorticity ζ = −2Ω, while when there is no angular velocity (Ω = 0) the fluid velocity describes the velocity field of a Dedekind ellipsoid with vorticity ζ = −( Rx Ry + Ry Rx )Ω.
Figure 3: Left panel: Solid (dotted) colored lines are density isocontours at t/P = 30.6(t/P = 0) for the model C010s17-A075.The velocity field (red arrows for t/P = 30.6 and black ones for t/P = 0) is also shown.Right panel: The m = 2 mode amplitude for triaxial model C010s17 as well as for all its velocity perturbed models.Dashed lines are linear fits.
Since models C010s17 and C019s08 are the analogues of Jacobi ellipsoids in general relativity with a compressible equation of state, we replace their velocity flow field by the corresponding one of Dedekind ellipsoids.By substituting in Eq. ( 7) Ω = 0 and Λ = Ω, we find that the star significantly destabilizes, hence we follow the procedure below.First we identify the contours of constant rest-mass density, and construct their tangential directions.We then assign at each point a velocity whose direction is the one computed from the tangent to the isocontours and its magnitude is given by |Λ|((yR x /R y ) 2 + (xR y /R x ) 2 ) 1/2 , where Λ is a free parameter.In this way we ensure that the velocity field is consistent with the density gradients and only its magnitude can cause significant deformations.
Setting A = Λ/Ω we find that for the model C010s17 and A = 0.7, 0.75, 0.8, the rotational motion of the triaxial figure is greatly reduced but significant nonaxisymmetric oscillations are present.We refer to these evolutions as C010s17-A070, C010s17-A075, and C010s17-A080 respectively.In the left panel Fig. 3 dotted colored lines show the density isocontours at t = 0 with the velocity field constructed in the way explained above for the C010s17-A075 model.Also shown are the isocontours and velocity field at the end of our simulations at t/P = 30.6.In the right panel of Fig. 3 we plot the m = 2 nonaxisymmetric mode amplitude for all four models (C010s17 and it perturbed models) normalized by the rest mass of the system.For the nonperturbed case, C010s17, the amplitude of C 2 is monotonically decreasing until the end of our simulations.This is due to the emission of gravitational waves and its loss of energy and angular momentum, which make the neutron star more axisymmetric.Notice that a nonrotating triaxial ellipsoid that preserved its shape (as an equilibrium Dedekind ellipsoid) would have a constant C 2 for all times in a simulation.The perturbed cases C010s17-A070, C010s17-A075, and C010s17-A080 show an initial increased bar mode which decays in different ways.In all four cases we also plot with dashed lines linear fits.As we can see all three perturbed models evolve towards axisymmetry with different rates.The model that clearly shows the least amount of triaxiality change (for t/P ≳ 10) is C010s17-A070.
In order to appreciate the overall motion of these ellipsoidal figures, we plot in Fig. 4, left column, the density profile of the nonperturbed case C010s17 at a select number of times t/P = 0, 0.25, 0.5, 0.75, 1.0, 20.0.In the middle column we plot for the same times model C010s17-A070, while at the right column we plot model C010s17-A075.As can seen from the first five instances (t/P ≲ 1) where the nonperturbed model C010s17 makes one complete revolution, models C010s17-A070 and C010s17-A075 barely rotate while they exhibit nonaxisymmetric deformations.By the end of our simulations at t ∼ 30P all models remain nonaxisymmetric (although less than at t = 0) and continue to oscillate mildly.
In the left panel of Fig. 5 we plot the gravitational wave strain h × for the C010s17 models following [65].It is obvious that the Dedekind-like velocity flow decreased the gravitational wave signature of the triaxial figures by less than half, even at early times.Consistent with Fig. 4 and the left panel of Fig. 3 we see that the model with the least amount of gravitational wave content is C010s17-A070 (orange line) whose strain exhibits a decrease at ∼ 20% of the original C010s17 case (blue line).Note that since the perturbed models are barely rotating these high frequency gravitational waves are due to the hydrodynamical flow perturbations in the neutron stars, which showcases the importance of accurate hydrodynamical modeling for the understanding of a physical system through gravitational waves.In the right panel of Fig. 5 we plot the power spectrum of the C010s17 models scaled for a 1.4 M ⊙ triaxial neutron star mass at 50 Mpc, along with the noise curves for Advanced LIGO's design sensitivity [66] and the ET-B configuration of Einstein Telescope [67].The peak frequency of the unperturbed C010s17 model (blue curve) at ∼ 840 Hz is consistent with its orbital angular frequency (Ω/π) and in principle detectable by Advanced LIGO.Consistent with the left panel, the power spectrum of the perturbed models is weaker but still detectable from next generation observatories such as the Einstein Telescope.Once again, the similarities between the different curves show that the star's rotation can be degenerate with the fluid flow in the frequency domain.
The computational experiment performed with model C010s17 was repeated for the more compact model C019s08.We found that for almost any value of the parameter Λ that we used the star was highly destabilized.For a select set of values (e.g.A = 0.3, 0.4) where the star survived, its rotation rate was unaffected and the its gravitational wave content was not reduced (actually the gravitational waves became more complicated due to the induced oscillations).Thus we were unable to create Dedekind-type of flows for these highly relativistic and compressible objects.One way to probably improve our treatment is to use the relativistic magnitude of the velocity γ ij v i v j instead of the Newtonian one used here.We plan to examine this problem in the future.

Conclusions
We constructed constant rest-mass sequences of triaxial uniformly rotating neutron stars with a compressible equation of state in general relativity.We examined the stability of the most triaxial members of these sequences founding them stable against radial and nonaxisymmetric perturbations.These quasiequilibria are the analogs of the incompressible Jacobi ellipsoids in Newtonian gravity.Jacobi ellipsoids are congruent to their Dedekind counterparts with no internal motion.A Jacobi ellipsoid with angular velocity Ω has the same principal axes as the Dedekind ellipsoid with vorticity ζ = ( Rx Ry + Ry Rx )Ω.In general relativity this picture may be different.Here we simulated a Dedekind-type of flow in an Jacobi-type relativistic figure of quasiequilibrium.We found that for small compactness (0.1) the triaxial neutron star evolves to a Riemann-S-type of ellipsoid with minimal rotation and gravitational wave emission.On the other hand, our high compactness (0.19) triaxial model although similarly dynamically stable, was unstable to almost all Dedekind-type of flows that we tried.An important product of this investigation is the influence of a hydrodynamical fluid flow on the generation of gravitational waves and therefore the parameter estimation of a certain physical system.

Acknowledgements
The simulation and analysis tasks for this manuscript were performed at NCAR-Wyoming Supercomputing Center (NWSC) via WRAP allocation WYOM0144 "Numerical simulations of rotating neutron stars with Einstein Toolkit".This work used Stampede2 at TACC through allocation PHY160053 from the Extreme Science and Engineering Discovery Environment (XSEDE), which was supported by National Science Foundation grant number #1548562 [68].This work used Stampede2 at TACC through allocation PHY160053 from the Advanced Cyberinfrastructure Coordination Ecosystem: Services & Support (ACCESS) program, which is supported by National Science Foundation grants #2138259, #2138286, #2138307, #2137603, and #2138296 [69].Plots were produced using matplotlib [70,71] and kuibit [72].

Figure 1 :
Figure 1: Left panel: Mass versus rest-mass density for the spherical (black line) and mass-shedding (red line) limits.Also plotted are sequences of constant rest-mass (green and orange lines) for compactions (M/R) s = 0.19 and 0.1.Magenta and blue circles denote triaxial models.Right panel: T /W versus the eccentricity e for sequences of Mclaurin (ML-orange, ML-green) and Jacobi (JB-blue, JB-magenta) type ellipsoids for compactions (M/R) s = 0.1, 0.19.Solid magenta and blue circles are models C010s17 and C019s08 respectively.

r a = 0 :
Radial coordinate where the radial grids start.r b = 10 6 : Radial coordinate where the radial grids end.r c = 1.25 : Radial coordinate between r a and r b where the radial grid spacing changes.N r = 384 : Number of intervals ∆r i in r ∈ [r a , r b ].N f r = 128 : Number of intervals ∆r i in r ∈ [r a , 1].N m r = 160 : Number of intervals ∆r i in r ∈ [r a , r c ]. N θ = 96 : Number of intervals ∆θ j in θ ∈ [0, π].N ϕ = 96 : Number of intervals ∆ϕ k in ϕ ∈ [0, 2π].L = 12 : Order of included multipoles.

Figure 2 :
Figure 2: Behavior of maximum density for models C010s17 and C019s08 under a 5% pressure depletion.

Table 2 :
Summary of grid parameters used by cocal to produce the triaxial models.Note that N f r = 128 is the number of points across the largest star radius.