Gravitational wave content and stability of uniformly, rotating, triaxial neutron stars in general relativity

TL;DR

This study simulates uniformly rotating, triaxially deformed neutron stars in general relativity to assess their stability and gravitational wave emission, finding they are dynamically stable and evolve towards axisymmetry.

Contribution

It provides the first detailed evolution analysis of triaxial neutron stars in general relativity, including stability and gravitational wave predictions.

Findings

All models are dynamically stable.

Stars emit gravitational waves with strain about one-tenth that of binary mergers.

Stars evolve towards axisymmetry while maintaining rotation and energy profiles.

Abstract

Targets for ground-based gravitational wave interferometers include continuous, quasiperiodic sources of gravitational radiation, such as isolated, spinning neutron stars. In this work we perform evolution simulations of uniformly rotating, triaxially deformed stars, the compressible analogues in general relativity of incompressible, Newtonian Jacobi ellipsoids. We investigate their stability and gravitational wave emission. We employ five models, both normal and supramassive, and track their evolution with different grid setups and resolutions, as well as with two different evolution codes. We find that all models are dynamically stable and produce a strain that is approximately one-tenth the average value of a merging binary system. We track their secular evolution and find that all our stars evolve towards axisymmetry, maintaining their uniform rotation, kinetic energy, and angular momentum profiles while losing their triaxiality.

Tables

Table 1. Table 1: Characteristic quantities for the maximum mass spherical solutions of the two EoSs considered in this work. First line refers to simple polytrope models G4C010, G4C025, while second line to piecewise models pwC023, pwC024, pwC026. To convert to cgs units multiply mass, density and pressure by 1.989 × 10 33 ​ g 1.989 superscript 10 33 g 1.989\times 10^{33}\ {\rm g} , 6.173 × 10 17 ​ g / cm 3 6.173 superscript 10 17 g superscript cm 3 6.173\times 10^{17}\ {\rm g/cm^{3}} , and 5.548 × 10 38 ​ g / ( cm ​ sec 2 ) 5.548 superscript 10 38 g cm superscript sec 2 5.548\times 10^{38}\ {\rm g/(cm\ sec^{2})} , respectively.

$\Gamma$	${\left(P/{\rho }_0\right)}_c$	${\rho }_c$	${\left({\rho }_0\right)}_c$	$M_0$	$M$	$M/R$
$4$	$1.334$	$0.004658$	$0.003224$	$2.882$	$2.250$	$0.3552$
$(4,2.5)$	$0.5674$	$0.006175$	$0.004536$	$1.960$	$1.657$	$0.2871$

Table 2. Table 2: Models G4C010, G4C025 have Γ = 4 Γ 4 \Gamma=4 throughout, while models pwC023, pwC024, pwC026 have ( Γ 1 , Γ 2 ) = ( 4 , 2.5 ) subscript Γ 1 subscript Γ 2 4 2.5 (\Gamma_{1},\Gamma_{2})=(4,2.5) and are supramassive. Here ρ 0 subscript 𝜌 0 \rho_{0} , ρ 𝜌 \rho , R i subscript 𝑅 𝑖 R_{i} , e z = 1 − ( R y / R x ) 2 subscript 𝑒 𝑧 1 superscript subscript 𝑅 𝑦 subscript 𝑅 𝑥 2 e_{z}=\sqrt{1-(R_{y}/R_{x})^{2}} , Ω Ω \Omega , P 𝑃 P , M 𝑀 M , M 0 subscript 𝑀 0 M_{0} , J 𝐽 J , ( M / R ) s subscript 𝑀 𝑅 𝑠 (M/R)_{s} , T / | W | 𝑇 𝑊 T/|W| , I 𝐼 I , ε z subscript 𝜀 𝑧 \varepsilon_{z} , E ˙ ˙ 𝐸 \dot{E} , J ˙ ˙ 𝐽 \dot{J} , h ℎ h , t d subscript 𝑡 𝑑 t_{d} , t s subscript 𝑡 𝑠 t_{s} are the rest-mass density, the total energy density, the coordinate radii, the eccentricity with respect to the z-axis, the angular velocity, the period, the ADM mass, the rest-mass, the ADM angular momentum, the corresponding spherical compactness, the parameter β 𝛽 \beta , the moment of inertia, the ellipticity with respect to the z-axis (Eq. 21 ), the luminosity, the angular momentum loss rate, the GW maximum amplitude, the dynamical timescale, and the secular timescale, respectively. To convert to geometric G = c = 1 𝐺 𝑐 1 G=c=1 or cgs units, use the fact that 1 = 1.477 ​ km = 4.927 ​ μ ​ s = 1.989 × 10 33 ​ g 1 1.477 km 4.927 𝜇 s 1.989 superscript 10 33 g 1=1.477\ {\rm km}=4.927\ {\rm\mu s}=1.989\times 10^{33}\ {\rm g} .

	Initial data models
	G4C010	G4C025	pwC023	pwC024	pwC026
${\rho }_0(\times {10}^{-3})$	$1.005$	$1.565$	$1.902$	$1.991$	$2.226$
$\rho (\times {10}^{-3})$	$1.019$	$1.644$	$2.065$	$2.176$	$2.477$
$R_x$	$7.677$	$7.429$	$7.774$	$7.625$	$7.266$
$R_z/R_x$	$0.4727$	$0.4957$	$0.4977$	$0.5015$	$0.5108$
$R_y/R_x$	$0.7500$	$0.9063$	$0.9219$	$0.9375$	$0.9688$
$e_z$	$0.8812$	$0.8685$	$0.8673$	$0.8652$	$0.8597$
$\Omega M$	$0.01823$	$0.08043$	$0.07850$	$0.08237$	$0.09138$
$P$ (Period)	$193.8$	$138.3$	$140.8$	$137.7$	$130.4$
$M$	$0.5623$	$1.771$	$1.760$	$1.805$	$1.896$
$M_0$	$0.5900$	$2.012$	$1.989$	$2.047$	$2.169$
$J/M^2$	$1.109$	$0.8516$	$0.8279$	$0.8202$	$0.8003$
${\left(M/R\right)}_s$	$0.1000$	$0.2500$	-	-	-
$M/R_x$	$0.07324$	$0.2383$	$0.2264$	$0.2367$	$0.2610$
$T/\left|W\right|$	$0.1543$	$0.1773$	$0.1666$	$0.1661$	$0.1633$
$I$	$10.81$	$58.77$	$57.46$	$58.51$	$59.70$
${\epsilon }_z$	$0.2320$	$0.05581$	$0.02771$	$0.0191$	$0.006200$
$\delta M(\times {10}^{-4})$	$0.8237$	$0.9893$	$1.129$	$1.063$	$1.007$
VE$(\times {10}^{-4})$	$12.13$	$5.463$	$8.047$	$7.753$	$7.546$
	Quadrupole estimates
$\dot{E}(\times {10}^{-8})$	$2.846$	$14.01$	$2.918$	$1.548$	$0.2017$
$\dot{J}(\times {10}^{-7})$	$8.778$	$30.84$	$6.540$	$3.391$	$0.4186$
$rh/M(\times {10}^{-3})$	$14.63$	$7.357$	$3.441$	$2.389$	$0.7774$
	Timescales
$t_d/M$	$50$	$10$	$10$	$10$	$10$
$t_s/M$	${10}^5$	${10}^5$	${10}^6$	${10}^6$	${10}^7$

Table 3. Table 3: Summary of grid parameters used by cocal to produce the five models. Note that N r f = 128 superscript subscript 𝑁 𝑟 f 128 {N_{r}^{\rm 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^f=128$	:	Number of intervals $\Delta r_i$ in $r\in [r_a,1]$.
$N_r^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.

Table 4. Table 4: Grid parameters used for the evolution of each model. Parameter N 𝑁 N corresponds to the number of points used to cover the largest radius of the star. Parameter d ​ x 𝑑 𝑥 dx is the step interval in the coarser level. To convert to physical units multiply by 1 = 1.477 ​ km . 1 1.477 km 1=1.477\ {\rm km}.

Model	$x_{\min }$	$x_{\max }$	$y_{\min }$	$y_{\max }$	$z_{\min }$	$z_{\max }$	Grid hierarchy	$dx$	$N$
G4C010	$-304$	$304$	$-304$	$304$	$0$	$304$	$\{9.5,19.0,38.0,76.0,152.0,304.0\}$	$2.5\overline{3}$	$96$
G4C025	$-304$	$304$	$-304$	$304$	$0$	$304$	$\{9.5,19.0,38.0,76.0,152.0,304.0\}$	$2.5\overline{3}$	$93$
pwC023	$-304$	$304$	$-304$	$304$	$0$	$304$	$\{9.5,19.0,38.0,76.0,152.0,304.0\}$	$2.5\overline{3}$	$98$
pwC024	$-304$	$304$	$-304$	$304$	$0$	$304$	$\{9.5,19.0,38.0,76.0,152.0,304.0\}$	$2.5\overline{3}$	$96$
pwC026	$-288$	$288$	$-288$	$288$	$0$	$288$	$\{9.0,18.0,36.0,72.0,144.0,288.0\}$	$2.4$	$96$
	$-288$	$288$	$-288$	$288$	$0$	$288$	$\{9.0,18.0,36.0,72.0,144.0,288.0\}$	$2.88$	$80$
	$-1152$	$1152$	$-1152$	$1152$	$0$	$1152$	$\{9.0,18.0,36.0,72.0,144.0,288.0,576.0\}$	$2.88$	$96$

Equations

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

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

$$u^{\alpha }=u^t{k}^{\alpha }=u^t(t^{\alpha }+v^{\alpha })=\alpha u^t(n^{\alpha }+U^{\alpha }).$$

$$\frac{h}{u^t}=\mathcal{E},$$

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

$$I^{ij}:={\int }_{{\Sigma }_t}{\rho }_0{u}^{\alpha }{x}^i{x}^{j}d{S}_{\alpha }$$

$$P=P_{\mathrm{cold}}+P_{\mathrm{th}}=P_{\mathrm{cold}}+({\Gamma }_{\mathrm{th}}-1){\rho }_0(ϵ-{ϵ}_{\mathrm{cold}})$$

$${ϵ}_{\mathrm{cold}}=-\int P_{\mathrm{cold}}d(1/{\rho }_0)=\frac{k}{\Gamma -1}{\rho }_0^{\Gamma -1}+const..$$

$$\frac{t_d}{M}\sim \frac{1}{\Omega M}\sim {\left(\frac{M}{R}\right)}^{-3/2},$$

$$\dot{E}\propto \dot{J}\propto (rh{)}^2.$$

$$\frac{{\tau }_{\mathrm{bar}}^{\mathrm{GW}}}{M}\sim 2\times 10^{-3}{\left(\frac{M}{R}\right)}^5(M\Omega {)}^{-6}(\beta -{\beta }_{\sec }{)}^{-5},$$

$$\frac{{\tau }_{\mathrm{r}}^{\mathrm{GW}}}{M}\sim 10{\left(\frac{M}{R}\right)}^4(M\Omega {)}^{-6},$$

$$h_{ij}(t,x^a)=\frac{2}{r}{\left[\frac{d^2{to0.0ptto5.18054pt\text{hss}--\text{hss}\text{hss}I}_{ij}^{TT}}{d{t}^2}\right]}_{\mathrm{ret}}$$

$$\frac{dE}{dt}=\frac{1}{5}⟨{\stackrel{...}{to0.0ptto5.18054pt\text{hss}--\text{hss}\text{hss}I}}_{ij}{\stackrel{...}{to0.0ptto5.18054pt\text{hss}--\text{hss}\text{hss}I}}_{ij}⟩,\frac{d{J}^i}{dt}=\frac{2}{5}{ϵ}^{ijk}⟨{\ddot{to0.0ptto5.18054pt\text{hss}--\text{hss}\text{hss}I}}_{ja}{\stackrel{...}{to0.0ptto5.18054pt\text{hss}--\text{hss}\text{hss}I}}_{ka}⟩,$$

$$\frac{d}{dt}{I}_{ij}={\int }_{\Sigma }{\rho }_0{u}^{\alpha }(x_i{v}_j+x_j{v}_i)d{S}_{\alpha },$$

$$\frac{d}{dt}{\int }_{V_t}{\rho }_{\ast }QdV={\int }_{V_t}{\rho }_{\ast }\frac{DQ}{Dt}dV,$$

$$\frac{d}{dt}{\int }_{{\Sigma }_t}Q{\rho }_0{u}^{\alpha }d{S}_{\alpha }={\int }_{{\Sigma }_t}{\mathcal{L}}_{\mathbf{k}}Q{\rho }_0{u}^{\alpha }d{S}_{\alpha }.$$

$$f^{\prime }(0)$$

$$\frac{D{x}^A}{Dt}$$

$$\frac{D{v}^A}{Dt}$$

$$\frac{D{a}^A}{Dt}$$

$${\dot{I}}^{ij}(0)$$

$${\ddot{I}}^{ij}(0)$$

$${\stackrel{...}{I}}^{ij}(0)$$

$$[h_{AB}]$$

$$h_{(+,\times )}=\frac{4{\Omega }^2}{r}(I_1-I_2)(\cos (2\Omega t),\sin (2\Omega t)),$$

$$|\dot{E}|=\frac{32}{5}(I_1-I_2{)}^2{\Omega }^6,$$

$$|\dot{J}|=\frac{32}{5}(I_1-I_2{)}^2{\Omega }^5.$$

$${\epsilon }_z:=\frac{|I_{11}-I_{22}|}{I_{11}+I_{22}}.$$

Full text

I Introduction

The discovery of gravitational waves abbott2016a from a binary black-hole system was a triumph that initiated a new era in astronomy and astrophysics. Although the prime candidates for the ground-based interferometers are binary systems, gravitational waves from isolated neutron stars also can be detected and help reveal the nature of these objects. Out of the $\sim 2500$ currently known pulsars in our Galaxy, approximately $90\%$ are isolated. Many of these single rotating stars may be promising sources of gravitational waves EatH ; Aasietal2013 ; Papaetal2016 .

A single neutron star can become an emitter of gravitational waves (GWs) as long as it has a nonspherical time changing quadrupole moment. The lack of symmetry can arise in various scenarios CA ; A2010 ; K03 . For example a pulsar can have a “small mountain” that could develop following a starquake in the NS ST83 ; HJA2006 , or it can exhibit different kinds of nonspherical oscillations FS14 . Another possibility is binary neutron star mergers, which are themselves prime candidates for the production of gravitational radiation. When the two component stars do not have large masses the remnant may not undergo “prompt” collapse, but instead form a hypermassive star and undergo “delayed collapse”, or form a spinning neutron star that is dynamically and secularly stable BR2016 . At formation such remnants may be nonaxisymmetric and strong GW emitters. A third scenario arises in gravitational stellar collapse, where the bouncing core can be rotating fast enough so that nonaxisymmetric instabilities set in and deform the star into an ellipsoid LS1995 . Fallback accretion onto newly born magnetars also supports the existence of triaxial deformations and the efficient production of gravitational waves PO2011 .

Despite the enormous amount of work done in the field of rotating stars FS2013 ; PS2016 full general relativistic (GR) numerical simulations that investigate the stability and accurately quantify the GW signature of single, uniformly rotating, triaxial stars have not been performed. One of the reasons is the scarcity of accurate initial models needed to study their evolution. Typically these objects are probed in the context of binary mergers or collapse scenarios, which involve a substantial amount of computational resources and make difficult a systematic parameter study. In these cases one typically ends up with a differentially rotating object while for single, isolated neutron stars one is often interested in uniformly rotating stars, the GR analogues of Jacobi ellipsoids in Newtonian theory. Such solutions have been obtained for the first time by Nozawa in his PhD thesis N1997 , where he allowed for azimuthal dependence in the spacetime metric, but restricted it to an axisymmetric form. Using a different method, triaxial quasiequilibrium models have been computed without such a restriction in the conformal flatness approximation HMSU2008 and in the waveless approximation UTGHSTY2016 as part of the cocal code.

The ab initio computation of such nonaxisymmetric objects presents a number of challenges. First, these objects are not stationary equilibria, since they emit GWs, and therefore an approximate scheme has to be applied in order to find quasistationary solutions. This choice has be compatible with the fact that the radiated energy within one rotational period is much smaller than the binding energy of the star. Second, such models are known to exist only for stiff equations of state. If we assume a polytropic law $P=k\rho_{0}^{\Gamma}$, where $\rho_{0}$ is the rest-mass density and $k,\ \Gamma$ are constants, then $\Gamma$ needs to be larger than $2.24$ in the Newtonian limit James64 . For softer equations of state mass shedding appears at lower angular velocity than the one needed to reach the triaxial state. GR increases the critical value of the polytropic index by a small amount (to $\Gamma\sim 2.8$) BFG1996 . Third, uniformly rotating, nonaxisymmetric solutions exist only for high spin rates, i.e. $\beta:=T/|W|$ larger than $0.14$ in the Newtonian case Ch69 . Here $T$ is the kinetic energy and $W$ the gravitational potential energy. In GR this critical value is higher SZ1996 ; BFG1998 ; GV2002 ; Saijo2006 ; RG2002 ; SL1996 ; YE1997 . The combination of the above characteristics imply that an evolutionary follow-up will also be nontrivial, since the first challenge described above will imply the presence of junk initial radiation, which must be controlled, while the second and third challenges require higher resolution than for slowly rotating stars. Since the GW timescale to radiate the rotational energy is $t_{\rm GW}/M\gtrsim(M/R)^{-4}$ only highly compact objects can be evolved to their endpoint state, while lower compacteness stars can be studied only partially. High compactness requires higher resolution, which increases the computational demands even more.

The dynamical stability of the quasiequilibrium solutions obtained in HMSU2008 ; UTGHSTY2016 is not yet known. If these objects are dynamically unstable, do they undergo prompt collapse to black holes, or do they evolve to significantly different, stable, axisymmetric equilibria by rearranging their mass and/or angular momentum profiles? If they are dynamically stable, their secular fate is still unknown. Being nonaxisymmetric and rotating they will generate GWs, which will radiate both energy and angular momentum. Will this lead to delayed collapse to a black hole, or will it lead to the formation of a Dedekind-like configuration, or something less exotic?

In BSS00 ; SBS00 the dynamical stability of axisymmetric, differentially rotating stars (even including an initial perturbation) has been studied numerically in GR and it was found that they are stable against quasiradial collapse and bar-mode formation for sufficiently small $\beta$. GR enhances the dynamical bar-mode formation since the critical value for $\beta=\beta_{\rm dyn}$ above which the stars become dynamically unstable was found to be $\sim 0.24$, slightly less than the corresponding Newtonian value $0.27$ for incompresible Maclaurin spheroids. A precise determination of the threshold for the dynamical instability, the effects of stellar compactness on that, as well as the timescale of the persistence of the bar deformation have been studied in BPMR07 ; MBPR07 . In SKE02 ; KE03 ; SKE03 linear stability analysis and simulations have been performed to analyze the occurrence of the dynamical instability against nonaxisymmetric bar mode deformation for differentially rotating stars. It was found that when differential rotation is high, the stars are dynamically unstable even when $\beta$ was of order of $0.01$. This dynamical instability does not create spiral arms Centrella01 ; Saijo03 ; Ott05 ; Ou06 ; Paschalidis15 ; East16 or fragmentations, but drives the star into a quasi-stationary ellipsoid that emitts GWs.

The secular bar-mode instability induced by gravitational radiation with a polytropic ($\Gamma=2$) equation of state (EoS) in the 2.5 post-Newtonian framework for rapidly rotating stars with $\beta\sim 0.2-0.25$ has been investigated in SK04 111The critical $\beta$ for instability in Newtonian Maclaurin spheroids is $\beta_{\rm sec}\sim 0.14$, but decreases in GR as the compaction increases SF98 . They tracked the evolution of the bar-mode up until the final object was a deformed ellipsoid which was still emitting GWs (therefore was not a Dedekind-like star). At the same time the nonlinear developent of the secular bar-mode instability using a stiffer EoS ($\Gamma=3$) and similarly including post-Newtonian terms for the gravitational radiation reaction was investigated in OTL04 . Although they were able to reach a “Dedekind-like” state, this was destroyed after ten dynamical times. According to the authors the reason could be either the nonlinear coupling of various oscillatory modes in the star, or an “elliptic flow” instability which manifests itself when the fluid flow is forced along elliptic streamlines.

In a previous work UTea2016 we computed for the first time triaxial supramassive neutron stars (uniformly rotating models with rest-mass higher than the maximum rest-mass of a nonrotating star, but lower than the maximum rest-mass when allowing for maximal uniform rotation), by using a piecewise polytropic EoS. In this work we perform the first evolutions of such stars and try to investigate their stability and gravitational wave content. Following UTea2016 we carefully construct five such models: two normal ones (uniformly rotating but not supramassive) with compactions $0.1$ and $0.25$222These are the corresponding compactenesses of the spherical solutions. adopting a stiff, $\Gamma=4$, EoS, and three supramassive models with compactions $0.23,\ 0.24,\ 0.26$ adopting a two-piece polytrope that has a soft core. Although these EoSs are rather extreme, our goal is to prove a matter of principle rather than focus on realistic EoSs. For a single polytrope a stiffer EoS can sustain a larger triaxial deformation, and hence the maximum mass of the triaxial star relative to that of the spherical star is expected to be larger. However, for the two-piece polytropic EoS, the maximum mass of the triaxial star relative to the spherical counterpart increases, even though the overall averaged stiffness of the EoS is softer. If the mass difference between the maximum axisymmetric and triaxial solutions is $\sim 10\%$ or less, then that implies that the EoS of high density matter becomes substantially softer in the core of neutron stars UTea2016 .

We were able to follow the evolution of these objects for more than twenty rotation periods, proving that they are dynamically stable. After an initial short period of time where junk radiation in the initial data propagates away, the neutron star evolves along quasiequilibrium states that satisfy the first law, $dM=\Omega dJ$. Along this trajectory the orbital angular velocity remains constant inside the neutron star, whose triaxial shape evolves toward axisymmetry. During this period the GW amplitude decreases significantly, especially in the highly compact models. The question that arises is: are we probing the secular fate of the stars or is this clear monotonic amplitude decrease an artifact of numerical dissipation.

We do not think that the decrease of the GW amplitude is due to numerical viscosity. We performed a resolution study which did not alter the main description above. We discuss the trigger for the declining amplitude below.

If our models are imagined to sample bar-mode perturbations of an axisymmetric configuration with $\beta>\beta_{\rm sec}$ then according to well-known results SF98 , our stars should be secularly unstable. We weren’t able to find any growth of a bar-mode. As in DLSS04 , where evolutions of models with $\beta$ larger than $\beta_{\rm sec}$ with an initial bar-mode perturbation were performed, we find the decay of the initial perturbation.

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

II Methods and physical parameters

The numerical methods used here are those implemented in the cocal and Illinois GRMHD codes, and have been described in great detail in our previous works TU07 ; UT12 ; UTG12 ; TU13 ; TU15 ; TMGRU16 ; EFLSTB08 ; ELS10 ; EPLS12 ; PLES10 ; PLES11 , so we only summarize the most important features here.

II.1 Initial data

Our initial rotating star spacetimes posses a helical Killing vector, $k^{\alpha}$,

[TABLE]

whereby the fluid variables are Lie-dragged along $k^{\alpha}$,

[TABLE]

Here $u^{\alpha}$ is the 4-velocity of the fluid, $\rho_{0},\ h,\ s$ are the rest-mass density, enthalpy, and the entropy per unit rest-mass. We have $\rho_{0}h=\rho+P$, where $\rho$ is the total energy density and $P$ is the pressure. The 4-velocity of the fluid will be along the helical Killing vector, i.e. there is a scalar $u^{t}$ such that

[TABLE]

In the above, $v^{i}=\Omega\phi^{i}=\Omega(-y,x,0)$ is the velocity with respect to the inertial frame, while $U^{\alpha}$ is the spatial velocity with respect to normal observers (those with 4-velocity $n^{\alpha}$). In the last equality, $\alpha$ is the lapse function, that normalizes the normal vector to the spacelike hypersurfaces which foliate the spacetime, $n_{\alpha}=-\alpha\nabla_{\alpha}t$.

For a perfect gas stress-energy tensor and an isentropic initial configuration the equations of motion yield a first integral,

[TABLE]

where $\mathcal{E}$ is a constant. The two constants that appear in our equations $\{\Omega,\mathcal{E}\}$ are determined via an iterative scheme. For the gravitational fields we use the Isenberg-Wilson-Mathews (IWM) approximation I1980 ; WM1989 which assumes a flat conformal metric and maximal slicing. The resulting five elliptic equations are solved together with Eq. (4) and coupled to a piecewise EoS as described in HMSU2008 ; UTGHSTY2016 .

A number of diagnostics are used to describe the initial solutions and explicit formulae are given in the appendix of UTGHSTY2016 and will not be repeated here. Since the IWM formulation is used, we have that $\gamma_{ij}=\psi^{4}f_{ij}$, where $\psi$ is the conformal factor and $f_{ij}$ the flat metric in spherical coordinates. The angular momentum $J=J_{\rm ADM}$ [where $J_{\rm ADM}$ is the Arnowitt-Deser-Misner (ADM) angular momentum] is computed via a surface integral at infinity or a volume integral over the spacelike hypersurface. The kinetic energy 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 the gravitational potential energy is defined as $W:=M_{\rm ADM}-M_{\rm P}-T$. Here $M_{\rm ADM}=M$ is the (ADM) mass and $M_{\rm P}$ is the rest-mass plus internal energy of the star (see e.g. BS10 ). These expressions are used then to compute the rotation parameter $\beta$. Also the moment of inertial is defined as $I:=J/\Omega$. As a measure of accuracy of the initial data we provide two diagnostics: The first one is the difference between the Komar and ADM mass,

[TABLE]

For stationary and asymptotically flat spacetimes $M_{\rm K}=M_{\rm ADM}$333Although for nonaxisymmetric systems the helical Killing vector (stationarity in the rotating frame) is incompatible with asymptotic flatness GS84 , one can define an approximate asymptotic region in which the gravitational wave energy is small compared with the total energy of the system. The same argument holds for the existence of the Komar mass that is associated with a timelike Killing field $t^{\alpha}$. Beig78 . The second diagnostic is the relativistic virial equation (VE) GB94 .

The initial-data gravitational wave diagnostics involve the second mass moments

[TABLE]

with $dS_{\alpha}=\nabla_{\alpha}t\sqrt{-g}d^{3}x$. In Appendix A we have derived some useful quantities such as the quadrupole approximation for the luminosity and the gravitational wave amplitude, that can be computed on a spacelike hypersurface in the presence of a helical Killing vector. However, full GW output, including the “junk” radiation inherent in the initial data, is computed in full GR as part of the integration of the field equations via the Baumgarte-Shapiro-Shibata-Nakamura (BSSN) formalism SN95 ; BS98 .

As in UTea2016 we employ the same “benchmark” EoSs. The first one is a simple $\Gamma=4$ polytrope, while the second is a piecewise-polytropic EoS with two pieces and a soft core, where $\{\Gamma_{1},\Gamma_{2}\}=\{2.5,4\}$. Characteristics of the maximum mass solutions for spherical stars using these two EoSs are reported in Table 1. The adiabatic constant $k$, is chosen so that the value of the rest-mass becomes $M_{0}=1.5$ (in units of solar mass) at the compactness $M/R=0.2$. By choosing different values of $k$ one can attain larger maximum masses. A well-known fact that relates the maximum masses of those models, is that a stiffer EoS can sustain a larger maximum mass (see below). The same result holds for the maximum masses of the axisymmetric solutions. The values of $\Gamma$ used are simply to prove a point of principle, rather than address physical EoS parameters: stiffness is necessary in order for these triaxial solutions to exist. A higher value of $\Gamma$ satisfies the necessary conditions for uniformly rotating triaxial solutions to exist, and this is the main reason behind such a choice. As discussed in UTea2016 the softening of the core enables us to compute for the first time supramassive, triaxially deformed, uniform rotating stars, without increasing further the maximum polytropic exponent. This was made possible from the following counter-intuitive fact which does not depend on the values of the specific polytropic indices: Assume a simple (any $\Gamma>2.24$) polytrope which in most cases does not support supramassive triaxial solutions. Then consider a second two piece polytropic EoS $\{\Gamma_{1},\Gamma_{2}\}$, with $\Gamma_{2}=\Gamma$ and a soft core with $\Gamma_{1}<\Gamma$. This second EoS is effectively softer than the first. Thus one expects that the piecewise EoS does not exhibit triaxial solutions with mass larger than the maximum-mass spherical solution. This was proven not to be the case UTea2016 , and if the relative difference between the maximum triaxial and axisymmetric masses is less than $10\%$444The maximum mass of triaxial solutions is always smaller than the maximum mass of axisymmetric ones. it provides strong evidence of softening in the core of the compact object.

In order to investigate the stability and gravitational wave signature of such solutions we consider five models, G4C010, G4C025, pwC023, pwC024, and pwC026 whose characteristics are reported in Table 2 555As we mentioned in the introduction all quantities reported are in $G=c=M_{\odot}=1$ units. This means that if one wants to convert mass to geometric $G=c=1$ units one has to multiply by $1=1.477\ {\rm km}$. For the angular velocity $\Omega$, one divides by $1=1.477\ {\rm km}$. Similarly to get $\Omega$ in cgs units again one divides be $1=4.927\ {\rm\mu s}$.. The last three columns are supramassive solutions while the others are normal ones. The triaxiality666Triaxiality is not used in any quantative way in this paper. It can be defined in various ways, like $R_{y}/R_{x},\ e_{z}$, or $\varepsilon_{z}$ (see Table 2) and signifies the departure from axisymmetry. In GW detection studies, triaxiality is measured by the ellipticity $\varepsilon_{z}$. Notice that the ellipticities of the models we consider here are larger than typical limits set by LIGO Abbottetal17 . However, as isolated pulsars are dim and hard to find, there could exist a population of undetected pulsars that LIGO has not probed, yet. is larger for the first column model and diminishes as we move to more compact stars. This means that the amplitude of the gravitational wave will be larger for the first model and smaller for the last one.

The models have been computed with the cocal code, a second-order finite-difference code whose methods are explained in UT12 ; UTG12 ; TU13 ; TU15 ; TMGRU16 . For single compact objects it emloys a single spherical patch $(r,\theta,\phi)$ with $r\in[r_{a},r_{b}]$, $\theta\in[0,\pi]$, and $\phi\in[0,2\pi]$, where $r_{a}=0$, $r_{b}=O(10^{6}M)$, and $M$ the total mass of the system (no compactification used). 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 3, along with the specific values used to obtain the quasi-equilibrium solutions of this work.

II.2 Evolution

For the evolution we use the Illinois GRMHD code777We do not use IllinoisGRMHD, which is the version of the code embedded in the Einstein Toolkit EPHMS ., which solves the Einstein field equations in the BSSN formalism SN95 ; BS98 ; BS10 . The code is built on the cactus cactus infrastructure and uses carpet 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 $\phi$, ($\gamma_{ij}=e^{4\phi}\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  Baker06 ; Campanelli06 ; NRAR13 , which are part of the family of gauge conditions using an advective “1 + log” slicing for the lapse, and a “Gamma-driver” for the shift  ABDKPST03 .

The equations of hydrodynamics are solved in conservation-law form adopting high-resolution shock-capturing methods EPLS12 ; ELS10 . The primitive, hydrodynamic matter variables are the rest mass density, $\rho_{0}$, the pressure $P$ and the coordinate three velocity $v^{i}=u^{i}/u^{0}$. The enthalpy is written as $h=1+\epsilon+P/\rho_{0}$, and therefore the stress energy tensor is $T_{\alpha\beta}=\rho_{0}hu_{\alpha}u_{\beta}+Pg_{\alpha\beta}$. Here, $\epsilon$ is the specific internal energy888This should not be confused with the ellipticity $\varepsilon_{z}$..

To close the system an EoS needs to be provided and for that we follow PLES10 ; PLES11 where the pressure is decomposed as a sum of a cold and a thermal part,

[TABLE]

where

[TABLE]

Here $k,\Gamma$ are the polytropic constant and exponent of the cold part (same as the initial data EoS) and $\Gamma_{\rm th}=5/3$. The constant that appears in the formula above is zero for a single EoS, but takes different values in a piecewise polytrope where one has to account for the continuity of pressure at the join between the different pieces.

The grid structure used in these evolutions is summarized in Table 4. Typically we use six refinement levels with the innermost level half-side length being approximately $\sim 1.25$ times larger than the radius of the star in the initial data ($R_{x}$). We use $240\times 240\times 120$ points for the innermost refinement level, which means that we have approximately 190 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 $\Delta x\sim 0.0791\bar{6}=117\ {\rm m}$. This number of points was necessary in order to have accurate evolutions of such stiff EoS ($\Gamma=4$) which present a challenge for any evolution code.

For the last model pwC026 we have done two extra simulations, as the compactness in this case was very high and the triaxiality very low. In this model the GW signal was very weak ($rh/M\sim 10^{-4}$) and therefore we wanted to corroborate our findings by using different resolution and box size for the outer boundary conditions. On the last two lines of Table 4 the lower resolution simulation has the same outer boundary distance ($288$) but $80$ points across the star radius, while we have also a simulation with seven refinement levels and the outer boundary at much larger distance ($1152.0$) than all other cases.

III Results

The main purpose of this work was to probe the stability of uniformly rotating, triaxial configurations and estimate their gravitational wave emission. Although remnants from neutron star mergers have an asymmetric shape and can be dynamically stable, they are differentially rotating. Is it possible for a triaxial star that rotates uniformly to be also dynamically stable? And if that is possible what is the secular fate of this configuration?

III.1 Dynamic stability

The dynamical timescale for our stars is

[TABLE]

and values range from $\sim 10$ for the most compact cases to $\sim 50$ for G4C010, the least compact model (see Table 2). We find that all of the models considered are dynamically stable. Figure 1 shows typical contour plots at $t=0$ (dashed colored lines) as well as the same contour plots after ten rotation periods (solid colored lines). The black dashed line signifies the initial data surface of the star in the xy-plane as calculated by cocal. Of particular importance are the lowest density contours at $\rho_{0}=0.0002=9.16788\times 10^{-15}\ {\rm cm^{-2}}$ (magenta colored). The choice of this particular value can be considered as one of the largest densities that follow closely the initial data profile (black dashed line). By following the evolution of this contour one can have an accurate picture of the surface of the star. After the junk radiation has propagated away the stars still retain their triaxiality. But by $t=10\ P$, all contours tend to circularize (the one of the highest density is initially circular). All these contours contract in the x-direction and expand in the y-direction. The amount of contraction/expansion diminishes as one moves towards the center of the star. Thus the star becomes more axisymmetric. After ten periods the x-axis has lost $9-8\%$ of its length. For the supramassive models this picture still holds, although since the ellipticities there are much smaller the amount of contraction/expansion is somewhat diminished. For the most supramassive model, pwC026, after ten rotation periods the decrease is $\sim 4\%$ and the object is essentially axisymmetric. While density contours are not gauge-invariant, they yield a qualitative picture that agrees with the gravitational wave signature that we discuss in the next section.

The constant angular velocity profile is well preserved (Fig. 2). The angular velocity across the x-axis (bottom panels) and the y-axis (top panels) is plotted for the G4C010 and G4C025 models. Red curves correspond to $\Omega$ after one rotation period while blue curves after ten rotation periods. Vertical brown dashed curves denote the initial data star radii, and the green curve is the Keplerian limit $\Omega_{\rm K}=(M/r^{3})^{1/2}$. The less compact the star the closest to the Kepler limit is the “atmospheric tail” outside the surface of the star. Although the y-axis starts shorter than the x-axis after ten rotation periods it has “closed the gap” and the two axes have essentially identical angular velocity profiles (this gap is the space between vertical brown dashed and gray dotted lines). This effect is more evident in the G4C010 model but can be clearly seen in the other most compact cases like G4C025.

III.2 Secular fate

Although dynamic stability was straightforward to establish, that has not been the case with secular stability. After evolving for more than twenty rotational periods one can see in Fig. 3 the major characteristics of GW emission. The frequency of the dominant ($l=m=2$) GW mode is twice the rotational frequency, and has amplitude approximately one-tenth of the average value of a merging binary system. The quadrupole-approximation prediction for the GW strain based on the initial configurations is given by Eq. 18, and is shown in the plots as a dashed horizontal green line. This approximate value for the strain is about $50-60\%$ of the maximum amplitude found in the evolution (see also Fig. 4). The GW amplitude in the more compact models (G4C025) experiences a more rapid decrease (almost ten times more) than in the G4C010 case, which has the smallest compactness ($0.1$). Similar behaviour is exhibited in the luminosity and radiated angular momentum plots. In all cases, after an initial period that lasts a little over $500\ M_{\odot}$ ($M_{\odot}=4.927\ {\rm\mu s}$) $\dot{E}$ and $\dot{J}$ intersect the predictions from the quadrupole approximation based on the initial-data (in the plots these are denoted by the horizontal blue and red dashed lines, Eqs 14). However, $\dot{E}$ and $\dot{J}$ undergo exponential decay in close agreement with the corresponding exponential decay in the GW amplitude, i.e.,

[TABLE]

In Fig. 3 we denote the exponential fits for all relevant functions. The evolutionary path of these rotating objects occurs along quasiequilibrium states as seen in the left panel of Fig. 5, which shows $dE=\Omega dJ$. After an initial period of $\sim 500\ M_{\odot}$, this law is satisfied in all cases, apart from a small perturbation at $1050\ M_{\odot}$ for the most supramassive case pwC026. During that period the ratio of the kinetic to gravitational potential energy remains essentially constant and equal to the initial value (see right panel of Fig. 5).

The thermal energy generated by shocks was also measured in these simulations by inspecting the entropy parameter $K:=P/P_{\rm cold}$, where $P_{\rm cold}=k\rho_{0}^{\Gamma}$. With $\epsilon_{\rm th}=\epsilon-\epsilon_{\rm cold}=(K-1)(\Gamma-1)\epsilon_{\rm cold}/(\Gamma_{\rm th}-1)$, then $K>1$ implies shock heated gas ELSB09 . Since we don’t have any mergers in our problem we didn’t expect any shocks, and this was the case for the bulk of the stars ($K\sim 1.0$).

Although we clearly see that triaxially deformed stars evolve in a quasiequilibrium manner towards axisymmetric objects, the key question is whether this evolution is due mainly to GW emission or to a hydrodynamical reconfiguration? Using the exponential fitting functions in Fig. 3 we read off the GW decay timescales. These are $5000M_{\odot}\sim 10^{4}\ M$ for the G4C010 and $900M_{\odot}\sim 500\ M$ for the G4C025 models. The GW driven bar-mode instability occurs for stars rotating with $\beta>\beta_{\rm sec}$ and $\beta_{\rm sec}\approx 0.14$ in the Newtonian incompressible limit. This value decreases in GR as the compactness increases SF98 . The two models discussed here have $\beta=0.15,0.18$ (see Table 2), and are thus greater the Newtonian critical value $\beta_{\rm sec}$. The GW timescale is FS75

[TABLE]

where $\beta_{\rm sec}$ may be approximated by $\beta_{\rm sec}=0.115-0.048M/M_{\rm max}^{\rm sph}$ MSB99 . Here $M_{\rm max}^{\rm sph}$ is the maximum spherical mass for the given EoS. For our cases values are taken from Table 1, which imply that $\tau_{\rm bar}^{\rm GW}\sim 10^{8}\ M$ for the G4C010 and $\tau_{\rm bar}^{\rm GW}\sim 10^{5}M$ for the G4C025 models, respectively. Supramassive models pwC023, pwC024, pwC026 have timescales $\tau_{\rm bar}^{\rm GW}\sim 10^{5}M$ too. We also note that the GW timescales as calculated from the crude quadrupole estimate, $t_{s}\sim T/|\dot{E}|$, and reported in Table 2, are in most cases (except for G4C010) longer than the timescales obtained from Eq. 11. Moreover, our configurations do not evolve toward Dedekind-like ellipsoids as in the case of the bar-mode unstable Newtonian configurations LS1995 ; Miller73 . It is possible that the nonlinear growth of the instability is halted by mode-mode coupling, as our triaxial configuration contains modes beyond $m=2$.

Another possibility for a gravitational wave driven secular instability is the nonaxisymetric r-mode. For the $l=m=2$ mode the timescale is LOM98

[TABLE]

which implies $\tau_{\rm r}^{GW}\sim 10^{7}\ M$ for the G4C010 and $\tau_{\rm bar}^{GW}\sim 10^{5}M$ for the G4C025 respectively. These timescales again are much longer than the timescales found numerically. Also, in this case the wave frequency $f_{\rm GW}=4/3f_{\rm rot}$ therefore this possibility is also ruled out by our data, for which $f_{\rm GW}=2f_{\rm rot}$.

Numerical viscosity, although nonzero, can in principle be responsible. The presence of viscosity can damp a GW radiation reaction-induced bar-mode instability LD77 , although it needs to be properly tuned. However, we evolved with two different resolutions and found no change in the behaviour which might have been expected if numerical viscosity were significant. Also we repeated the calculation with the Whisky code Betal05 ; GR07 ; DAPRG13 ; RRG14 and got very similar results. It may be that even a small numerical viscosity over time is sufficient to damp the mode, given the long timescale ($>>t_{\rm dyn}$) for GW emission. If we modeled numerical viscosity by a turbulent viscosity $\nu\sim\alpha Rc_{s}\sim\alpha(R/M)^{1/2}$, where $c_{s}$ is the sound speed, then a damping time of $10^{4}\ M$ associated with this would only require $\alpha\sim 10^{-3}$ to be effective. Such a small value might go unnoticed by a modest resolution study. On the other hand, if viscosity were to dominate GW dissipation, one still expects that the bar mode will be triggered above $\beta=\beta_{\rm sec}$, since viscosity alone can drive the instability, and the triaxiality would grow Shapiro04 , but this is not observed. Hence we conclude that although our triaxial stars evolve towards axisymmetry, it is not the bar or r-mode secular effects that are mainly responsible for this fate but rather a hydrodynamical reconfiguration of the initial data.

IV Discussion

In this work we investigated the stability properties and gravitational wave signatures of uniformly rotating, triaxial neutron stars in GR. Using the cocal code we have constructed normal as well as supramassive solutions in quasiequilibrium and we evolved them for the first time with the Illinois GRMHD code.

All five solutions that we considered are dynamically stable and evolve secularly towards an axisymmetric configuration. Although we monitored the evolution for more than twenty rotation periods we were unable to probe the final (secular) fate of these stars, which is orders of magnitude longer. We coroborated our findings by using different resolutions, placement of outer boundary conditions, atmospheric treatments, and simulations with a different (Whisky) code.

According to SF98 a perturbed axisymmetric star with $\beta>\beta_{\rm sec}$ will be secularly unstable and develop a bar mode. In our case the initial models already contain a bar perturbation and are rotating beyond the secular bar-mode instability limit, but we found no further growth of a bar mode in the time frame of our simulations, which was shorter than the predicted, theoretical secular timescale. On the contrary we observed the decay of the star’s triaxiality, which is in accordance with previous investigations DLSS04 .

Acknowledgements.

We thank Cecilia Chirenti, Nikolaos Stergioulas and Enping Zhou, for useful discussions. This work was supported by NSF Grants PHY-1300903, PHY-1602536, PHY-1607449, NASA Grants NNX13AH44G, NNX16AR67G (Fermi), and JSPS Grants-in-Aid for Scientific Research (C) No. 26400274, and 15K05085. VP gratefully acknowledges support from the Simons foundation. This work made use of the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number TG-MCA99S008. This research is part of the Blue Waters sustained-petascale computing project, which is supported by the National Science Foundation (awards OCI-0725070 and ACI-1238993) and the State of Illinois. Blue Waters is a joint effort of the University of Illinois at Urbana-Champaign and its National Center for Supercomputing Applications. Some numerical computations were carried out on the XC30 system at the Center for Computational Astrophysics (CfCA) of the National Astronomical Observatory of Japan.

Appendix A Quadrupole formulae in helical symmetry

For an estimate of the GWs one can compute the time derivatives of the quadrupole mass moments. Typically the quadrupole formula reads

[TABLE]

where \mathrel{\mathchoice{{\hbox to0.0pt{\hbox to5.18054pt{\hss--\hss}\hss}\hbox{I}}}{{\hbox to0.0pt{\hbox to5.18054pt{\hss--\hss}\hss}\hbox{I}}}{\scriptsize{\hbox to0.0pt{\hbox to3.62637pt{\hss--\hss}\hss}\hbox{I}}}{\tiny{\hbox to0.0pt{\hbox to2.59027pt{\hss--\hss}\hss}\hbox{I}}}}_{ij}:=I_{ij}-\frac{1}{3}f_{ij}I_{kk}, and \mathrel{\mathchoice{{\hbox to0.0pt{\hbox to5.18054pt{\hss--\hss}\hss}\hbox{I}}}{{\hbox to0.0pt{\hbox to5.18054pt{\hss--\hss}\hss}\hbox{I}}}{\scriptsize{\hbox to0.0pt{\hbox to3.62637pt{\hss--\hss}\hss}\hbox{I}}}{\tiny{\hbox to0.0pt{\hbox to2.59027pt{\hss--\hss}\hss}\hbox{I}}}}_{ij}^{TT} is the transverse traceless reduced quadrupole moment MTW . The second time derivatives are computed at a retarded time. The corresponding gravitational wave luminosity and the angular momentum carried away per unit time are

[TABLE]

where $\langle\cdot\rangle$ denote an average over several wavelengths. In full dynamical spacetimes there is no unique definition of the quadrupole moment but typically one uses Eq. (6) as a generalized integral over the hypersurface $\Sigma$, SS03 which can be thought as an Euclidean integral over a weighted density $\rho_{\ast}=\rho_{0}u^{t}\sqrt{-g}$. Its time derivative

[TABLE]

can be obtained by using the conservation of rest mass $\partial_{t}\rho_{\ast}+\partial_{i}(\rho_{\ast}v^{i})=0$, and integration by parts FE90 .

Another way to obtain the same result is to employ the transport theorem that says that for any density $\rho_{\ast}$ that satisfies the above continuity equation and any function $Q(t,x^{i})$, we have

[TABLE]

where $\frac{DQ}{Dt}=\partial_{t}Q+v^{i}\partial_{i}Q$ is the Lagrangian derivative of $Q$. For a fluid velocity $v^{i}=\Omega\phi^{i}$, we have $DQ/Dt=\mathcal{L}_{\bf k}Q$, and thus we can write a fully 4-dim version of the classical theorem as

[TABLE]

A straightforward proof of Eq. 17 can be obtained if we consider $f(t)=\int_{\Sigma_{t}}Q\rho_{0}u^{\alpha}dS_{\alpha}$. Let $\Sigma=\Sigma_{0}$ and $\Sigma_{t}=\psi_{t}(\Sigma)$, where $t^{\alpha}$ is the generator of the diffeomorphism family $\psi_{t}$. Then

[TABLE]

To obtain the last line we converted the second integral in the previous line over a divergence, to a surface integral that vanishes, and also used the continuity equation in the form $\mathcal{L}_{\bf k}(\rho_{0}u^{t}\sqrt{-g})=0$.

For the computation of Eq. (14) we need to compute the third material derivatives of $x^{i}x^{j}$. We denote by $\phi^{i}=(\phi^{A},0)$ where capital letters take values in $\{1,2\}$. Then $\phi^{A}=-\epsilon^{AB}x_{B}$ and the nonzero components are

[TABLE]

Setting $\varpi^{i}=(x^{A},0)$ we have

[TABLE]

Using the derivatives of the multiple moments above one can compute the luminosity or the angular momentum radiated from Eq. 14. For the GW strain, assuming rotation around the z-axis, we have

[TABLE]

For the case of an exact triaxial ellipsoid the two elliptical polarization modes for head on observation along the z-axis, we set

[TABLE]

where $I_{k}$ are the principal moments of inertia. Then the emitted power and angular momentum will be,

[TABLE]

A parameter which is often mentioned is called ellipticity of the source is defined as $\varepsilon:=|I_{1}-I_{2}|/I_{3}$. Although there is no rigorous counterpart in GR we can generalize as

[TABLE]

This is the quantity that is reported in Table 2.