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

Abstract

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.

1. 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 are 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 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 (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 $\beta :=T/\left|W\right|$, where T is the rotational kinetic energy and $\left|W\right|$ 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 $\beta ={\beta }_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 $\beta ={\beta }_d>{\beta }_s$, the star becomes dynamically unstable to a bar mode deformation. This instability emerges regardless of any possible dissipation, and its growth is set by the dynamical timescale. For incompressible stars in Newtonian gravity, ${\beta }_s^{\mathrm{Newt}}=0.1375$ and ${\beta }_d^{\mathrm{Newt}}=0.2738$ [19]. Although the values of $\beta$ at these critical points can change in general relativity, with compressible equations of state and differential rotation, the overall idea (the existence of distinct secular and dynamical instability points) remains (Nonaxisymmetric instabilities for values of $\beta$ as low as 0.01 have also been found [20,21]. These so-called shear instabilities depend on $\beta$ and the degree of differential rotation [22]).

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 a 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\ne R_z$. At the point of secular instability, when $\beta =0.1375$, two families of triaxial ($R_x\ne R_y\ne R_z$) solutions emerge: (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 (this does not mean that the evolution along the Dedekind sequence does not produce gravitational waves).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{\rho }^{\Gamma }$, the triaxial sequence exists only if $\Gamma ⪆2.24$ [31]. In general relativity, the critical adiabatic index does not change significantly but slightly increases to ∼2.4 [32,33,34,35]. At the same time, the critical parameter ${\beta }_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 $\beta <0.1375$ [37,38] so that the two instabilities no longer occur at the same value of ${\beta }_s$.

Sequences of triaxial solutions in general relativity were 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 were 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\%$ of the one coming from the original triaxial neutron star as it evolves towards axisymmetry, and this highlights the importance of the fluid flow in accurate gravitational wave analysis.

Here, we employ geometric units in which $G=c=M_{\odot }=1$, unless stated otherwise. Greek indices denote spacetime dimensions (0, 1, 2, 3), while Latin indices denote spatial ones (1, 2, 3).

2. 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 important features of our initial data models and their evolutions.

2.1. Initial Data

We construct uniformly rotating triaxial neutron stars having angular velocity $\Omega$ and velocity with respect to the inertial frame $v^i=\Omega {\varphi }^i=\Omega (-y,x,0)$. The fluid’s 4-velocity can be written as

$$u^{\alpha }=u^t{k}^{\alpha }=u^t(t^{\alpha }+v^{\alpha }),$$

(1)

where $u^t$ is a scalar. The spacetime of the rotating star possesses a helical Killing vector, $k^{\alpha }$, where

$$k^{\alpha }=t^{\alpha }+\Omega {\varphi }^{\alpha },$$

(2)

with the fluid variables being Lie-dragged along $k^{\alpha }$,

$${\mathcal{L}}_{\mathbf{k}}\left(h{u}_{\alpha }\right)={\mathcal{L}}_{\mathbf{k}}\rho ={\mathcal{L}}_{\mathbf{k}}s=0.$$

(3)

Here, $\rho ,h,s$ are the rest mass density, enthalpy, and the entropy per unit rest mass. We have $\rho 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 $\Gamma =4$. For the polytropic constant, we choose $k=1$ in $G=c=M_{\odot }=1$. Similar to [41,43], the value of $\Gamma$ used is simply to prove a point of principle, rather than to address physical EOS parameters.

The models are 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 [50,51,52,53,54]. Therefore the 3-metric is ${\gamma }_{ij}={\psi }^4{f}_{ij}$, where $\psi$ is the conformal factor and $f_{ij}$ the flat metric. The unknown gravitational variables in the 3 + 1 formulation are the lapse $\alpha$, the shift ${\beta }^i$, and the conformal factor $\psi$. In the cocal code, the full system of equations (waveless 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 the following: (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:=\frac{1}{2}J\Omega$ (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 then used to compute the rotation parameter $\beta =T/\left|W\right|$. We also define the moment of inertia as $I:=J/\Omega$.

To measure the accuracy of the initial data, we use two diagnostics: The first one is the difference between the Komar and ADM mass [40],

$$\delta M=\frac{|M_{\mathrm{K}}-M|}{M_{\mathrm{K}}}.$$

(4)

For stationary and asymptotically flat spacetimes, $M_{\mathrm{K}}=M$ [56]. The second diagnostic is the relativistic virial equation [57]. In our calculations, both diagnostics are $\sim O\left({10}^{-4}\right)$.

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$) colors in Figure 1, are the analogues of the Maclaurin spheroids in Newtonian gravity [19]. In the left panel of Figure 1, the mass is plotted against the central density, while in the right panel we plot the the rotation parameter $\beta$ versus the eccentricity $e:=\sqrt{1-{({\overline{R}}_z/{\overline{R}}_x)}^2}$ with respect to the proper radii. Note that the mass and the density can by rescaled to any number using the polytropic constant k (in geometric units, $k^{n/2}$, where $\Gamma =1+1/n$ and n is the polytropic index, has units of length); hence, the axes in the left panel are normalized accordingly. In the left panel of Figure 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 a 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 ($\beta$), a second branch of solutions appears. These are triaxial solutions ($R_x\ne R_y\ne 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 Figure 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/\left|W\right|$ 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 $\Omega$, 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 and hence less total energy $T+W=M-M_p$ ($T+W$ is more negative for the triaxial solution). Thus, it is the preferred figure of equilibrium.

From the right panel of Figure 1, we notice that the bifurcation point happens at larger $\beta$ 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

$${\beta }_s={\beta }_s^{\mathrm{Newt}}+0.126{\left(\frac{M}{R}\right)}_s\left(1+{\left(\frac{M}{R}\right)}_s\right)$$

(5)

where ${\beta }_s^{\mathrm{Newt}}=0.1375$ at eccentricity $e_s^{\mathrm{Newt}}=0.8127$. Equation (5) predicts that ${\beta }_s=0.15$ at ${(M/R)}_s=0.1$, while for ${(M/R)}_s=0.19$, it is ${\beta }_s=0.166$, which is in broad agreement with the right panel of Figure 1. Notice also that the IWM formulation slightly overestimates ${\beta }_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 Figure 1 as blue and magenta dots. They constitute the most triaxial solutions along the corresponding sequences of constant rest mass. In

Table 1. Initial data models C010s17 and C019s08 used in this work. Here, $\rho$, $R_i$, $e=\sqrt{1-{({\overline{R}}_z/{\overline{R}}_x)}^2}$, $\Omega$, M, $M_0$, ${(M/R)}_s$, $T/\left|W\right|$, and 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, we use the fact that $1=1.477\mathrm{km}=4.927\mathsf{\mu }\mathrm{s}=1.989\times {10}^{33}\mathrm{g}$.

Model	$\mathit{\rho }$	${\mathit{R}}_{\mathit{x}}$	${\mathit{R}}_{\mathit{z}}/{\mathit{R}}_{\mathit{x}}$	${\mathit{R}}_{\mathit{y}}/{\mathit{R}}_{\mathit{x}}$	e	$\mathbf{\Omega }\mathit{M}$
C010s17	$0.3430$	$0.4421$	$0.4444$	$0.6875$	$0.8918$	$0.01808$
C019s08	$0.4583$	$0.4476$	$0.4619$	$0.7813$	$0.8778$	$0.05098$
Model	$\mathit{M}$	${\mathit{M}}_{\mathbf{0}}$	$\mathit{J}/{\mathit{M}}^{\mathbf{2}}$	${(\mathit{M}/\mathit{R})}_{\mathbf{s}}$	$\mathit{T}/\left|\mathbf{W}\right|$	$\mathit{I}(\times {\mathbf{10}}^{-\mathbf{3}})$
C010s17	$0.03042$	$0.03193$	$1.117$	$0.1$	$0.1547$	$1.740$
C019s08	$0.06888$	$0.07578$	$0.9011$	$0.19$	$0.1676$	$5.781$

, these two solutions are dubbed as C010s17 (magenta corresponds to compactness $0.1$) and C019s08 (blue corresponds to compactness $0.19$).

We employed the single-star module of the cocal code to compute the quasiequilibrium solutions of this work. This module uses the KEH method [58,59] on a single spherical patch $(r,\theta ,\varphi )$ with $r\in [r_a,r_b]$, $\theta \in [0,\pi ]$, and $\varphi \in [0,2\pi ]$, where $r_a=0$, $r_b=O\left({10}^{6}M\right)$ (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. Summary of grid parameters used by cocal to produce the triaxial models. Note that $N_r^{\mathrm{f}}=128$ is the number of points across the largest star radius.

$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 $\Delta r_i$ in $r\in [r_a,r_b]$.
$N_r^{\mathrm{f}}=128$	: Number of intervals $\Delta r_i$ in $r\in [r_a,1]$.
$N_r^{\mathrm{m}}=160$	: Number of intervals $\Delta r_i$ in $r\in [r_a,r_c]$.
$N_{\theta }=96$	: Number of intervals $\Delta {\theta }_j$ in $\theta \in [0,\pi ]$.
$N_{\varphi }=96$	: Number of intervals $\Delta {\varphi }_k$ in $\varphi \in [0,2\pi ]$.
$L=12$	: Order of included multipoles.

, along with the specific values used here.

2.2. 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 ${\tilde{\gamma }}_{ij}$, the conformal factor $\varphi$, (${\gamma }_{ij}=e^{4\varphi }{\tilde{\gamma }}_{ij}$), the conformally rescaled, tracefree part of the extrinsic curvature, ${\tilde{A}}_{ij}$, the trace of the extrinsic curvature, K, and three auxiliary variables ${\tilde{\Gamma }}^i=-{\partial }_j{\tilde{\gamma }}^{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 2nd-order “Gamma-driver” for the shift [63].

The equations of hydrodynamics are solved in conservation-law form, adopting high-resolution shock-capturing methods [64]. The primitive, hydrodynamic matter variables are the rest mass density, $\rho$, the pressure p, and the three-coordinate velocity $v^i=u^i/u^0$. The enthalpy is written as $h=1+e+p/\rho$, and therefore the stress energy tensor is $T_{\alpha \beta }=\rho h{u}_{\alpha }{u}_{\beta }+p{g}_{\alpha \beta }$. Here, e is the specific internal energy (this should not be confused with the eccentricity in Table 1).

The grid structure used in these evolutions is summarized in

Table 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 at the coarser level.

Model	Grid Hierarchy	$\mathit{dx}$	N
C010s17	$\{8.49,4.24,2.12,1.06,0.531\}$	$8.84\times {10}^{-2}$	80
C019s08	$\{8.59,4.30,2.15,1.07,0.537\}$	$8.95\times {10}^{-2}$	80

. Typically, we use five refinement levels, with the innermost level’s half-side length being approximately ∼1.5 times larger than the radius of the star in the initial data ($R_x$). We use $240\times 240\times 240$ cells for the innermost refinement level, which means that we have approximately 160 points across the neutron star’s 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 $\Delta x\sim 5.53\times {10}^{-3}$ (C010s17) and $\Delta x\sim 5.60\times {10}^{-3}$ (C019s08). This number of points is necessary in order to have accurate evolutions of such stiff equations of state ($\Gamma =4$), which present a challenge for any evolution code.

3. 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

$$\frac{t_d}{P}=\frac{t_d}{M}\frac{M}{P}\sim \frac{1}{\Omega M}\frac{M}{P}\sim 0.16.$$

(6)

In Figure 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 Figure 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

$$v^i=\left(\left(\frac{R_x}{R_y}\Lambda -\Omega \right)y,\left(-\frac{R_y}{R_x}\Lambda +\Omega \right)x,0\right),$$

(7)

where $\Omega$ is the angular velocity of the ellipsoid and $\Lambda$ the angular frequency of the internal fluid circulation [19,27]. When there is no internal fluid circulation ($\Lambda =0$), the fluid velocity describes the velocity field of a Jacobi ellipsoid with vorticity $\zeta =-2\Omega$, while when there is no angular velocity ($\Omega =0$) the fluid velocity describes the velocity field of a Dedekind ellipsoid with vorticity $\zeta =-(\frac{R_x}{R_y}+\frac{R_y}{R_x})\Omega$.

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 with the corresponding one of Dedekind ellipsoids. By substituting in Equation (7) $\Omega =0$ and $\Lambda =\Omega$, 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 $|\Lambda |{({(y{R}_x/R_y)}^2+{(x{R}_y/R_x)}^2)}^{1/2}$, where $\Lambda$ 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. Although the constraint equations are not solved after the new velocity profile is imposed, they remain sufficiently small ($\lesssim {10}^{-5}$) at least initially.

Setting $A=\Lambda /\Omega$, 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 of Figure 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 Figure 3, we plot the $m=2$ nonaxisymmetric mode amplitude for all four models (C010s17 and its perturbed models)

$$C_m=\int \rho u^t{e}^{im\varphi }\sqrt{-g}{d}^{3}x,$$

(8)

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 preserves its shape (as an equilibrium Dedekind ellipsoid) would have a constant $C_2$ at 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 linear fits with dashed lines. 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\gtrsim 10$) is C010s17-A070.

In order to appreciate the overall motion of these ellipsoidal figures, we plot in Figure 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 in the right column we plot model C010s17-A075. As can be seen from the first five instances ($t/P\lesssim 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\sim 30P$, all models remain nonaxisymmetric (although less than at $t=0$) and continue to oscillate mildly.

In the left panel of Figure 5, we plot the gravitational wave strain $h_{\times }$ for the C010s17 models following [65]. It is obvious that the Dedekind-like velocity flow decreases the gravitational wave signature of the triaxial figures by less than half, even at early times. Consistent with Figure 4 and the left panel of Figure 3, we see that the model with the least amount of gravitational wave content is C010s17-A070 (orange line), whose strain exhibits a decrease of ∼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 Figure 5, we plot the power spectrum of the C010s17 models scaled for a $1.4M_{\odot }$ triaxial neutron star mass at $50\mathrm{Mpc}$, along with the noise curves for Advanced LIGO’s design sensitivity [66] and the ET-B configuration of the Einstein Telescope [67]. The peak frequency of the unperturbed C010s17 model (blue curve) at ∼840$\mathrm{Hz}$ is consistent with its orbital angular frequency ($\Omega /\pi$) and, in principle, is 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 $\Lambda$ 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 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 flows for these highly relativistic and compressible objects. One way to probably improve our treatment is to use the relativistic magnitude of the velocity $\sqrt{{\gamma }_{ij}{v}^i{v}^j}$ instead of the Newtonian one used here. We plan to examine this problem in the future.

4. 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, finding 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 $\Omega$ has the same principal axes as the Dedekind ellipsoid with vorticity $\zeta =(\frac{R_x}{R_y}+\frac{R_y}{R_x})\Omega$. 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 flows that we tried. This does not mean that general relativistic analogues of Dedekind ellipsoids do not exist, but a self-consistent calculation is necessary to address their existence and stability. 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.