# Uniformly rotating, axisymmetric and triaxial quark stars in general relativity

###### Abstract

Quasi-equilibrium models of uniformly rotating axisymmetric and triaxial quark stars are computed in general relativistic gravity scenario. The Isenberg-Wilson-Mathews (IWM) formulation is employed and the Compact Object CALculator (cocal) code is extended to treat rotating stars with finite surface density and new equations of state (EOSs). Besides the MIT bag model for quark matter which is composed of de-confined quarks, we examine a new EOS proposed by Lai and Xu that is based on quark clustering and results in a stiff EOS that can support masses up to in the case we considered. We perform convergence tests for our new code to evaluate the effect of finite surface density in the accuracy of our solutions and construct sequences of solutions for both small and high compactness. The onset of secular instability due to viscous dissipation is identified and possible implications are discussed. An estimate of the gravitational wave amplitude and luminosity based on quadrupole formulas is presented and comparison with neutron stars is discussed.

## I Introduction

The recent gravitational-wave (GW) event GW170817 together with accompanying electromagnetic emission observations The LIGO Scientific Collaboration and The Virgo Collaboration (2017); The LIGO Scientific Collaboration et al. (2017a) from a binary neutron star (BNS) merger has opened a brand new multi-messenger observation era for us to explore the Universe. Apart from enriching our knowledge on origins of short gamma-ray bursts LIGO Scientific Collaboration et al. (2017) and nucleosynthesis associated with BNS mergers The LIGO Scientific Collaboration et al. (2017b); Baiotti and Rezzolla (2017), it also provides an effective way for us to constrain the equation of state (EOS) of neutron stars (NS). In addition to systems such as binary black-hole mergers and BNS mergers, rapidly rotating compact stars have also been considered as important candidates of GW sources Andersson et al. (2011), which could be detected by ground-based GW observatories Abramovici et al. (1992); Punturo et al. (2010); Accadia et al. (2011); Kuroda and LCGT Collaboration (2010); Aso et al. (2013) and help us understand the nature of strong interaction of dense matter.

It has been long since the equilibrium models of self-gravitating, uniformly rotating, incompressible fluid stars were systematically studied in a Newtonian gravity scheme Chandrasekhar (1969). Depending on the rotational kinetic energy, the configuration could be axisymmetric Maclaurin ellipsoids as well as nonaxisymmetric ellipsoids, such as Jacobian (triaxial) ellipsoid. For compact stars that we are interested in for GW astronomy, however, general relativity is required to replace Newtonian gravity. The field of relativistic rotating stars has been studied for many years Meinel et al. (2008); Friedman and Stergioulas (2013).

A rotating NS will spontaneously break its axial symmetry if the rotational kinetic energy to gravitational binding energy ratio, exceeds a critical value. This instability can either be of secular type Shapiro and Zane (1998); Bonazzola et al. (1996); Gondek-Rosińska and Gourgoulhon (2002); Bonazzola et al. (1998) or dynamical Houser and Centrella (1996); Pickett et al. (1996); Brown (2000); New et al. (2000); Liu (2002); Watts et al. (2005); Baiotti et al. (2007); Manca et al. (2007); Corvino et al. (2010), depending on the process driving the instability and with only small modifications if a magnetisation is present Camarda et al. (2009); Franci et al. (2013); Muhlberger et al. (2014) (see Andersson (2003) for a review). A high ratio can also be reached for a newly born rotating compact star during a core collapse supernova or for a NS which is spun up by accretion Lai and Shapiro (1995); Bildsten (1998); Woosley and Janka (2005); Watts et al. (2008); Piro and Thrane (2012).

Quasi-equilibrium figures of triaxially rotating NSs have also been created and studied in full general relativity Huang et al. (2008); Uryū et al. (2016). In this case, the bifurcation from an axisymmetric to triaxial configuration happens very close to the mass shedding limit, and, for soft NS EOSs or for NSs with large compactness, the triaxial sequence could totally vanish James (1964); Bonazzola et al. (1996, 1998); Hachisu and Eriguchi (1982); Lai et al. (1993). As a result, it is presently unclear whether triaxial configurations of NSs can actually be realized in practice.

On the other hand, it is worth noting that the EOS of compact stars is still a matter of lively debate since astronomical observations are not sufficient to rule out many of the nuclear-physics EOSs that are compatible with the observations. As a result, besides the popular idea of NSs, other models for compact stars are possible and have been considered in the past. A particularly well developed literature is the one concerned with strange quark stars (QSs), since it was long conjectured that strange quark matter composed of de-confined up, down and strange quarks could be absolutely stable Bodmer (1971); Witten (1984). There is also possible observational evidence indicating the existence of QSs (for a recent example, see Dai et al. (2016)). Additionally, the small tidal deformability of QSs is favoured by the observation of GW170817 Lai et al. (2017) and possible models with QS merger or QS formations are also suggested to explain the electromagnetic counterparts for a short gamma-ray burst (c.f. Li et al. (2016); Lai et al. (2017)).

Following this possibility, a large effort has been developed to calculate equilibrium configurations of QSs, starting from the first attempts Itoh (1970); Alcock et al. (1986); Haensel et al. (1986). At present, both uniformly rotating Gondek-Rosińska et al. (2000); Gourgoulhon et al. (1999); Stergioulas et al. (1999) and differentially rotating QSs Szkudlarek et al. (2012) have been studied in full general relativity. Unlike NSs, which are bound by self-gravity, QSs are self-bound by strong interaction. Consequently, rotating QSs can reach a much larger ratio compared with NSs and the triaxial instability could play a more important role Gondek-Rosinska et al. (2000); Gondek-Rosińska et al. (2000, 2001). The triaxial bar mode (Jacobi-like) instability for MIT bag-model EOS has been investigated in a general relativistic framework Gondek-Rosińska et al. (2003).

We here use the Compact Object CALculator code, cocal, to build general-relativistic triaxial QS solution sequences using different EOS models. cocal is a code to calculate general-relativistic equilibrium and quasi-equilibrium solutions for binary compact stars (black hole and NS) as well as rotating (uniformly or differentially) NSs Uryū and Tsokaros (2012); Uryū et al. (2012); Tsokaros and Uryū (2012); Tsokaros et al. (2015); Huang et al. (2008); Uryū et al. (2016). The part of the cocal code handling the calculation of the EOS was originally designed for piecewise polytropic EOSs. We have here extended the code to include polynomial type of EOSs, as those that can be used to describe QSs. In doing so the trivial relationship between the thermodynamic quantities for a piecewise polytrope (e.g., see Eqs. (64)–(68) in Tsokaros et al. (2015)), is lost and now one has to apply root finding methods. Another issue is related to the surface fitted coordinates that are used in cocal to track the surface of the star. For NSs, the surface was identified as the place where the rest-mass density goes to zero or where the specific enthalpy becomes one. This is no longer generally true for a self-bound QS and a different approach needs to be developed. The nonlinear algebraic system that determines the angular velocity, the constant from the Euler equation, and the renormalization constant of the spherical grid has to be modified in order to accommodate the arbitrary surface enthalpy.

We here compute solutions for both axisymmetric and triaxial rotating QSs with the new code, as well as sequences with various QS EOS and different compactnesses. We checked our new implementation for those cases were previous studies have been possible Huang et al. (2008), and we confirm the accuracy of our new code. We discuss the astrophysical implications of the quantities of rotating QSs at the bifurcation point. For instance, the spin frequency at the bifurcation point could be a more realistic spin up limit for compact stars rather than the mass shedding limit, which relates to the fastest spinning pulsar we might be able to observe. The GW strain and luminosity estimates for our models are given, while full numerical simulations are left for the future (see Tsokaros et al. (2017) for recent simulations involving triaxial NSs).

The structure of this paper is organized as follows: In Sec. II we discuss the formulation we used and the field equations (Sec. II.1), the hydrodynamics (Sec. II.2) and the EOS part (Sec. II.3). In order to test the behavior of the modified code, we have performed convergence tests with five resolutions and compared with rotating NS solutions built by the original cocal code. These tests can be found in Sec. III. Triaxially deformed rotating QS sequences for different compactnesses, are presented in Sec. IV, while the implications for the astrophysical observations of this work are presented in Sec. VI. Hereafter we use units with unless otherwise stated; a conversion table to the standard cgs units can be found, for instance, in Rezzolla and Zanotti (2013).

## Ii Formulation and numerical method

### ii.1 Field equations

In order to solve the field equations numerically, the Isenberg-Wilson-Mathews (IWM) formulation Isenberg and Nester (1980); Isenberg (2008); Wilson and Mathews (1989) is employed. In a coordinate chart , the decomposition of the spacetime metric gives

(1) |

where are the lapse and shift vector (the kinematical quantities), while is the IWM approximation for the three-metric.

The extrinsic curvature of the foliation is defined by

(2) |

and a maximal slicing condition is assumed.

Decomposing the Einstein equations with respect to the normal of foliation, we get the following 5 equations in terms of the five metric coefficients on the initial slice :

(3) | |||

(4) | |||

(5) |

where the first and second equations are the Hamiltonian and momentum constraints, respectively. Here is the projection tensor onto the spatial slices. These equations can be written in the form of elliptic equations with the non-linear source terms, respectively,

(6) | |||

(7) | |||

(8) |

where , and the source terms of matter are defined by , , and .

The above set of equations must be supplied with boundary conditions at infinity. Since we are working in the inertial frame and we impose asymptotic flatness, we must have

(9) |

### ii.2 Hydrostatic equilibrium

The hydrostatic equation for a perfect fluid in quasi-equilibrium can be derived from the relativistic Euler equation Rezzolla and Zanotti (2013)

(10) |

where is the 4-velocity of the fluid, , and is the specific enthalpy defined by ( is the rest-mass density and the total energy density).

When the symmetry along a helical Killing vector is imposed for the fluid variables, which is approximately true also in the case for a rotating nonaxisymmetric star in quasi-equilibrium, the integral of the Euler equation becomes

(11) |

where is a constant. From the normalization of the four velocity , one obtains

(12) |

where . The fluid sources of Eqs. (6)–8, i.e., , and , are defined in terms of the energy momentum tensor in the previous section. In terms of the fluid and field variables they can be written as Rezzolla and Zanotti (2013)

(13) | |||

(14) | |||

(15) |

in which is the relativistic analogue of the Emden function. Here is related to through Eq. (11). Therefore in order to close the system, an additional relationship is needed between the specific enthalpy, the pressure and the rest-mass density of the fluid, i.e., an EOS. Once such a relation is available, to solve the field equations [Eqs. (6)–8] and the hydrostatic equation [Eq. (11)] one has to find the two constants that appear in all of them. This procedure is described in detail for example in Ref. Tsokaros et al. (2015).

### ii.3 Equation of State

In this work, we have considered two types of EOS for QSs. One of them is the MIT bag-model EOS Chodos et al. (1974), since it is the most widely used EOS for QSs. In the case when strange quark mass is neglected, the pressure is related to total energy density according to

(16) |

where are two constants, the second being the total energy density at the surface. Related to is the so called bag constant, . In this work, and following Limousin et al. (2005), the simplest MIT bag-model EOS has been employed, where and .

Besides the MIT bag-model EOS, we have also considered another QS EOS
suggested by Lai and Xu Lai and Xu (2009), which we will refer to as the LX
EOS hereafter. Unlike the conventional QS models (e.g., the MIT bag-model
EOS) which are composed of de-confined quarks, Lai and Xu Lai and Xu (2009)
suggested that quark clustering is possible at the density of a cold
compact star since the coupling of strong interaction is still decent at
such energy scale. Due to the non-perturbative effect of strong
interaction at low energy scales and the many-body problem, it is very
difficult to derive the EOS of such a quark-cluster star^{1}^{1}1Such a
quark-cluster star has also been named a strangeon star in
Ref. Lai and Xu (2017). from first principles.

Lai and Xu attempted to approach the EOS of such quark-cluster star with phenomenological models, i.e., to compare the intercluster potential with the interaction between inert molecules (a similar approach has also been discussed in Guo et al. (2014)). They also take the lattice effects into account as the potential could be deep enough to trap the quark clusters. Combining the inter-cluster potential and the lattice thermodynamics, they have derived an EOS in the following form:

(17) |

where is the reduced Planck constant. The parameters in this expression, and , are the depth of the potential and characteristic range of the interaction, respectively. The EOS is also dependent on the number of quarks in each cluster () since it relates the energy density () and rest-mass density () to the number-density of quark clusters ( in Eq. (17)). Similarly to the MIT bag-model EOS case, we use the rest-mass density parameter, which is

(18) |

where is the atomic mass unit. While several different choices of parameters are considered in Ref. Lai and Xu (2009), in our work we restrict our attention to and for . We also note that although it is not as obvious as for the MIT bag-model EOS, the LX EOS also has a nonzero surface density since expression Eq. (17) has a unique zero root when the number density is positive.

Being a stiff EOS, the LX EOS is favored by the discovery of massive pulsars Lai and Xu (2011); Demorest et al. (2010); Antoniadis et al. (2013). The rest-mass density and mass-radius relationships for spherical models can be seen in Fig. 1, and the characteristics of the maximum mass models are reported in Table 2. The LX EOS has also been discussed in relation with the possibility of understanding some puzzling observations related to compact stars, such as the energy release during pulsar glitches Zhou et al. (2014), the peculiar X-ray flares Xu and Liang (2009) and the optical/UV excess of X-ray-dim isolated NSs Wang et al. (2017). Particularly, a solid QS model has been suggested in order to understand those observations Xu (2003). However, as pointed out in Zhou et al. (2014), the critical strain of such a star is very small. A starquake will be induced when the relative difference in ellipticity is , for most, between the actual configuration of the star and the configuration as if the star is a perfect fluid. This is consistent with the pulsar glitch observations on Vela. Therefore, we find it a good approximation to calculate the quasi-equilibrium configuration of such a star with perfect fluid assumption.

Understanding the models and properties of the QS EOSs that we want to consider, we can modify the EOS part of the simulation code, which was originally designed for NS models, accordingly.

As mentioned above, in the case of NSs, a piecewise-polytropic EOS is usually assumed to describe the EOS Read et al. (2009); Rezzolla and Zanotti (2013). In each piece, the pressure and rest-mass density are related as

(19) |

For QSs, due to the nonzero surface density and nonzero energy density integration constant, we will assume that the EOS is generally a polynomial

(20) |

Given the relationship between and , one can apply the first law of thermodynamics to obtain other quantities such as the energy density and the specific enthalpy. In the zero-temperature case, the first law of thermodynamics can be expressed as Rezzolla and Zanotti (2013)

(21) |

which can be integrated to obtain the total energy density. The integral constant is usually chosen to be 1, since when there is no internal energy, the energy density and the rest-mass density coincide (apart from the square of the speed of light). However, for QS EOSs (like in Guo et al. (2014) and the MIT bag-model EOS), the integral constant is different from unity and needs to be properly taken into account.

MIT bag-model EOS | ||

LX EOS | ||

In this case, the energy density and specific enthalpy are related to the rest-mass density by

(22) |

(23) |

Here is the integral constant we mentioned above. It is usually taken to be zero for NS models. Here it is introduced again in order to accommodate stars that have different surface limit for the thermodynamic variables.

The nonzero surface density of QSs requires a different boundary condition in our simulation. For typical NSs when we adjust the position of the surface, we are actually locating the points where the specific enthalpy is . For QSs, the surface identification will be at values of the specific enthalpy different from unity and consistent with Eq. (23). What we typically use as input parameter is the surface rest-mass density , from which we then calculate using Eq. (23).

As an example, for MIT bag-model EOS the first law at zero temperature implies

(24) | |||||

(25) | |||||

(26) |

where is a constant of integration. The above EOS is of the form (20) with

(27) |

and . Having all thermodynamical variables in terms of the rest-mass density is convenient from the computational point of view since this is one of the fundamental variables used in the cocal code, therefore the modifications with respect to the EOS will be minimal.

Given a fixed choice of and for the MIT bag-model EOS, one can obtain a unique solution of the field equations under hydrostatic equilibrium (see Secs. II.1 and II.2). In other words, the relationship between the gravitational mass versus the central energy density, the mass-radius relationship and the spacetime metric will not depend on the coefficient , as it is eliminated out from Eq. (16) (similar argument can be found in Li et al. (2017); Bhattacharyya et al. (2016)). At the same time, it will indeed affect the rest mass hence the binding energy of the QS since it relates the rest-mass density and the number density of the components. Hence a reasonable choice for will still be helpful although it will not affect anything that we are interested in for this work. Here we chose such that the EOS corresponds to the model as in Alford et al. (2005) and we assume a rest mass of for each baryon number (, where is the number density of quarks). Any other choices for are in principle possible and they will not affect our solution except for the rest mass of the star. Actual values of these constants can be found in Table 1 in cgs units.

In Eqs. (24) and (25) we have employed the relationship which is quite similar to the explicit form of MIT bag model. By factoring out the rest mass, those two equations can be rewritten as a function of number density instead of rest-mass density. However, one clarification we want to discuss at this point is that, the choice of using Eq. (24) and (25) is not essential. Moreover, Eqs. (24) and (25) are related by the first law of thermodynamics [see Eq. (21)], which is not essential as well. We can describe the MIT EOS Eq. (16), in a parametric form with an arbitrary parameter , as long as the functions and satisfy Eq. (16). In doing so we choose to satisfy Eq. (21) and therefore arrive at Eq. (24), (25). As one can see from Eqs. (13)–(15), the fluid terms that appear are and . Thus, for the MIT bag-model EOS, the only thermodynamic variable that appears in the field equations is . Every model thus calculated will be uniquely defined by a deformation parameter and the central total energy density. The scaling constant which is analogous to scaling as for NSs Rezzolla and Zanotti (2013), is here .

EOS | |||||
---|---|---|---|---|---|

MIT | 0.2940 | 2.217 | 0.2706 | ||

LX | 2.326 | 3.325 | 0.3956 |

Before concluding this discussion on the EOSs we should mention that although a polynomial-type EOS is implemented for the calculation of QSs with a finite surface density, the developments in the new version of cocal allow us to calculate any compact star with an EOS that can be described by a polynomial function, including NSs and hybrid stars. For instance, some phenomenological approaches suggested recently in Ref. Baym et al. (2017) also result in a polynomial-type EOS, which can be computed straightforwardly with the new code.

## Iii Code tests

The working properties of the cocal code for single rotating stars is
presented in detail in previous works, e.g., Huang et al. (2008); Uryū et al. (2016)^{2}^{2}2See
Uryū and Tsokaros (2012); Uryū et al. (2012); Tsokaros and Uryū (2012); Tsokaros et al. (2015) for the general
binary case., so that here we will only mention the most important
quantities that are used in our simulations. The method has its origins
in the works of Ostriker and Marck Ostriker and Mark (1968), who used it to
compute Newtonian stars, and the works of Komatsu, Eriguchi, and Hachisu
Komatsu et al. (1989), who devised a stable numerical algorithm and obtained
first axisymmetric general relativistic rotating stars. From this latest
work the method is commonly referred as the KEH method
and consists of an integral representation of the Poisson equation
commonly referred as the representation formula. Since we have only one
computational domain with trivial boundary conditions at infinity, the
Green’s function is and is expanded as a series of the
associated Legendre polynomials and trigonometric functions. The maximum
number of the terms included in this expansion is given as in Table
3.

This approach is used to compute the gravitational fields
while hydrostatic equilibrium is achieved
through Eq. (11). At every step in the iteration to reach
the solution at the desired accuracy, three constants need to be
computed. The first one is the angular velocity of the star, ,
while the second one is the constant from the Euler integral,
, in Eq. (11), and, finally, the third
constant is , a normalization factor for the whole domain where the
equations are solved.^{3}^{3}3We recall that cocal uses normalized
variables and the quantities listed in Table
4 refer to those and should be denoted
by a hat. For simplicity, however, we have omitted these hats in the
Table. At every step during the iteration the nonlinear equation with
respect to these three constants is solved typically by evaluating
Eq. (11) at three points in the star. For axisymmetric
configurations, we use the center of the star and two points on the
surface, one on the positive -axis and one at the North pole of the
stellar model. An axisymmetric equilibrium is achieved by setting the
ratio between the polar axis over the equatorial radius along the
-axis. For triaxial configurations, on the other hand, the three
points are the center of the star together with two points again on the
surface, one at the positive -axis, and one on the positive
-axis. Each triaxial solution has a fixed ratio of the radius on the
-axis over the radius of the -axis.

From a numerical point of view, cocal is a finite-difference code
that uses spherical coordinates ^{4}^{4}4Note that the field equations for the
shift vector Eq. (7) are expressed in Cartesian
coordinates. and the basic parameters are summarized in Table
3. In the angular directions the
discretization is uniform, i.e., , and
. In the radial direction the grid is uniform
until point with , and in the interval
the radial grid in non-uniform and follows a geometric series
law Huang et al. (2008). While field variables are evaluated at the
gridpoints, source terms under the integrals are evaluated at midpoints
between two successive gridpoints since the corresponding integrals use
the midpoint rule. For integrations in and we use a
second-order midpoint rule. For integrations in we use a
fourth-order midpoint rule. This was proven necessary to keep second
order convergence at the region of maximum field strength
Uryū et al. (2012). Derivatives at midpoints are calculated using
second-order rule for the angular variables , and
third-order order rule for the radial variable (again for keeping
second-order convergence at same regions Uryū and Tsokaros (2012)). Derivatives
evaluated at gridpoints always use a fourth-order formula in all
variables.

: | Radial coordinate where the grid starts. | |

: | Radial coordinate where the grid ends. | |

: | Radial coordinate between and where the | |

grid changes from equidistant to non-equidistant. | ||

: | Total number of intervals between and . | |

: | Number of intervals in . | |

: | Number of intervals in . | |

: | Total number of intervals for . | |

: | Total number of intervals for . | |

: | Number of multipole in the Legendre expansion. |

Type | |||||||||
---|---|---|---|---|---|---|---|---|---|

H2.0 | 0 | 1.25 | 192 | 80 | 64 | 48 | 48 | 12 | |

H2.5 | 0 | 1.25 | 288 | 120 | 96 | 72 | 72 | 12 | |

H3.0 | 0 | 1.25 | 384 | 160 | 128 | 96 | 96 | 12 | |

H3.5 | 0 | 1.25 | 576 | 240 | 192 | 144 | 144 | 12 | |

H4.0 | 0 | 1.25 | 768 | 320 | 256 | 192 | 192 | 12 |

It is worth noting that, since for QSs the relationship between the specific enthalpy and the rest-mass density follows a general polynomial function [cf., Eq. (23)], a root-finding method needs to be employed when calculating the thermodynamical quantities from the enthalpy. The regular polynomial expression of the specific enthalpy with respect to rest-mass density allows us to use a Newton-Raphson method as the derivative can also be expressed easily. In view of this, the computational costs with a QS EOS are not significantly larger than those with NS EOS. However, in order to guarantee a solution of rest-mass density when the specific enthalpy is given, a bi-section root-finding method needs to be employed if the Newton-Raphson method does not converge sufficiently rapidly. In this case, the initial range of the bi-section method is set to be the specific enthalpy corresponding to the rest-mass densities at the stellar center and at the surface, respectively.

In previous works Huang et al. (2008); Uryū et al. (2012, 2016); Tsokaros et al. (2016), the cocal code has been extensively tested, both with respect to its convergence properties, as well as with respect in actual evolutions with other well established codes Tsokaros et al. (2016).

In what follows we report the convergence tests we have performed in order to investigate the properties of the code under these new conditions. Before doing that, we note that special care is needed when using a root-finding method to calculate thermodynamical quantities for a given specific enthalpy in the case of rotating QSs. In particular, it is crucial to consider what is the accuracy set during the root-finding step. Of course, it is possible to require that the accuracy in the root-finding step is much higher than the other convergence criteria in the code to guarantee an accurate result. In this way, however, the computational costs will increase considerably, since the thermodynamical quantities need to be calculated at every gridpoint and at every iteration. We found that an accuracy of for the thermodynamic variables solutions neither compromises the accuracy of the solutions, nor slows down the code significantly.

### iii.1 Comparison with rotating NSs

Although the newly developed code presented here is intended for QS EOSs, it can be also used to produce rotating NSs if one restricts the EOS to a single polytrope. This can be accomplished by setting the polynomial terms to be only one, the surface rest-mass density and the energy integral constant . In this case, Eq. (20) becomes

(28) |

and the relationship between the energy density, the specific enthalpy and the rest-mass density [Eqs. (22) and 23] will be exactly the same as that for a polytropic NS.

We choose a stiff EOS with and produce axisymmetric solution sequences for small and high compactness , and is the corresponding Arnowitt-Deser-Misner (ADM) mass for a nonrotating model, both with the original rotating-NS and the modified rotating-QS solver. The grid-structure parameters used are , , , , , , , and (see Table 3). Overall, we have found that the relative difference in all physical quantities is of the order of , which is what is expected since the criteria for convergence in cocal are that the relative difference in metric and fluid variables between two successive iterations is less than .

### iii.2 Convergence test for rotating QSs

For the convergence analysis in this work we use the five resolutions shown in Table 4. The outer boundary of the domain is placed at , while the surface of the star is always inside the sphere . The radius along the -axis is exactly in the normalized variables. There are exactly intervals along the radii in the , , and directions. The number of Legendre terms used in the expansions is kept constant () in all resolutions since convergence with respect of those has been already investigated in Tsokaros and Uryū (2007). When going from the low-resolution setup H2.0 to the high-resolution one H4.0, the spacings decrease as .

As a result, if we denote as , a quantity evaluated at two different resolutions, then

(29) |

where is a constant and is the grid separation at resolution . Choosing the combinations , , and so that we have and normalizing by we plot in Fig. 2 the relative error with respect to the grid spacing for , , , , and the eccentricity , both for the LX EOS (left panel) and for the MIT bag-model EOS (right plot). The deformation is kept at for both EOSs, while the central densities are , and for LX and MIT bag-model EOS respectively. The dashed black line reports a reference first-order convergence, while the solid black line refers to second-order convergence.

Note that quantities like the ADM mass, the angular velocity and the eccentricity converge to second order, while quantities like the angular momentum, the ratio , and the rest-mass converge to an order that is closer to first. Furthermore, Fig. 2 shows that some quantities (e.g., the ADM mass of the LX EOS) shows a convergence order that is larger than second, but this is an artefact of the specific deformation. In general, we found second-order convergence in and at least first order for .

Note also that the two panels in Fig. 2 are very similar, even though the EOSs are quite different, with the MIT bag-model EOS being relatively soft (i.e., ), while the LX EOS is comparatively stiff (i.e., ). Hence, the overall larger error that is reported in Fig. 2 when compared to the corresponding Fig. 1 in Ref. Huang et al. (2008), is mostly due to the finite rest-mass density at the stellar surface. In the original rotating-NS code, in fact, the surface was determined through a first-order interpolation scheme. This approach, however, is not sufficiently accurate for rotating QSs and would not lead to the desired convergence order unless the surface finder scheme was upgraded to second order.

## Iv Triaxial solutions

The onset of a secular instability to triaxial solutions for the MIT bag-model EOS stars has been studied previously via a similar method in Ref. Gondek-Rosińska et al. (2003). Surface-fitted coordinates have been used to accurately describe the discontinuous density at the surface of the star, and a set of equations similar to the one of the conformal flat approximation used here was solved. In order to find the secular bar-mode instability point, the authors of Ref. Gondek-Rosińska et al. (2003) performed a perturbation on the lapse function of an axisymmetric solution and build a series of triaxial quasi-equilibrium configurations to see whether this perturbation is damped or grows.

Here, we build quasi-equilibrium sequences with constant rest mass (axisymmetric and triaxial) for both the MIT bag-model EOS and the LX EOS. We begin with the axisymmetric sequence in which we calculate a series of solutions with varying parameters, i.e., the parameters that determine the compactness (e.g., the central rest-mass density ) and the rotation (). In doing so, we impose axisymmetry as a separate condition and manage to reach eccentricities as high as for and compactness . In order to access the triaxial branch of solutions, we recompute the above sequence of solutions but this time without imposing axisymmetry. As the rotation rate increases ( decreases) the triaxial deformation () is spontaneously triggered, since at large rotation rate the triaxial configuration possesses lower total energy and is therefore favoured over the axisymmetric solution. This approach is different from the approach followed in Ref. Gondek-Rosińska and Gourgoulhon (2002), where the triaxial perturbation was triggered after a suitable modification of a metric potential.

We keep decreasing to reach the mass shedding limit with the triaxial configuration. We can then move along the triaxial solution sequence by increasing which now acts as the new rotating parameter. The sequence is then terminated close to the axisymmetric sequence. The bifurcation point can be found by extrapolating this triaxial sequence towards the axisymmetric solutions. The largest triaxial deformation calculated in this work, for both the MIT bag-model EOS and the LX EOS, is for the case (a three-dimensional of the surface for this solution is shown in Fig. 3), for , and for . Similar with NSs Huang et al. (2007), the endpoint of the triaxial sequence happens in lower eccentricities as the compactness increases.

In Figs. 4 and 5, the relation between the ratio versus the eccentricity of the star has been plotted for three different compactnesses () for both the MIT bag-model EOS and the LX EOS. Unlike in a Newtonian incompressible star, for which the bifurcation to triaxial deformation happens at for any compactness, in general relativity the bifurcation point depends on the compactness. According to Gondek-Rosińska and Gourgoulhon (2002),

(30) |

This relation holds true not only for NSs but also for QSs with the MIT bag-model EOS (see Fig. 1 of Gondek-Rosińska et al. (2003)). The largest for the onset of secular instability is found to be for rotating QSs in the configurations that we considered for both the MIT bag-model EOS and the LX EOS, and it will be even larger for higher compactnesses. Both the LX EOS and the MIT bag-model EOS in our calculations follow this relationship within a maximum error of . This implies that the secular instability to a ”Jacobi type” ellipsoidal figure in general relativity is not particularly affected by the stiffness of the EOS for quark matter.

EOS | ||||||||||

MIT | 0.1 | 7.021 (8.077) | 0.5647 (0.5693) | 0.02808 | 0.6515 | 0.4580 | 0.1520 | 16.31 | 0.1343 | |

MIT | 0.15 | 7.962 (9.922) | 0.5565 (0.5640) | 0.02985 | 1.169 | 1.293 | 0.1609 | 43.30 | 0.2231 | |

MIT | 0.2 | 8.415 (11.43) | 0.5478 (0.5590) | 0.03199 | 1.731 | 2.649 | 0.1706 | 82.83 | 0.3308 | |

LX | 0.1 | 5.698 (6.557) | 0.5644 (0.5689) | 0.03451 | 0.5312 | 0.3039 | 0.1518 | 8.805 | 0.1343 | |

LX | 0.15 | 6.515 (8.130) | 0.5566 (0.5639) | 0.03630 | 0.9686 | 0.8850 | 0.1607 | 24.38 | 0.2251 | |

LX | 0.2 | 6.972 (9.528) | 0.5469 (0.5574) | 0.03838 | 1.481 | 1.932 | 0.1715 | 50.34 | 0.3401 | |

0.1 | 6.624 (7.634) | 0.5634(0.5693) | 0.03180 | 0.5841 | 0.3708 | 0.1507 | 11.66 | 0.1328 | ||

0.2 | 7.312 (9.979) | 0.5394(0.5535) | 0.03926 | 1.435 | 1.835 | 0.1688 | 46.72 | 0.3311 | ||

0.1 | 10.30 (11.83) | 0.5461(0.5536) | 0.02197 | 0.8416 | 0.7644 | 0.1493 | 34.80 | 0.1281 |

It is worth noting that when compared with the rotating NSs calculated in Ref. Uryū et al. (2016), rotating QSs have longer triaxial sequences. In another word, the triaxial sequence of rotating QSs terminates at larger eccentricity as well as larger triaxial deformation (in another word, smaller ratio). A rotating NS with and compactness bifurcates from axisymmetry at and can rotate as fast as to reach eccentricities (see Fig. 6 in Uryū et al. (2016)). For the QS models considered here and both EOSs, we have a bifurcation point at and the mass shedding limit at . For more compact NSs with , the bifurcation point happens at and the mass shedding limit at . The corresponding compactness QS models bifurcate at and rotate as fast as .

A few remarks are useful to make at this point. First, we note that these values of eccentricity are strictly valid under the assumption of the conformal flatness approximation, which is however accurate for smaller compactnesses. These estimates are less accurate when the compactness increases and are slightly different when adopting more accurate formulations, such as the waveless approximation (see Fig. 6 in Uryū et al. (2016)). Second, another difference between triaxial NSs and QSs is that for triaxial NSs, the ratio is essentially constant along the triaxial sequence, especially for higher compactnesses (for lower compactnesses there is an increase towards the mass shedding limit, but this is very slight). For rotating triaxial QSs, on the other hand, although this qualitative behaviour is still true, a greater curvature towards higher ratios can be seen. For example, for the models mentioned above, the difference between the critical value of and the one at mass-shedding limit is for NSs while it is for QSs. Third, in Ref. Gondek-Rosińska et al. (2003) it has been shown that at the bifurcation point the relation between the scaled angular frequency, where ), and the scaled gravitational mass , depends only very weakly on the bag constant. If we consider models with compactness model (see top line in Table 5), and , while the scaled bifurcation frequency for such a scaled mass model is roughly as deduced from Fig. 7 in Gondek-Rosińska et al. (2003). Similarly, for the models, the renormalized ADM mass and frequency are and , while the corresponding range is in Ref. Gondek-Rosińska et al. (2003); finally, for the case, the values are and , respectively, while the range is found in Gondek-Rosińska et al. (2003). Overall, the comparison of these three values shows a very good agreement with the results presented in Fig. 7 of Ref. Gondek-Rosińska et al. (2003).

In order to understand the rotation properties of the triaxial solutions, we also report quantities such as dimensionless spin and dimensionless angular momentum in Fig. 6- 8. The dimensionless spin as a function of the eccentricity for MIT bag-model and LX model are shown in Fig. 6 and 7, respectively. The top panel of Fig. 8 reports the spin angular momentum as a function of the eccentricity for the LX EOS. Similar as , the angular momentum increases with the eccentricity. The main difference is that the relative positioning of the curves as a function of compactness is reversed when compared with the plots. In other words for a given eccentricity the greatest angular momentum is achieved for the smallest compactness, while the greatest for the largest one. This is true both for axisymmetric and triaxial solutions. Also as we can see from the bottom panel of Fig. 8, more compact objects can reach greater rotational frequencies while less compact can reach larger angular momenta, which can exceed unity. According to Fig. 6 and 7 as well as the bottom panel of Fig. 8, triaxial sequences lose angular velocity and gain spin angular momentum as one moves towards the mass shedding limit.

An obviously interesting property of triaxially rotating compact stars is that they can act as strong sources of GWs. A full general-relativistic evolution needs to be employed in order to determine accurately the details of the GW emission from such triaxially rotating stars and this is beyond the scope of this paper (see however Ref. Tsokaros et al. (2017) for the case of NSs). At the same time, we can apply the quadrupole formula to make reasonable estimates using the quasi-equilibrium initial data we have computed. The relationship between the normalized GW strain and the eccentricity of the star is shown in the top panel of Fig. 9. Compared with the results of triaxially rotating NSs calculated in Tsokaros et al. (2017), we find that the GW strain for QSs are several times larger for similar values of the compactness. For example, model G4C025 in Tsokaros et al. (2017) with radiates GW with normalized strain 0.007357, while the corresponding amplitude for both the MIT bag-model EOS and the LX EOS is around 0.025 with same eccentricity (see Fig. 9). Also shown in Fig. 9 for the two EOSs considered, are the relations between the strain and the eccentricity, which are very similar and both essentially linear.

## V Triaxial supramassive solutions

Besides the constant rest-mass sequences mentioned above, we have also
built sequences with constant central rest-mass density for both the MIT
bag-model EOS and the LX EOS. We recall that when fixing the central
rest-mass density, the mass of the solutions will increase as the axis
ratio decreases. Furthermore, since we do not impose
axisymmetry, the triaxial deformation will be spontaneously triggered
when is large enough. Therefore, with such calculations we can
determine whether triaxial supramassive QSs, i.e., triaxial solutions with
ADM mass larger than the TOV maximum mass () exist and
the properties of such solutions^{5}^{5}5We recall that for NSs, a
universal relation has been found between and the
maximum mass that can be sustained by axisymmetric solutions in uniform
rotation, (see also Weih et al. (2018) for the case of
differentially rotating stars). More specifically, Breu and Rezzolla
Breu and Rezzolla (2016) found that for a large class of EOSs; we expect a similar
universal behaviour to be present also for QSs,
although the scaling between and
is likely to be different. (Note that all models
shown in Figs. 4–8 are not
supramassive).

According to Uryū et al. (2016), triaxial supramassive NS does exist for the case with polytropic EOS in the range . Furthermore, for the case with a two segments piecewise polytropic EOS, sequences of triaxial supramassive NS solutions become longer, and hence the existence of supramassive triaxial NS becomes evident, when the EOS of the lower density region is stiff (), and the higher density region is soft ()Uryū et al. (2016). Therefore, it is likely that the triaxial supramassive QS also exists because the QS EOS used in this paper has an analogous property, namely, the effective is smaller (softer) in the higher density region, and larger (stiffer) in the lower density region (for the MIT model, see Gourgoulhon et al. (1999)). Besides having an interest of its own, determining the existence of such solutions could be relevant to establish whether a BNS merger could lead to the formation of such an object. Based on current mass measurement constraints Demorest et al. (2010); Antoniadis et al. (2013) and on known BNS systems, the mass of the post-merger product will very likely be larger than .

In order to study this, we have fixed the central rest-mass density close to the value corresponding to and built rotating solution sequences for both the MIT bag-model EOS and the LX EOS, respectively. In this way, we were indeed able to find triaxial supramassive solutions for both EOSs, reporting in Table.6 the solutions with largest triaxial deformation, i.e., the smallest ratio .

Finally, we note that although such models have large compactnesses and we are aware that the IWM formalism becomes increasingly inaccurate for large compactnesses (i.e., with ), we also believe that the associated errors will not change the qualitative result, namely, that triaxial supramassive QS models exist for the EOSs considered here. At the same time, we plan to re-investigate this point in the future, when more accurate methods, such as waveless formulation, will be employed to compute QS solutions.

EOS | ||||||||||
---|---|---|---|---|---|---|---|---|---|---|

MIT | 0.4375 (0.4713) | 0.7657 (0.7938) | 9.978 (16.32) | 0.03870 | 2.862 | 6.847 | 0.1839 | 173.1 | 2.217 | |

LX | 0.4375 (0.4912) | 0.7586 (0.8104) | 7.660 (16.49) | 0.05001 | 3.727 | 11.30 | 0.1948 | 222.1 | 3.325 |

## Vi Conclusion

We have presented a new version of the cocal code to compute axisymmetric and triaxial solutions of uniformly rotating QSs in general relativity with two EOSs, i.e., the MIT bag-model and the LX EOS. Comparisons have been made with NSs as well. Overall, three main properties are found when comparing solution sequences of QSs with of NSs. Firstly, QSs generally have a longer triaxial sequences of solutions than NSs. In another words, QSs can reach a larger triaxial deformation (or smaller ratio) before terminating the sequence at the mass-shedding limit; this is mostly due to the larger ratio that can be attained by QSs. Secondly, when considering similar triaxial configurations, QSs are (slightly) more efficient GW sources; this is mostly due to the finite surface rest-mass density and hence larger mass quadrupole for QSs. Thirdly, triaxial supramassive solutions can be found for QSs as well; this is due again to the fact that larger values of the ratio can be sustained before reaching the mass-shedding limit.

Besides having an interest of its own within solutions of self-gravitating objects in general relativity, triaxially rotating compact stars are important sources for ground-based GW observatories. Our calculations have shown that for rotating QSs with different EOSs, the bifurcation point to triaxial sequence happens at a spin period of , so that the corresponding GW frequency is and hence within the band of GW observatories such as Advanced LIGO or Virgo. Indeed, exploiting the largest triaxial deformation solution obtained in our calculations, the GW strain amplitude can be as large as at a distance of .

Although this is an interesting prospect, it is still unclear whether such triaxial configurations can be produced in practice, since the radiation-reaction timescales needed for the triggering of the secular triaxial instability are still very uncertain, as are the other mechanisms that could contrast the instability. For example, if the triaxial deformation is induced in an isolated star, e.g., a newly born fast rotating star, GW radiation may take away the excess angular momentum very rapidly so that the star would go back to the axisymmetric sequence again after the ratio drops below the critical value. Similarly, when considering stars in binary systems, there is the prospect that an accreting system, such as the one spinning up pulsars, could drive the accreting compact star to exceed the critical ratio, hence leading to a break of axisymmetry. In this process, which is also known as forced GW emission, the triaxial deformation can be maintained via the angular momentum supplied by the accreted matter. Notwithstanding the large uncertainties involved with the details of this picture, such as the presence or not of bifurcation point or the realistic degree of deformation attained by the unstable stars, the fact that these details depend sensitively on the EOS Uryū et al. (2016), suggests that a detection of this type of signal could serve as an important probe for distinguishing the EOS of compact stars.

Finally, we note that the triaxial configurations could also be invoked to explain the spin-up limit for rotating compact stars, which is far smaller than the mass-shedding limit. The results presented here and in Ref. Huang et al. (2008) suggest that when triaxial deformations are taken into account, the rotational period of a compact star actually decreases as it gains angular momentum, e.g., by accretion, along the triaxial sequence. As a result, the “spin-up” process provided by the accretion of matter onto the pulsar can actually spin down the pulsar if the bifurcation point is reached. In this case, no accreting pulsar could spin up faster than the period at the bifurcation point. Of course, depending on the microphysical properties of the QS (e.g., the magnitude of the shear viscosity or of the breaking strain in the crust) it is also possible that other mechanisms of emission of gravitational waves, (e.g., other dynamical instabilities such as the barmode instability Watts et al. (2005); Corvino et al. (2010) or the -mode instability Andersson (1998); Friedman and Morsink (1998); Andersson et al. (1999), or nonzero ellipticities) could intervene at lower spinning frequencies and therefore before the onset of instability to a triaxial deformation is reached Bildsten (1998); Andersson et al. (2011). As a result, the search for fast spinning pulsars with more powerful radio telescopes, such as the Square Kilometre Array (SKA) and the Five-hundred-meter Aperture Spherical radio Telescope (FAST) Kramer and Stappers (2015); Nan et al. (2011) could provide important clues on the properties of pulsars and test the validity of the solid QS assumption Xu (2003).

###### Acknowledgements.

It is a pleasure to thank J. G. Lu, J. Papenfort, Z. Younsi, and all the members of the Pulsar group in Peking University and the Relastro group in Frankfurt for useful discussions. We will thank Dr. L. Shao and Dr. Y. Hu for their help with GW experiments. E. Z. is grateful to the China Scholarship Council for supporting the joint PhD training in Frankfurt. This research is supported in part by the ERC synergy grant ”BlackHoleCam: Imaging the Event Horizon of Black Holes” (Grant No. 610058), by ”NewCompStar”, COST Action MP1304, by the LOEWE-Program in the Helmholtz International Center (HIC) for FAIR, by the European Union’s Horizon 2020 Research and Innovation Programme (Grant 671698) (call FETHPC-1-2014, project ExaHyPE). This work is supported by National Key R&D Program (No.2017YFA0402600) and NNSF (11673002,U1531243). A. T. is supported by NSF Grants PHY-1662211 and PHY-1602536, and NASA Grant 80NSSC17K0070. K. U. is supported by JSPS Grant-in-Aid for Scientific Research(C) 15K05085.## References

- The LIGO Scientific Collaboration and The Virgo Collaboration (2017) The LIGO Scientific Collaboration and The Virgo Collaboration (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. Lett. 119, 161101 (2017).
- The LIGO Scientific Collaboration et al. (2017a) The LIGO Scientific Collaboration, the Virgo Collaboration, B. P. Abbott, R. Abbott, T. D. Abbott, F. Acernese, K. Ackley, C. aAdams, T. Adams, P. Addesso, and et al. (LIGO Scientific Collaboration and Virgo Collaboration), Astrophys. J. Lett. 848, L12 (2017a).
- LIGO Scientific Collaboration et al. (2017) LIGO Scientific Collaboration, Virgo Collaboration, F. Gamma-Ray Burst Monitor, and INTEGRAL, Astrophys. J. Lett. 848, L13 (2017), arXiv:1710.05834 [astro-ph.HE] .
- The LIGO Scientific Collaboration et al. (2017b) The LIGO Scientific Collaboration, the Virgo Collaboration, B. P. Abbott, R. Abbott, T. D. Abbott, F. Acernese, K. Ackley, C. Adams, T. Adams, P. Addesso, and et al., ArXiv e-prints (2017b), arXiv:1710.05836 [astro-ph.HE] .
- Baiotti and Rezzolla (2017) L. Baiotti and L. Rezzolla, Rept. Prog. Phys. 80, 096901 (2017), arXiv:1607.03540 [gr-qc] .
- Andersson et al. (2011) N. Andersson, V. Ferrari, D. I. Jones, K. D. Kokkotas, B. Krishnan, J. S. Read, L. Rezzolla, and B. Zink, General Relativity and Gravitation 43, 409 (2011), arXiv:0912.0384 [astro-ph.SR] .
- Abramovici et al. (1992) A. A. Abramovici, W. Althouse, R. P. Drever, Y. Gursel, S. Kawamura, F. Raab, D. Shoemaker, L. Sievers, R. Spero, K. S. Thorne, R. Vogt, R. Weiss, S. Whitcomb, and M. Zuker, Science 256, 325 (1992).
- Punturo et al. (2010) M. Punturo et al., Class. Quantum Grav. 27, 084007 (2010).
- Accadia et al. (2011) T. Accadia et al., Class. Quantum Grav. 28, 114002 (2011).
- Kuroda and LCGT Collaboration (2010) K. Kuroda and LCGT Collaboration, Class. Quantum Grav. 27, 084004 (2010).
- Aso et al. (2013) Y. Aso, Y. Michimura, K. Somiya, M. Ando, O. Miyakawa, T. Sekiguchi, D. Tatsumi, and H. Yamamoto, Phys. Rev. D 88, 043007 (2013), arXiv:1306.6747 [gr-qc] .
- Chandrasekhar (1969) S. Chandrasekhar, The Silliman Foundation Lectures, New Haven: Yale University Press, 1969 (1969).
- Meinel et al. (2008) R. Meinel, M. Ansorg, A. Kleinwächter, G. Neugebauer, and D. Petroff, Relativistic Figures of Equilibrium (Cambridge University Press, 2008).
- Friedman and Stergioulas (2013) J. L. Friedman and N. Stergioulas, Rotating Relativistic Stars, by John L. Friedman , Nikolaos Stergioulas, Cambridge, UK: Cambridge University Press, 2013 (2013).
- Shapiro and Zane (1998) S. L. Shapiro and S. Zane, Astrop. J. Supp. 117, 531 (1998).
- Bonazzola et al. (1996) S. Bonazzola, J. Frieben, and E. Gourgoulhon, Astrophys. J. 460, 379 (1996), gr-qc/9509023 .
- Gondek-Rosińska and Gourgoulhon (2002) D. Gondek-Rosińska and E. Gourgoulhon, Phys. Rev. D 66, 044021 (2002), gr-qc/0205102 .
- Bonazzola et al. (1998) S. Bonazzola, J. Frieben, and E. Gourgoulhon, Astron. Astrophys. 331, 280 (1998), gr-qc/9710121 .
- Houser and Centrella (1996) J. L. Houser and J. M. Centrella, Phys. Rev. D 54, 7278 (1996), gr-qc/9611033 .
- Pickett et al. (1996) B. K. Pickett, R. H. Durisen, and G. A. Davis, Astrophys. J. 458, 714 (1996).
- Brown (2000) J. D. Brown, Phys. Rev. D 62, 084024 (2000), gr-qc/0004002 .
- New et al. (2000) K. C. B. New, J. M. Centrella, and J. E. Tohline, Phys. Rev. D 62, 064019 (2000), arXiv:astro-ph/9911525 .
- Liu (2002) Y. T. Liu, Phys. Rev. D 65, 124003 (2002), arXiv:gr-qc/0109078 .
- Watts et al. (2005) A. L. Watts, N. Andersson, and D. I. Jones, Astrophys. J. 618, L37 (2005), astro-ph/0309554 .
- Baiotti et al. (2007) L. Baiotti, R. De Pietri, G. M. Manca, and L. Rezzolla, Phys. Rev. D 75, 044023 (2007), astro-ph/0609473 .
- Manca et al. (2007) G. M. Manca, L. Baiotti, R. DePietri, and L. Rezzolla, Class. Quantum Grav. 24, S171 (2007), arXiv:0705.1826 .
- Corvino et al. (2010) G. Corvino, L. Rezzolla, S. Bernuzzi, R. De Pietri, and B. Giacomazzo, Class. Quantum Grav. 27, 114104 (2010), arXiv:1001.5281 [gr-qc] .
- Camarda et al. (2009) K. D. Camarda, P. Anninos, P. C. Fragile, and J. A. Font, Astrophys. J. 707, 1610 (2009), arXiv:0911.0670 [astro-ph.SR] .
- Franci et al. (2013) L. Franci, R. De Pietri, K. Dionysopoulou, and L. Rezzolla, Journal of Physics Conference Series 470, 012008 (2013), arXiv:1309.6549 [gr-qc] .
- Muhlberger et al. (2014) C. D. Muhlberger, F. H. Nouri, M. D. Duez, F. Foucart, L. E. Kidder, C. D. Ott, M. A. Scheel, B. Szilágyi, and S. A. Teukolsky, Phys. Rev. D 90, 104014 (2014), arXiv:1405.2144 [astro-ph.HE] .
- Andersson (2003) N. Andersson, Class. Quantum Grav. 20, 105 (2003), arXiv:astro-ph/0211057 .
- Lai and Shapiro (1995) D. Lai and S. L. Shapiro, Astrophys. J. 442, 259 (1995), astro-ph/9408053 .
- Bildsten (1998) L. Bildsten, Astrophys. J. L 501, L89 (1998), astro-ph/9804325 .
- Woosley and Janka (2005) S. Woosley and T. Janka, Nature Physics 1, 147 (2005), astro-ph/0601261 .
- Watts et al. (2008) A. L. Watts, B. Krishnan, L. Bildsten, and B. F. Schutz, Mon. Not. R. Astron. Soc. 389, 839 (2008), arXiv:0803.4097 .
- Piro and Thrane (2012) A. L. Piro and E. Thrane, Astrophys. J. 761, 63 (2012), arXiv:1207.3805 [astro-ph.HE] .
- Huang et al. (2008) X. Huang, C. Markakis, N. Sugiyama, and K. Uryū, Phys. Rev. D78, 124023 (2008), arXiv:0809.0673 [astro-ph] .
- Uryū et al. (2016) K. Uryū, A. Tsokaros, F. Galeazzi, H. Hotta, M. Sugimura, K. Taniguchi, and S. Yoshida, Phys. Rev. D93, 044056 (2016).
- James (1964) R. A. James, Astrophys. J. 140, 552 (1964).
- Hachisu and Eriguchi (1982) I. Hachisu and Y. Eriguchi, Progress of Theoretical Physics 68, 206 (1982).
- Lai et al. (1993) D. Lai, F. A. Rasio, and S. L. Shapiro, Astrophys. J., Suppl. Ser. 88, 205 (1993).
- Bodmer (1971) A. R. Bodmer, Phys. Rev. D 4, 1601 (1971).
- Witten (1984) E. Witten, Phys. Rev. D 30, 272 (1984).
- Dai et al. (2016) Z. G. Dai, S. Q. Wang, J. S. Wang, L. J. Wang, and Y. W. Yu, Astrophys. J. 817, 132 (2016), arXiv:1508.07745 [astro-ph.HE] .
- Lai et al. (2017) X. Y. Lai, Y. W. Yu, E. P. Zhou, Y. Y. Li, and R. X. Xu, Research in Astronomy and Astrophysics , in press (2017), arXiv:1710.04964 [astro-ph.HE] .
- Li et al. (2016) A. Li, B. Zhang, N.-B. Zhang, H. Gao, B. Qi, and T. Liu, Phys. Rev. D 94, 083010 (2016), arXiv:1606.02934 [astro-ph.HE] .
- Itoh (1970) N. Itoh, Progress of Theoretical Physics 44, 291 (1970).
- Alcock et al. (1986) C. Alcock, E. Farhi, and A. Olinto, Astrophys. J. 310, 261 (1986).
- Haensel et al. (1986) P. Haensel, J. L. Zdunik, and R. Schaefer, Astron. Astrophys. 160, 121 (1986).
- Gondek-Rosińska et al. (2000) D. Gondek-Rosińska, T. Bulik, L. Zdunik, E. Gourgoulhon, S. Ray, J. Dey, and M. Dey, Astron. Astrophys. 363, 1005 (2000), astro-ph/0007004 .
- Gourgoulhon et al. (1999) E. Gourgoulhon, P. Haensel, R. Livine, E. Paluch, S. Bonazzola, and J.-A. Marck, Astron. Astrophys. 349, 851 (1999), astro-ph/9907225 .
- Stergioulas et al. (1999) N. Stergioulas, W. Kluźniak, and T. Bulik, Astron. Astrophys. 352, L116 (1999), astro-ph/9909152 .
- Szkudlarek et al. (2012) M. Szkudlarek, Gondek-Rosiń, D. ska, L. Villain, and M. Ansorg, in Electromagnetic Radiation from Pulsars and Magnetars, Astronomical Society of the Pacific Conference Series, Vol. 466, edited by W. Lewandowski, O. Maron, and J. Kijak (2012) p. 231.
- Gondek-Rosinska et al. (2000) D. Gondek-Rosinska, P. Haensel, J. L. Zdunik, and E. Gourgoulhon, in IAU Colloq. 177: Pulsar Astronomy - 2000 and Beyond, Astronomical Society of the Pacific Conference Series, Vol. 202, edited by M. Kramer, N. Wex, and R. Wielebinski (2000) p. 661, astro-ph/0009282 .
- Gondek-Rosińska et al. (2001) D. Gondek-Rosińska, N. Stergioulas, T. Bulik, W. Kluźniak, and E. Gourgoulhon, Astron. Astrophys. 380, 190 (2001), astro-ph/0110209 .
- Gondek-Rosińska et al. (2003) D. Gondek-Rosińska, E. Gourgoulhon, and P. Haensel, Astron. Astrophys. 412, 777 (2003), astro-ph/0311128 .
- Uryū and Tsokaros (2012) K. Uryū and A. Tsokaros, Phys. Rev. D 85, 064014 (2012), arXiv:1108.3065 [gr-qc] .
- Uryū et al. (2012) K. Uryū, A. Tsokaros, and P. Grandclement, Phys. Rev. D86, 104001 (2012), arXiv:1210.5811 [gr-qc] .
- Tsokaros and Uryū (2012) A. Tsokaros and K. Uryū, Journal of Engineering Mathematics 82, 133 (2012).
- Tsokaros et al. (2015) A. Tsokaros, K. Uryū, and L. Rezzolla, Phys. Rev. D 91, 104030 (2015), arXiv:1502.05674 [gr-qc] .
- Tsokaros et al. (2017) A. Tsokaros, M. Ruiz, V. Paschalidis, S. L. Shapiro, L. Baiotti, and K. Uryū, Phys. Rev. D 95, 124057 (2017), arXiv:1704.00038 [gr-qc] .
- Rezzolla and Zanotti (2013) L. Rezzolla and O. Zanotti, Relativistic Hydrodynamics (Oxford University Press, Oxford, UK, 2013).
- Isenberg and Nester (1980) J. Isenberg and J. Nester, in General Relativity and Gravitation. Vol. 1. One hundred years after the birth of Albert Einstein. Edited by A. Held. New York, NY: Plenum Press, p. 23, 1980, Vol. 1, edited by A. Held (1980) p. 23.
- Isenberg (2008) J. A. Isenberg, International Journal of Modern Physics D 17, 265 (2008), arXiv:gr-qc/0702113 .
- Wilson and Mathews (1989) J. R. Wilson and G. J. Mathews, “Relativistic hydrodynamics.” in Frontiers in Numerical Relativity, edited by C. R. Evans, L. S. Finn, and D. W. Hobill (1989) pp. 306–314.
- Chodos et al. (1974) A. Chodos, R. L. Jaffe, K. Johnson, C. B. Thorn, and V. F. Weisskopf, Phys. Rev. D 9, 3471 (1974).
- Limousin et al. (2005) F. Limousin, D. Gondek-Rosińska, and E. Gourgoulhon, Phys. Rev. D 71, 064012 (2005), gr-qc/0411127 .
- Lai and Xu (2009) X. Y. Lai and R. X. Xu, Mon. Not. R. Astron. Soc. 398, L31 (2009), arXiv:0905.2839 [astro-ph.HE] .
- Lai and Xu (2017) X. Y. Lai and R. X. Xu, in Journal of Physics Conference Series, Journal of Physics Conference Series, Vol. 861 (2017) p. 012027, arXiv:1701.08463 [astro-ph.HE] .
- Guo et al. (2014) Y.-J. Guo, X.-Y. Lai, and R.-X. Xu, Chinese Physics C 38, 055101 (2014).
- Lai and Xu (2011) X.-Y. Lai and R.-X. Xu, Research in Astronomy and Astrophysics 11, 687 (2011), arXiv:1011.0526 [astro-ph.SR] .
- Demorest et al. (2010) P. B. Demorest, T. Pennucci, S. M. Ransom, M. S. E. Roberts, and J. W. T. Hessels, Nature 467, 1081 (2010), arXiv:1010.5788 [astro-ph.HE] .
- Antoniadis et al. (2013) J. Antoniadis, P. C. C. Freire, N. Wex, T. M. Tauris, R. S. Lynch, M. H. van Kerkwijk, M. Kramer, C. Bassa, V. S. Dhillon, T. Driebe, J. W. T. Hessels, V. M. Kaspi, V. I. Kondratiev, N. Langer, T. R. Marsh, M. A. McLaughlin, T. T. Pennucci, S. M. Ransom, I. H. Stairs, J. van Leeuwen, J. P. W. Verbiest, and D. G. Whelan, Science 340, 448 (2013), arXiv:1304.6875 [astro-ph.HE] .
- Zhou et al. (2014) E. P. Zhou, J. G. Lu, H. Tong, and R. X. Xu, Mon. Not. R. Astron. Soc. 443, 2705 (2014), arXiv:1404.2793 [astro-ph.HE] .
- Xu and Liang (2009) R. Xu and E. Liang, Science in China: Physics, Mechanics and Astronomy 52, 315 (2009), arXiv:0811.0234 .
- Wang et al. (2017) W. Wang, J. Lu, H. Tong, M. Ge, Z. Li, Y. Men, and R. Xu, Astrophys. J. 837, 81 (2017), arXiv:1603.08288 [astro-ph.HE] .
- Xu (2003) R. X. Xu, Astrophys. J. Lett. 596, L59 (2003), astro-ph/0302165 .
- Read et al. (2009) J. S. Read, B. D. Lackey, B. J. Owen, and J. L. Friedman, Phys. Rev. D 79, 124032 (2009), arXiv:0812.2163 .
- Li et al. (2017) A. Li, Z.-Y. Zhu, and X. Zhou, Astrophys. J. 844, 41 (2017).
- Bhattacharyya et al. (2016) S. Bhattacharyya, I. Bombaci, D. Logoteta, and A. V. Thampan, Mon. Not. R. Astron. Soc. 457, 3101 (2016), arXiv:1601.06120 [astro-ph.HE] .
- Alford et al. (2005) M. Alford, M. Braby, M. Paris, and S. Reddy, Astrophys. J. 629, 969 (2005), nucl-th/0411016 .
- Baym et al. (2017) G. Baym, T. Hatsuda, T. Kojo, P. D. Powell, Y. Song, and T. Takatsuka, ArXiv e-prints (2017), arXiv:1707.04966 [astro-ph.HE] .
- Ostriker and Mark (1968) J. P. Ostriker and J. W.-K. Mark, Astrophys. J. 151, 1075 (1968).
- Komatsu et al. (1989) H. Komatsu, Y. Eriguchi, and I. Hachisu, Mon. Not. R. Astron. Soc. 237, 355 (1989).
- Tsokaros et al. (2016) A. Tsokaros, B. C. Mundim, F. Galeazzi, L. Rezzolla, and K. Uryū, arXiv:1605.07205 (2016), arXiv:1605.07205 [gr-qc] .
- Tsokaros and Uryū (2007) A. A. Tsokaros and K. Uryū, Phys. Rev. D 75, 044026 (2007).
- Huang et al. (2007) L. Huang, M. Cai, Z.-Q. Shen, and F. Yuan, Mon. Not. R. Astron. Soc. 379, 833 (2007), astro-ph/0703254 .
- Weih et al. (2018) L. R. Weih, E. R. Most, and L. Rezzolla, Mon. Not. R. Astron. Soc. 473, L126 (2018), arXiv:1709.06058 [gr-qc] .
- Breu and Rezzolla (2016) C. Breu and L. Rezzolla, Mon. Not. R. Astron. Soc. 459, 646 (2016), arXiv:1601.06083 [gr-qc] .
- Uryū et al. (2016) K. Uryū, A. Tsokaros, L. Baiotti, F. Galeazzi, N. Sugiyama, K. Taniguchi, and S. Yoshida, Phys. Rev. D 94, 101302 (2016), arXiv:1606.04604 [astro-ph.HE] .
- Andersson (1998) N. Andersson, Astrophys. J. 502, 708 (1998), arXiv:gr-qc/9706075 .
- Friedman and Morsink (1998) J. L. Friedman and S. M. Morsink, Astrophys. J. 502, 714 (1998), arXiv:gr-qc/9706073 .
- Andersson et al. (1999) N. Andersson, K. D. Kokkotas, and N. Stergioulas, Astrophys. J. 510, 854 (1999).
- Kramer and Stappers (2015) M. Kramer and B. Stappers, Advancing Astrophysics with the Square Kilometre Array (AASKA14) , 36 (2015), arXiv:1507.04423 [astro-ph.IM] .
- Nan et al. (2011) R. Nan, D. Li, C. Jin, Q. Wang, L. Zhu, W. Zhu, H. Zhang, Y. Yue, and L. Qian, International Journal of Modern Physics D 20, 989 (2011), arXiv:1105.3794 [astro-ph.IM] .