Evolution of bare quark stars in full general relativity: Single star case

We introduce our approaches, in particular the modifications of the primitive recovery procedure, to handle bare quark stars in numerical relativity simulations. Reliability and convergence of our implementation are demonstrated by evolving two triaxially rotating quark star models with different mass as well as a differentially rotating quark star model which has sufficiently large kinetic energy to be dynamically unstable. These simulations allow us to verify that our method is capable of resolving the evolution of the discontinuous surface of quark stars and possible mass ejection from them. The evolution of the triaxial deformation and the properties of the gravitational-wave emission from triaxially rotating quark stars have been also studied, together with the mass ejection of the differentially rotating case. It is found that supramassive quark stars are not likely to be ideal sources of continuous gravitational wave as the star recovers axisymmetry much faster than models with smaller mass and gravitational-wave amplitude decays rapidly in a timescale of $10\,$ms, although the instantaneous amplitude from more massive models is larger. As with the differentially rotating case, our result confirms that quark stars could experience non-axisymmetric instabilities similar to the neutron star case but with quite small degree of differential rotation, which is expected according to previous initial data studies.


I. INTRODUCTION
The detection of gravitational-wave (GW) signals from binary neutron star (BNS) merger events (GW170817 and GW190425) [1,2], as well as the electromagnetic (EM) counterparts [3], have provided significant insights in the equation of state (EoS) of neutron stars (NSs). On the one hand, the inspiral GW signal puts a constraint on the tidal deformability of NSs [4][5][6][7], which translates into a constraint on the radius of NSs. On the other hand, the EM counterparts indicate the fate of the merger remnant, which is tightly related to the maximum mass of cold and non-spinning NSs (or Tolman-Oppenheimer-Volkov maximum mass, M TOV ) and the total mass of the binary. Various constraints have thus been put on the EoS of NSs accordingly [8][9][10][11][12].
Nevertheless, the state of matter at density as high as in NSs, particularly the role of quark matter, is still in open debate. Recent studies have revealed that the existence of deconfined quark matter is possible in the high density core part of massive NSs, confronted with the most recent astrophysical observations [13]. Many attempts aimed at searching for evidence of such a strong interaction phase transition (from conventional hadronic matter to quark matter) inside NSs from astrophysical observations have also been made [14,15]. As another possibility, 3-flavor quark matter could be the absolutely stable state of matter at high density [16,17]. This may result in the existence of the so-called strange quark star (QS). Such object could exist above a certain NS mass (i.e., the twofamily scenario [18,19]), and hence could be formed during the merger of BNSs. Alternatively, there are also suggestions that QSs could exist even for the normal NS mass range [20] and following this assumption the binary QS (BQS) merger scenario has also been studied in terms of the tidal deformability and mass ejections [21][22][23].
Figuring out the differences between the merger events involving QSs and NSs and verifying them with multi-messenger observations could greatly enrich our understanding in the nature of strong interaction. This requires the modeling of QSs with the tool of general relativistic hydrodynamics as a first step. Previously, various efforts have been made in the calculation of the quasi-/equilibrium structures of uniformly and differentially rotating QSs in general relativity (GR) [24][25][26][27][28][29] as well as the configuration of BQSs in the last orbits [30]. Despite all these attempts in constructing initial data for QSs, the progress in dynamical evolution of them is quite limited. In fact, there has been only one successful simulation of BQSs [22,23], in which only approximate GR and the smooth particle hydrodynamics method are adopted. It is both necessary and helpful for us to improve our understanding by pushing forward to evolve QSs in full GR grid-based hydrodynamics, which has yet been performed, yet.
The problem in directly evolving QSs results from the structural difference between QSs and NSs. QSs are selfbound by strong interaction and the surface density is nonzero (hence also referred to as a 'bare' QS). Handling this discontinuity on the surface of QSs is a challenge in numerical simulation. In addition, the specific enthalpy of QSs calculated in the conventional way does not go to unity as pressure goes to zero [28], and it may cause a problem when we recover fundamental thermodynamic variables from conserved variables which are directly evolved in numerical relativity (NR) (the so-called primitive recovery procedure). Realizing this, we have tried to modify the NR code SACRA-MPI [31,32] to resolve these two problems and achieved successful evolution of single QSs. In this paper, we report the details of our methods to handle QSs in NR and the first results of dynamical evolution of two triaxially rotating QS models and a differentially rotating QS model.
There are several motivations for us to tackle the problem of evolving triaxially and differentially rotating QSs as a first step. First, these systems provide us a good test-bed problem for evolving BQSs, for which we have to handle the mo-tion of the sharp density discontinuity at the surface. The non-zero density surface of the triaxial QSs naturally moves in the computational domain during the evolution. This is also the case for differentially rotating stars with sufficiently large values of T /|W | ratio (the ratio of kinetic energy, T , to gravitational potential energy, W ) to be unstable for dynamical bar modes [33,34]. Thus, an essential technical element needed to handle the merger of BQSs can be tested by the simulation for these systems. Secondly, previous studies indicated that triaxial deformation could play a more important role for QSs than for NSs, as it could reach higher values of T /|W | [28,35]. Specifically, supramassive triaxial QSs (i.e., those with mass larger than M TOV ) could exist for a large parameter space, while normally no such solution for NSs is present unless the EoS is extremely stiff in the crust and softer in the core [28,36]. Therefore, a newly-born QS in supernova or binary merger events may be a source for ground based GW detectors and is worthy to be investigated. Thirdly, determining the properties of mass ejection from binary merger events is a key for understanding the EM counterparts and we would like to make sure that our approach is capable of capturing the ejected mass. Mass ejection is expected from differentially rotating compact stars during the growth of the dynamical barmode instability [34] if T /|W | is large enough. Thus we can confirm the ability of our methods by evolving such unstable differentially rotating configurations.
The paper is organized as follows: In Sec. II we describe the EoS model used in this work and how we treat it in our numerical implementation. Details about the setup of our initial data solver and dynamical evolution code as well as the convergence tests are presented in Sec. III. The results for the evolution of the QS models, as well as the physical interpretation of the results are reported in Sec. IV. In the end, we briefly conclude and discuss future prospects of evolving QSs in Sec. V. More quantitative tests on the code performance are found in the Appendix. Throughout this paper, c and G denote the speed of light and the gravitational constant, respectively.

II. QUARK STAR EQUATION OF STATE
The most widely used EoS for QSs is the MIT bag model, which is initially used as a phenomenological model for hadrons [37]. Modern versions of the model take into account more nuclear physics details, such as perturbative quantum choromodynamics (QCD) corrections due to the gluonmediated interaction between quarks as well as the finite mass of strange quarks [38,39]. As our first attempt to evolve QSs, we neglect these details and stick with the simplest form of the MIT bag model, which could be written as: where p is the pressure for given energy density e, and B is the bag constant which is related to the finite surface density. Ideally, the pressure of quark matter vanishes at e = 4B and density drops to 0 discontinuously across the surface; that is, Eq. (1) is valid only for e ≥ 4B. However, matter with density lower than the surface density exists in nature. Even for the case of 'bare' QS (which is exactly the model considered in this paper) that is not supposed to possess a low density crust consisting of conventional hadronic matter, such low density matter is still relevant when mass ejection is present (for instance, in the cases of dynamically unstable models and BQS mergers). Theoretically, the exact form and fate of the material ejected from a BQS merger are unclear. Some previous studies suggested that the ejection from QSs is likely to be in the form of strange quark nuggets [40,41], which are also self-bound as the bare QS itself (i.e., the pressure is also zero on the surface of the ejected nuggets). The main source of pressure from the low density ejecta of QSs is therefore the thermal contribution. Nevertheless, more recent studies indicate that the possibility of ejected quark matter decaying into normal nucleon matter does exist [42,43]. In this paper, we will not consider more complicated nuclear physical process at the moment and model such matter as ideal gas (the details of which is found in Sec. III B).
To obtain the rest-mass density and specific enthalpy which are variables necessary for relativistic hydrodynamics, we need the more explicit form of the MIT bag model as where ρ could be regarded as the rest-mass density 1 and K is constant. The choice of K is usually determined by details of interactions between quarks in the system [38,39,44,45], in which a 4 is a model parameter that characterizes the perturbative QCD corrections. However, as discussed in many previous studies [28,44,46], the choice of a 4 (hence K) has only negligible influence on the EoS model once the bag constant is fixed. In particular, in the case that we neglect the mass of strange quarks (which is exactly the case assumed in this work), the choice of K has no impact on the EoS model at all. This could be easily understood if one eliminates ρ from Eq. (2) to recover Eq. (1), regardless of whatever value for K is used. Consequently, we could employ a value of K which is most convenient for us to implement numerically, and ρ becomes a purely auxiliary variable 2 .
With the expression of Eq. (2), the specific enthalpy is written as On the surface of the QSs, the density ρ s is determined by the fact that the pressure p s should vanish:  [47,48] as well as the upper limit of tidal deformability [1].
which reads Hence the specific enthalpy on the surface of QSs in this model is Now it becomes obvious that normally with a choice of K that is motivated by nuclear physics, h s does not approach c 2 . This could bring another source of discontinuity, because for the case of the normal matter for which the internal energy is much smaller than the rest-mass energy, the specific enthalpy approaches c 2 . In particular, most of the physically motivated choice of K for the quark matter results in h s < c 2 . This implies that the specific enthalpy is not only discontinuous across the surface but also not monotonic, and ρ as a function of h is not even a single-valued function. This could make the primitive recovery procedure extremely complicated, if ever possible.
Taking into account the fact that changing K does not matter much to the EoS modeling for our purpose, it is then quite straightforward to make a choice for K such that h s = c 2 , which requires This implies that for every choice of the bag constant, we could determine a unique choice for Eq. (2) that resolves the issue of an ill-behaved function of h(ρ). The EoS parameters adopted in this work and also properties of the spherical QSs are listed in Table I and the mass-radius relation is shown in Fig. 1.
Although shock heating does not play an important role for the evolution of a single star, we still employ the Γ th prescription to take into account finite temperature effects as we ultimately will try to deal with BQS systems. In this model, the pressure is expressed as a sum of two components, the cold part (p cold (ρ)) determined by Eq. (2) and a thermal part related to the specific internal energy density as where th is the thermal part of specific internal energy density defined by Here we follow the normal definition such that specific internal energy density = e/ρ − c 2 . It is worth noting that only with the choice of K according to Eq. (7), are we able to guarantee that cold approaches zero on the surface of the star. In this work, Γ th is chosen to be 4/3 which makes the primitive recovery process much easier to do for QSs (which we will demonstrate later in Sec. III B), although in principle any other values are also possible to be employed at a higher computational cost.

A. Initial Data
The quasi-equilibrium configurations of rotating QSs are calculated by the initial data solver COCAL (Compact Object CALculator) [49,50]. The implementation of QS EoS models as well as convergence and accuracy tests are found in our previous papers for axisymmetrically uniformly/differentially rotating and triaxially uniformly rotating cases [28,29].
For the triaxial rotation case, as mentioned before, one significant difference between NSs and QSs is that supramassive triaxially rotating QSs could exist for a vast space of parameters. To explore both the supramassive case (which could be formed during a binary merger event) as well as a relatively standard mass case (which could be formed in a supernova), we have built two solution sequences with different values of the central density. The way of constructing the solution sequences is the same as in [28]: We construct the sequence of solutions for a fixed value of the central density with decreasing R z /R x ratio without imposing axisymmetry during the calculation. Here the direction of z-axis corresponds to the angular momentum direction of the rotating star and R z refers to the coordinate radii of the star along z-axis. Since we are considering triaxial configuration here, R x refers to the coordinate radii of the longest axis in the equatorial plane and R y to the shortest. As a consequence, in the presence of the solution with R z /R x small enough to possess a sufficiently high value of the T /|W | ratio for bifurcating into the triaxial solution branch, a star with R y /R x < 1 is spontaneously obtained. We then pick up a triaxial solution with mass that is astrophysically interesting and use it as our initial data.
As with the differential rotation case, the widely used jconst law [51][52][53][54][55][56][57]: is chosen to construct the solutions, in which A is a parameter that characterizes the degree of differential rotation and a normalized versionÂ = A/r e is also often used where r e is the equatorial coordinate radius of the star (for smaller values of A, the differential rotation degree is higher). Previous initialdata studies [27,29] show that the properties of differentially rotating QSs for a given value ofÂ is quite different from that of NSs. In particular, continuous transition to toroidal sequences (i.e., type C solution according to the classification of [52]) happens at lower degree of differential rotation (at A ∼ 3 for QSs while atÂ = 1 for NSs). In previous studies for NSs [34,58], theÂ parameter is typically chosen to be unity or smaller for the case of NSs to explore the dynamical bar-mode instability. Taking the difference between QSs and NSs mentioned above into account, we chooseÂ = 3.0 for the MIT bag model and pick up a solution with T /|W | large enough as the initial data. In addition, we evolve a differentially rotating NS model with the APR4 EoS [59] andÂ = 1.0 for comparison purpose. The quantities estimated for the initial-data models we consider in this work are listed in Table II. In the following, the triaxial model with lower mass will be referred to as MIT148 and the supramassive triaxial model as MIT265, according to their ADM mass. The differential rotation models will be referred to as MIT275dr and APR206dr for the QS and NS cases, respectively.

B. Dynamical Evolution
The dynamical evolution of the initial data is performed with our numerical relativity code SACRA-MPI. To solve Einstein's evolution equation, a moving puncture version of the Baumgarte-Shapiro-Shibata-Nakamura formalism [61][62][63][64] is employed in the code. A constraint propagation prescription similar to Z4c scheme [65] is also implemented in part of the computational domain (cf. [66] for details). More details about SACRA-MPI, such as the finite differencing schemes as well as the adaptive/fixed mesh refinement (AMR/FMR) setups, are found in [32].
In numerical hydrodynamics, we employ a high-resolution shock capturing scheme. For this case, one always has to impose a small but non-zero density for the exterior of the star to recover the 4-velocity (which is a primitive variable) from the momentum density (which is an evolved variable). In this work, this density is set to be 10 −12 times of the surface density of our QS model 3 . Nevertheless, the situation is a bit more complicated than NS cases, even if we manage to choose a model according to Eq. (7) such that the specific enthalpy is now well-behaved across the surface of QSs. For NS cases, Eq. (8) is valid in the entire computational domain, with p cold calculated with exactly the same EoS both for the NS and the region outside the NS. However, it is obvious that this is not the case for QSs, because p cold vanishes for ρ < ρ s . As a consequence, only the thermal component in Eq. (8) exists for the region outside the QSs: and we choose Γ th = 4/3 for ρ < ρ s , same as the thermal contribution of the matter in the QSs such that the EoS model is consistent throughout the inside and outside of the star.
Having to use two different EoS models for the QS and the other brings another problem. In relativistic hydrodynamics, the conserved variables are evolved in every time step, and the primitive variables (such as the rest-mass density, 4-velocity, and specific internal energy) have to be determined from the conserved variables with the help of EoS and normalization relation of the 4-velocity (i.e., by the so-called primitive recovery procedure). However, since we have no idea about the value of the rest-mass density before the primitive recovery, we do not know whether the part of matter whose primitive variables to be recovered is inside the QS or not. As a result, we do not know which EoS model is adopted for the primitive recovery procedure in advance. Therefore, we have to modify the primitive recovery part to incorporate QS EoS models. We will briefly explain our strategy in the following.
For our implementation of relativistic hydrodynamics, the conserved variables which are directly evolved are Model Rx is the coordinate (proper) equatorial radius and Rz/Rx is the ratio of coordinate radius along zand x-axis and Ry/Rx is for yto x-axis ratio. ρm is the maximum rest-mass density inside the star. MADM, M b , J, T /|W | , Pc andÂ are Arnowit-Deset-Misner mass, baryonic mass, angular momentum, ratio between kinetic and gravitational potential energy, the central rotation period, and the differential rotation parameter in the j-const law. Definitions can be found in the Appendix of [60].Ė is the luminosity of the GW at t = 0 as estimated by the quadrupole formula which is a good approximation for the instantaneous GW strain in the beginning of the dynamical evolution.
where √ γ is the determinant of the 3-metric on the spacelike hypersurface, γ ij , and is solved during the time evolution, u i is the spatial component of the 4-velocity of the fluid element, p, ρ, h and w are the total pressure (including the thermal contribution), rest-mass density, specific enthalpy, and Lorentz factor of the fluid, which we would like to recover. Including the three components in u i , we have in total 7 unknown variables and only 5 relations between them (Eqs. (12)- (14)), and hence, additional two relations are needed. For this, one is the EoS model and the other is the normalization relation of the 4-velocity, which is rewritten to the definition of Lorentz factor as where Eq. (13) was used for the second equality. We here define q 2 ≡ γ ij S i S j /ρ 2 * (which is purely determined during the evolution) in the following to simplify the expressions. Another given quantity is defined using Eqs. (12) and (14) as Note that EoS provides us p as a function of w (since ρ * and √ γ are obtained during evolution, ρ can be obtained once w is determined according to Eq. (12)) and h. Therefore, Eqs. (15) and (16) are two equations for two unknowns h and w. Thus, the primitive recovery, in essence, is to find a solution for h and w using the above two algebraic equations. For the region outside the QSs, for which the EoS is Eq. (11), the following algebraic equation for h is obtained by substituting out w from Eqs. (15) and (16): with which we obtain the following equation for w: Here, we defined s = e 0 − B √ γ/ρ * . A solution can be obtained analytically 4 for w 2 as long as s 4 − (3/4)q 2 s 2 is larger than 0. We can then take the root for w 2 which is larger than 1, if exists, and recover the value for ρ to see whether it is indeed larger than ρ s . In the case that the recovered value of ρ is actually smaller than ρ s or there exists no root for w 2 that satisfies w ≥ 1, we treat the fluid as the non-quark matter, and solve Eq. (17) with the Newton-Raphson method to find a solution. We succeeded in evolving the two triaxially rotating and one differentially rotating QS models with this implementation. In particular, it is worth mentioning that for the case of triaxial solutions and dynamically unstable differentially rotating solution, the surface of a QS is moving in the computational domain. This implies that there exists certain grid points on which the fluid is sometimes inside the QS and sometimes outside. The fact that the triaxial solutions can be evolved with this method verifies its reliability. We will focus on the tests of the code performance in the next subsection.

C. Code Tests
We employ three different grid setups in this work to explore the accuracy and convergence behavior of the code in evolving QSs. Nine FMR levels, which are all centered at the QS core, are constructed for all the three grid configurations, with every level doubling the size and grid intervals of its finer level. In every refinement level, uniform and vertex-centered Cartesian coordinates are used to cover the x, y, and z directions with 2N + 1, 2N + 1 and N + 1 points (equatorial plane symmetry is assumed) and we choose N = 80, 120, and 160 for three different setups. The distance from the coordinate origin to the boundary of the finest level along each axis is set to be 16 GM /c 2 (≈ 23.6 km) for all the three setups; i.e., an outer boundary of the entire FMR domain is as large as 16 × 2 8 = 4096 GM /c 2 ≈ 6048 km. The boundary is much larger than the typical GW wavelength of the models we considered here, λ gw ∼ π/Ω, which is of the order of 100 GM /c 2 ≈ 147.7 km according to the rotational period found from the initial data. It is straightforward to find the grid resolution in the finest level as ∆ x = ∆ y = ∆ z = 295, 197, and 148 m for three different resolutions.
We monitor both the volume-averaged and densityweighted Hamiltonian constraint violations for all the three resolutions up to t = 30 ms and show them in Figs. 2 and 3 for models MIT148 and MIT265, respectively. The initial constraint violation, which is mainly a result of the interpolation of the initial data from the surface fit coordinate employed in COCAL to the Cartesian coordinates in SACRA-MPI, is significantly reduced at the early stage of the dynamical evolution, due to the constraint violation propagation scheme. After that, the constraint violation stays at a relatively low level for the entire duration of the run.
The convergence behavior is also studied by analyzing this result. A certain quantity evaluated at one grid resolution f N could be expressed as where f exact is the exact solution if one has infinitely fine resolution, ∆ N is the grid interval associated with the resolution N , A f is a constant, and ζ is the convergence order of the code. For constraint violations, one should expect zero violation at infinite resolution. Hence, we could directly divide the constraint violation value at a given time in different resolutions and figure out ζ, given the ratio of grid intervals between two resolutions. By doing so, we have found approximately linear order convergence (ζ ≈ 1.0) for volumeaveraged constraint violation for any given time in the simulation. This agrees with our expectation: On one hand, we have already shown in our previous study that in the case of QSs, many quantities obey first-order convergence [28] for constructing initial data due to the presence of the density discontinuity at the QS surface; on the other hand, as a minimum modulus flux-limiter function (for details, cf. [67]) is used in SACRA-MPI, the order of the interpolation accuracy during reconstruction could decrease to approximately linear order if we encounter strong discontinuity to avoid unphysical oscillations: This is exactly the case on the surface of a QS. The density-weighted result shows a convergence order of ≈ 1.4 for MIT265 and 1.5 for MIT148, which further supports this understanding: The inner part of the QS, for which no discontinuity is present and thus the reconstruction is made with higher order interpolation, has larger density than the surface, and hence, the density-weighted result should show a higher convergence order. Apart from the relatively-low convergence order, we emphasize that our method for resolving the discontinuous density across the surface of the QS shows no sign of the problem. To clarify this point, we generate contour plots for the restmass density ρ on the x-y plane together with the 3-velocity field (Figs. 4 and 5). Figures illustrate that the primitive recovery scheme we introduced works well for QSs. First it is found that the density drops from the surface density of the QS (which is approximately 3.737×10 14 g cm −3 ) for two orders of magnitude within ∼ 3 grid points; this is indicated by the rapid color change from orange to green, blue, and eventually white color within 3 grid points near the surface. We note that this sharp discontinuity is still resolved by ≈ 3 grid points after 20 ms as in the beginning of the simulations for both models (see the lower panels of Figs. 4 and 5). No sign of the diffusion of this finite surface density is observed. Secondly, the velocity field, which is related to the Lorentz factor w obtained by the primitive recovery, shows a smooth configuration as found from the figures; the velocity field inside the QSs follows a rigid bulk rotation and the outside does not (since we did not put any initial velocity for the atmosphere, as found from the upper panels of Figs. 4 and 5). Especially, if we focus on the points near the surface of the QSs, it is found that the velocity on any points that are inside the surface follows the bulk motion of the QS while on the points one grid outside deviate. Since the surface of the QSs is moving in the domain (i.e., points inside/outside the QSs at one time step might move to their outside/inside at the next time step), the fact that the velocity field is well recovered for the main body of the QSs indicates the success of our implementation.

A. Triaxial rotation case
It has been suggested previously that a fast rotating compact star formed during collapse of a massive star or merger of a binary system could actually be in a configuration similar to those of Jacobi ellipsoids [68]. The evolution of such configurations as well as the GW emission properties have been investigated both analytically and numerically [33,[68][69][70][71]. According to these studies, the general understanding is that due to angular momentum loss by GW radiation, such triaxial star settles towards an axisymmetric configuration. The frequency of GWs is twice the spin frequency. In this stage, the angular velocity of the star, i.e., the GW frequency, increases as a result of the decrease of the moment of inertia during the reconfiguration of the triaxial star. When the star approaches the axisymmetric bifurcation point, other secular instabilities such as Dedekind instability and Chandrasekhar-Friedman-Schutz instability could be induced [68,70]. However, such instabilities only take place in a timescale much longer than the dynamical one as where β and β sec are the T /|W | ratio of the star and that of the onset of the secular instability, respectively. Although QS models indeed could reach higher T /|W | ratio than NSs and so is the case of the two models considered in this paper, this timescale is still as long as ∼40 s for MIT265 and ∼2000 s for MIT148, which are much longer than the timescale that can be covered by NR simulations. Hence, the simulation for secularly unstable stars is beyond the topic of this paper. On the other hand, the dynamical timescale is approximately written as which is the same order of magnitude as the spin period, 1 ms, of the models we considered and much shorter than the timescale that can be covered by NR simulations. Therefore, we can conclude that triaxial QSs are dynamically stable, if the QSs still maintain their structure after tens of rotation periods. Figures 4 and 5 indicate that the configuration of the QSs settles towards an axisymmetric one. Quantitatively speaking, after 10 ms, the length of the equatorial radius along the longer axis shrinks to 12.4 km (which is 11% shorter than the initial value of R x ) for model MIT148. On the other hand, for MIT265, R x shrinks to 10.6 km and becomes approximately identical to R y . This implies that the QSs lost nearly all the triaxial deformation in the first 10 ms. As the triaxial deformation decays during the evolution, the same occurs in the GW amplitude. As shown in Fig. 6, specifically, the GW strain amplitude decays exponentially with time. This indicates that the GW dissipation would be indeed the dominant mechanism for the reconfiguration of the triaxial QSs.
A similar numerical result was found previously for NSs [71]. The authors in [71] suggested that the timescale for reconfiguration of the stars might be incompatible with the estimation of the GW dissipation timescale that they assume: where T is the rotational kinetic energy of the star andĖ is the GW luminosity. They showed that this timescale is of the order of 1 s, and hence, it cannot account for the change in the stellar shape that proceeds in a timescale of 10 ms. However, Eq. (23) is not appropriate for estimating the relevant timescale and actually overestimates the timescale by a factor of ∼ 100. The reason for this is that the triaxial stars only need to lose the extra angular momentum (or extra rotational kinetic energy) relative to that of the bifurcation configuration, rather than all their angular momentum (or rotational kinetic energy), to settle to the axisymmetric configuration. Therefore, a better estimation should be given by where ∆J is the angular momentum difference between the initial configuration and that of the bifurcation point solution (see the illustration in Fig. 7) andJ is the angular momentum loss rate due to the GW emission which is related toĖ approximately asJ ≈Ė/Ω. According to the previous studies [28,72], triaxial stars can only obtain maximumly 10% extra angular momentum relative to that of the bifurcation point, J b . In particular, for the stars with the larger compactness, the smaller amount of the extra angular momentum can be reached (cf. Fig.8 in [28] or Table V in [72]). For the high compactness case, the extra angular momentum could be as small as ∼ 1% of J b . Therefore, it is reasonable that the triaxial stars settle to an axisymmetric configuration in a timescale of tens of milliseconds by the GW emission. In fact, according to the exponential fit shown in Fig. 6, this timescale is 17.85 ms for MIT148 and 4.35 ms for MIT265. These values are consistent with the argument that ∆J is smaller for higher-compactness QSs. The decay of the triaxial deformation results not only in the damping of the GW amplitude but also in the shift of the GW frequency, as the triaxial star spins up a little bit. We analyze the evolution of the GW frequency by performing Fourier transformation for the GW strain in different time intervals. The result is shown is Fig. 8. For model MIT148, the GW signal is approximately monochromatic with the dominant frequency at f = 2.14 kHz 5 up until 35 ms, although the power decreases by a factor of 5 in ∼ 30 ms. On the other hand, for model MIT265, the frequency of the dominant peak shifts from f = 2.64 kHz in the first time bin to f = 2.74 kHz in the second time bin. The comparison between models MIT148 and MIT265 indicates that due to angular momentum dissipation by GW radiation, the reconfiguration of the triaxial star toward axisymmetry is faster for models with higher mass, which is again a result of the fact that higher compactness stars contain much less extra angular momentum relative to that of the axisymmetric bifurcation point.
In addition, it is worth noting that, although the instantaneous GW amplitude in the beginning for model MIT265 is larger than that for model MIT148 (same occurs for the GW amplitude estimated in the initial data), it is unlikely that supramassive triaxial QSs are stronger GW sources. This is again due to the fact that the decay timescale of the GW amplitude is too short for the supramassive cases, and in the realistic detection, the signal to noise ratio could be less gained. 5 Note that the time window has a size of 10 ms for performing the Fourier transformation, the resulting frequency resolution is 100 Hz, and thus, in principle, we cannot resolve if the frequency shift is smaller than 100 Hz. . Results from three different resolutions are shown by the solid lines with different colors (blue for N = 160, red for N = 120, and green for N = 80 case). The black-dashed horizontal line indicates the GW strain estimated by the quadrupole formula according to the initial data to show consistency. We find that the norm of the GW strain (i.e., h 2 + + h 2 × ) can be well fitted with an exponential decay. We show the fit according to the best resolution results with the blue dashed curve. The best fitting exponential decay timescale is 17.85 ms for MIT148 and 4.35 ms for MIT265. The fixed frequency integration method [73] is used for obtaining the GW strain. The time shown here is the retarded time tret ∼ t − 1.5 ms as GW signal is extracted at r = 300 GM /c 2 . This fact can be directly understood from Fig. 8, as the power spectrum is indeed stronger at any given time bin for model MIT148 than for MIT265. At any rate, the signal still decays too fast to become ideal continuous GW sources even for MIT148 unless the triaxial configuration could be maintained by external angular momentum supplies (by accretion for ex- ample), known as the forced GW emission scenario [74,75].

B. Differential rotation case
It is analytically shown that an incompressible Maclaurin spheroid becomes dynamically unstable to bar formation in Newtonian gravity if T /|W | is larger than a critical value of ∼ 0.27 [33]. Such a high T /|W | ratio could be reached by rapidly spinning NSs and QSs only in the presence of differential rotation. Previous studies have shown that also in GR, the dynamical bar-mode instability could be induced for differentially rotating NSs with a sufficiently high value of T /|W |, and the critical value is found to be slightly smaller than the Newtonian value (∼ 0.25 in GR) depending weakly on the stiffness of the EoS models [34,58,76,77]. It has been shown that for such NSs, bar-like perturbation grows exponentially in the early stage until saturation is reached. Beyond the saturation, the evolution differs for models with different T /|W | values: spiral arms are formed and mass ejection subsequently occurs for relatively large values of T /|W |, while for models with T /|W | close to the critical value, no significant sign of the spiral arm structure and mass ejection is found [34]. In this paper, both MIT275dr and APR206dr possess T /|W | larger than 0.26. Thus, we expect that the spiral arm structure is formed. These models could be also useful for un- derstanding the similarity or difference between the bar-mode instability of QSs and NSs. The l = m = 2 mode GW strain, which is directly related to m = 2 bar-mode, during the evolution is shown in Fig. 9. It is worth noting that we do not impose any bar-like perturbation in the beginning of the evolution, which is different from the treatments in previous studies [34,58] as our main focus is not to explore the parameter space for un-/stable differentially rotating QSs, but to understand the capability of the code as well as the difference between QSs and NSs. In our evolution, bar-mode instability spontaneously sets in from a random perturbation that should present in any numerical simulation, and exponentially grows from the tiny perturbation spending ∼ 10-20 ms for both models. Saturation is achieved approximately at t ret = 19.75 ms for model MIT275dr and at t ret = 12.25 ms for model APR206dr, and then, the GW strain decreases to a relatively lower level.
The snapshot of density contours on the equatorial plane straightforwardly demonstrates the growth of the bar mode (see Figs. 10 and 11). The evolution of the density contours for model APR206dr is quite consistent with the previous results of differentially rotating compact stars with a relatively large value of T /|W | [34,58]: the star is significantly distorted when the bar-mode growth saturates (the middle panel, which we choose the time when the peak of the GW strain is reached) and the spiral arm structure forms afterwards (the lower panel). On the contrary, the result of model MIT275dr is similar to the previous results with T /|W | slightly larger than the critical value: the star adjusts to ellipsoidal structure without formation of spiral arms.
In spite of all the differences mentioned above, it is found that the growth rate of the bar-mode instability is quite sim- is achieved during the simulation (cf. Fig. 9) and at tret = 26.50 ms, the GW amplitude decreases to a relatively low level.  ilar for the two models considered in this work if we normalize the simulation time by the central rotation frequency (cf. Fig. 12). For both models, the GW amplitude grows for about 3 orders of magnitude in the duration of Ω c t = 100-200. More careful analysis yields an e-folding timescale of about 2.24 central rotation periods for model APR206dr and 2.04 for model MIT275dr. This result is also consistent with the range derived in previous studies [34,58], indicating that the instabilities found in this work are indeed the dynamical bar-mode one.
The mass ejection is also studied (as shown in Fig. 13) by applying the criterion that matter becomes unbound if u t is smaller than −1.0. For model APR206dr, a burst of mass ejection is found approximately at the same time when the GW strain reaches the maximum value and in total about 6 × 10 −4 M is ejected in the end. However, the situation is a bit more complicated for the QS case due to the artificial mass ejection in the beginning. The main reason for this artificial mass ejection is the existence of matter with density smaller than the surface density near the surface when initial data is interpolated to the grid points for dynamical evolution. Such matter has to be treated as ideal gas once the dynamical evolution started and this may drive part of those matter away from equilibrium. Nevertheless, a rise in the ejected mass is still found at t ∼ 20 ms, which corresponds to the duration of the growth in the GW strain. A total mass ejection of about 3.5 × 10 −4 M is found after subtracting the mass ejection before t = 18 ms (i.e., before the bar-mode begins to exponentially grow). A more careful study on this artificial ejec- tion is shown in Appendix B. It is found that, for the current implementation, the quantitative measurement of the ejecta is reliable for slowly rotating models but uncertainty will significantly increase for models with extremely large angular velocity (as is the case for MIT275dr) due to artificial ejection.
An interesting feature is that although model MIT275dr possesses a higher T /|W | value than model APR206dr (cf . Table II), the spiral arm formation and mass ejection are much less significant than the NS case, showing that the critical T /|W | value is likely to be higher for the QS model we considered here. There are several possible reasons for this. First of all, for QSs the density profile is much more uniform than for NSs, and hence, the T /|W | ratio is possible to be higher to approach the limit for incompressible stars. Secondly, previous studies concentrate in a differential rotation degree aroundÂ = 1.0 for NSs and found weak correlation between critical T /|W | andÂ parameter. However, for MIT275dr model considered here, a quite largeÂ parameter (3.0) is adopted and for NSs with such small differential rotation degree it is not possible to find dynamically unstable solutions. The parameter space for dynamical bar-mode instability for such largeÂ parameter might be different and we need to perform more systematic studies for differentially rotating QSs with different parameters to fully understand this in the future.

V. DISCUSSION AND CONCLUSION
In this paper, we introduced our approach to perform NR simulations for QSs which have a finite surface density. We made two major efforts for successfully implementing QSs described by the MIT bag model into our NR code SACRA-MPI. First, we appropriately chose the EoS parameters in such a way to make specific enthalpy continuous and monotonic across the surface of the QS, without changing any quantities of the QSs described by this model. Secondly, we developed a procedure for the primitive recovery suitable for QSs. Specifically, in our approach, the primitive recovery is performed separately for matter inside and outside QSs, and in addition, we introduced an analytical approach for the QS part. We then examined the implementation by performing simulations for two triaxial QS models with different mass and a differential rotation model which is dynamically unstable to the bar-mode deformation. The surface of the QSs with the finite density is found to be well resolved for the entire time during the evolution and the primitive recovery is correctly executed for matter both inside and outside the QSs. Convergence studies revealed a first-order convergence as in the case of initial data construction for QSs [28]. The success in capturing the evolving surface of QSs and mass ejection from it is also the important gain for us to move to the simulations for BQSs. Besides, this method can in principle be extended to perform simulations on compact stars with a first-order strong-interaction phase transition (for which there is a density jump inside the star). We plan to explore this topic in the future.
For the triaxially rotating QSs, the results are found to be qualitatively similar to the case of NSs [71]. The triaxial deformation of QSs decays as it radiates angular momentum by the GW emission in a timescale of 10 ms, which is found to be consistent with the GW dissipation timescale. As a result, the amplitude of GWs decays exponentially. It is also found that for higher compactness, the QS settles to an axisymmetric configuration faster, due to smaller amount of the extra angular momentum contained in comparison with that of the bifurcation point solution. We found that the GW emission of the model with typical mass of an individual compact star is approximately monochromatic, while the GW emission of the supramassive case shifts to a higher frequency within ∼ 10 ms as the star resumes axisymmetry at this time and the GW amplitude is significantly reduced. As a conclusion, in spite of the fact that rotating QSs could reach a T /|W | ratio higher than that of NSs, and thus, supramassive triaxially rotating QSs do exist, it is unlikely that such supramassive triaxial rotating QSs could be a better source for continuous GWs, as the signal from such objects does not last for more than ∼ 10 ms.
For the case of differentially rotating stars, we confirmed that QSs with much smaller differential rotation degree could still experience the bar-mode instability as long as the T /|W | ratios is large enough (i.e., 0.27), as expected in our previous initial data studies [29]. In particular, even without initially artificial bar-like perturbation, the bar-mode deformation of QSs still exponentially grows in a comparable timescale of the NS cases. The result is consistent with the previous results of NSs with T /|W | ratios not much larger than the critical value [34,58]. This might be a result of the difference between the structure of QS and NSs or the fact that theÂ parameter considered here for QSs is quite different from previous studies of NSs. To fully understand the parameter space of bar-mode instability of differentially rotating QSs, a more systematic investigation needs to be done in the future. In order to more quantitatively understand the validity of our numerical implementation, it is helpful to investigate models which are simple but the evolution of them could be predictable. Evolution of TOV solutions, i.e., non-rotational solutions which should be static, is one good choice to understand the accuracy and convergence behavior of our methods quantitatively.

ACKNOWLEDGMENTS
For this purpose, we have prepared evolution of a TOV solution with M ADM = 2.07 M , which is quite close to the TOV maximum mass of our EoS model, in 3 different resolutions. As introduced in Sec. III C, nine FMR domains are used for the simulation and the half size of the finest level being 8GM /c 2 ≈ 11.8 km for this test run. There different resolutions with ∆ x = ∆ y = ∆ z = 148, 98.4, and 73.8 m are chosen, which corresponds to N = 80, 120 and 160 grid points along one direction.
Due to numerical perturbation which is present in any code, the central density will oscillate in dynamical timescale even for TOV solutions of NSs. Since we are considering the TOV solution whose mass is very close to its maximum, the QS will become unstable and collapse to BH if the numerical method is not reliable. Hence, it is important to check the convergence behavior on the oscillation of the central density for this test. As can be seen in Fig. 14, the oscillation amplitude of the central density decreases with time for all three resolutions. The amplitude is ∼ 3% in the beginning and then stabilize to roughly 1% after 2 ms. Particularly, the oscillation amplitude is smaller in higher resolution and can be used to estimate the convergence order.
in which f 160 , f 120 and f 80 here represent the oscillation amplitude in three different resolutions, respectively. Since the oscillation phase could be different in 3 configurations and the amplitude could vary for every period, we have chosen roughly ten oscillation cycles from t = 3.5 ms to 5 ms when the oscillation is already more stabilized for all 3 resolutions and define f N as the difference between maximum and minimum value of the central density in this duration for the corresponding resolution. With such a definition, R is found to be 0.2516 which indicates that second order convergence is achieved for this quantity (ζ ∼ 2). This higher order convergence supports our explanation for the relatively lower convergence order found in the evolution of constraint violation (cf. the discussion in Sec. III C): indeed, the local hydrodynamic quantity follows higher order convergence if it is not close to the surface. Apart from local quantities, we have also investigated global density profile in radial direction, which could quantitatively demonstrate the capability of the implementation in resolving the discontinuous density across the surface of a QS. The comparison of the initial and final (at t ∼ 6.7 ms) density profile in the N = 160 resolution case is shown in Fig. 15. As can be seen, the density profile inside the QS is approximately unchanged during the entire evolution. The most obvious difference is the density discontinuity across the surface as can be seen from the inset in Fig. 15: the very steep density drop at the surface becomes smoother during the evolution. Nevertheless, the density drop from ρ s to 0 is resolved by ∼3 grid points in the end, which is similar to what can be seen from the rotating solution cases.
Besides local quantities as mentioned above, we have also examined the performance of our method by looking at global quantities such as the conservation of baryonic mass and mass that becomes unbound during the simulation. Throughout the simulation, the change in baryonic mass stays as low as ∼ 10 −8 M . The unbound matter is at a negligible amount of ∼ 10 −12 M even at its maximum, which is as expected for a static solution.

Appendix B: On the Artificial Mass Ejection
In Sec. IV B, it is mentioned that unphysical mass ejection of ∼ 10 −4 M is found in the very beginning of the simulation for model MIT275dr. A possible explanation for this artificial mass ejection is that when the initial data is interpolated to the SACRA-MPI grids for evolution, it is unavoidable that some part of the grids close to the surface of the QS might obtain a density which is finite but smaller than ρ s , as it is not possible for the Cartesian coordinates used in SACRA-MPI to cover exactly the surface of the QS. And hence this part of the matter will be treated as atmosphere and might be driven off equilibrium in the later evolution. Given certain amount of initial energy perturbation due to the reason above, it should always be easier to eject matter from the surface of a QS if the angular velocity is larger as the binding energy for those matter are smaller.
To verify this explanation, and more importantly to know the reliability of the current implementation in estimating the amount of ejecta for binary merger case in the future, we have prepared several uniformly rotating QS models with different angular velocity and investigated the mass ejection of those models in the early stage of the evolution. We . In addition, in order to verify that it is the angular velocity rather than mass that determines the final amount of the spurious mass ejection, we have constructed one additional model with smaller central density (ρ c = 6.304 × 10 14 g cm −3 ) and with mass and angular velocity of (M = 2.25 M , Ω = 7402.0 rad s −1 ). Furthermore, in order to estimate the maximum possible spurious mass ejection, we have constructed a model with mass and angular velocity very close to the mass shedding limit (M = 2.95 M , Ω = 8682.4 rad s −1 ) 6 . The result is shown in Fig. 16. The artificial mass ejection monotonically increases with the angular velocity. Particularly, the comparison between the dashed green curve and the other curves clearly show that the amount of artificial ejection is affected mostly by the angular velocity of the model instead of the mass. The dependence still holds even if we consider the TOV solution case as reported in Appdendix A and the MIT275dr model. As a consequence, the systematic 6 The uniform-rotation mass shedding limit of this EoS model is found at M = 3.02 M and Ω = 9149.6 rad s −1 . uncertainty in the quantitative measurement of physical ejection will be higher for models with larger angular velocity. A maximum uncertainty will be order of 10 −4 M as indicated by the results of the model M = 2.95 M . In addition, since the initial artificial mass ejection is related to the surface part of the star where both density and pressure gradient are discontinuous, the convergence behavior is expected to be worse than other part of the star. We have done the test on the model MIT275dr in three different resolutions and the result is shown in Fig. 17. As can be seen, the influence of the resolution on the amount of ejecta is very weak. Especially, the result is approximately the same for N = 120 and N = 160 cases. Further increasing the resolution within the current affordability of computational resources does not seem to be a viable way to resolve the problem.
For the purpose of extracting the amount of ejecta for binary quark star mergers, the relevant angular velocity should be order of ∼ 2000 rad s −1 for a typical 1.4 M -1.4 M binary at a separation of 40 km, for which the initial artificial mass ejection should be much less than the MIT275dr model considered in this paper. Our preliminary BQS merger simulation confirmed this and we will cover more details about the binary case in the future.