Fully general-relativistic simulations of isolated and binary strange quark stars
FFully general-relativistic simulations of isolated and binary strange quark stars
Zhenyu Zhu
1, 2 and Luciano Rezzolla
1, 3, 4 Institut f¨ur Theoretische Physik, Max-von-Laue-Strasse 1, 60438 Frankfurt, Germany Department of Astronomy, Xiamen University, Xiamen 361005, China Frankfurt Institute for Advanced Studies, Ruth-Moufang-Strasse 1, 60438 Frankfurt, Germany School of Mathematics, Trinity College, Dublin 2, Ireland (Dated: February 16, 2021)The hypothesis that strange quark matter is the true ground state of matter has been investigated for almostfour decades, but only a few works have explored the dynamics of binary systems of quark stars. This is partlydue to the numerical challenges that need to be faced when modelling the large discontinuities at the surfaceof these stars. We here present a novel technique in which the EOS of a quark star is suitably rescaled toproduce a smooth change of the specific enthalpy across a very thin crust. The introduction of the crust has beencarefully tested by considering the oscillation properties of isolated quark stars, showing that the response ofthe simulated quark stars matches accurately the perturbative predictions. Using this technique, we have carriedout the first fully general-relativistic simulations of the merger of quark-star binaries finding several importantdifferences between quark-star binaries and hadronic-star binaries with the same mass and comparable tidaldeformability. In particular, we find that dynamical mass loss is significantly suppressed in quark-star binaries.In addition, quark-star binaries have merger and post-merger frequencies that obey the same quasi-universalrelations derived from hadron stars if expressed in terms of the tidal deformability, but not when expressed interms of the average stellar compactness. Hence, it may be difficult to distinguish the two classes of stars ifno information on the stellar radius is available. Finally, differences are found in the distributions in velocityand entropy of the ejected matter, for which quark-stars have much smaller tails. Whether these differences inthe ejected matter will leave an imprint in the electromagnetic counterpart and nucleosynthetic yields remainsunclear, calling for the construction of an accurate model for the evaporation of the ejected quarks into nucleons.
I. INTRODUCTION
The detection of gravitational waves (GWs) from the coa-lescence of two compact objects with masses that can be as-sociated to compact stellar objects has been reported recentlyby the LIGO-Virgo Scientific Collaboration [1–3]. One ofthese GW detections, GW170817, was accompanied by anelectromagnetic counterpart [4, 5] and has therefore been as-sociated with the merger of a binary system of neutron starswhich have long been proposed to be behind the production ofshort gamma-ray bursts [6–9]. Furthermore, the detection of akilonova emission through the ejected material and the subse-quent r-process nucleosynthesis [10–18] has provided furtherevidence that the GW signal in GW170817 must have beenproduced by the merger of two compact matter objects.Combining the knowledge of nuclear physics, general rela-tivity, and numerical relativity simulations, this detection hascertainly helped to deepen our understanding of cold densematter equation of state (EOS) through a tight constraintson maximum mass, radii and tidal deformability of the neu-tron stars [19–31]. Besides the more conventional scenarioof the merger of purely hadronic compact stars leading to apurely hadronic merged object, other possibilities have beenexplored in great detail. One of them is the possibility thatthe merging objects were hybrid (or twin) stars [32–38], thata phase transition to quark matter could have taken place afterthe merger [39–41], or that the merger involved strange-quarkstars [42, 43]. All of these different scenarios are in princi-ple compatible with the GW signal of GW170817, which was necessarily limited to the inspiral only.The solution of the equations of general-relativistic hy-drodynamic (GRHD) or general-relativistic magnetohydro-dynamic (GRMHD) are indispensable tools for the accuratemodelling of these scenarios and the quantitative prediction ofthe signals produced by the inspiral and merger, be it throughthe gravitational radiation, the emission of neutrinos, or theejection of matter.Over the years, the scenario of the merger of fully hadroniccompact stars has been studied by GRHD and GRMHD sim-ulations in great detail [8, 13, 15, 18, 44–55]. Among the nu-merous results that have been obtained with these simulations(see [56, 57] for some reviews), two are particularly relevantfor the results presented in this work. The first one is aboutthe spectral properties of the GW signal – both during the in-spiral and after the merger – have been analysed in great detailand shown to follow quasi-universal relations in terms of thestellar tidal deformability or compactness [58–65]. The sec-ond one is instead related to the matter ejected at merger andafter the merger [13, 15, 51, 54, 55, 66–71], and the impactthat it has on r-process nucleosynthesis, on the lifetime of themerged object [72], and on the maximum mass of compactstars [19, 21, 31].While the fully hadronic scenario has been covered in greatdetail, the alternative scenarios involving matter in differentstates – either before or after the merger – has so far beenconsidered less extensively. The few investigations performedso far have in fact been limited to the analysis of the post-merger GW signal when a phase transition sets in after themerger and leads to clear signatures in the GW signal [39–41].In particular, these studies have highlighted that differences,sometimes significant, can be found in the post-merge GW a r X i v : . [ a s t r o - ph . H E ] F e b spectrum and that the universal relations found in the case ofhadronic stars are obviously broken if a considerable quarkcore is produced.The literature is even scarcer when it comes to simulationsexploring the inspiral and merger of quark-star binaries thatcomposed of pure strange quark matter (SQM). Indeed, thefirst and only works so far are more than a decade old andhave been obtained using smooth particle hydrodynamics anda conformally flat approximation to general relativity [73, 74].This is in great part due to the considerable additional difficul-ties that the numerical simulation of these objects implies andthat originates from the very sharp decrease in density and en-thalpy at the surface of the quark star.We recall that the SQM hypothesis was firstly proposed byWitten and suggesting that SQM, rather than nuclear matter,is the absolute ground state of matter [75]. In this hypothesis,SQM is mainly composed by up, down, and strange quarks,with a small fraction of electrons also present. Because ofthe self-bound properties of SQM, objects composed of thismatter can exist in any size, from scales as small as those ofnuclei, to scales as large as that of a compact star. Indeed,in this hypothesis, a quark star composed of SQM has beenproposed as a possible candidate of compact stars. The prop-erties of single quark stars have been studied by numerousworks [76–88], while the study of the binary quark-star merg-ers have been explored far less [43, 73, 74, 89–91].The detection of the kilonova signal AT2017gfo associatedwith the GW170817 event has been considered by some as astrong evidence against the existence of SQM. This is becausethe ejected material of a binary quark-star merger would becomposed by SQM that – when assumed to be the most sta-ble form of matter – cannot be an efficient source of r -processnucleosynthesis. At the same time, recent studies have high-lighted that the SQM scenario can be still conciliated with thesignal from AT2017gfo if the SQM can evaporate into nu-cleons as a result of the high temperatures reached after themerger. In this case, most of ejected SQM from the quark-starbinary would have evaporated into nucleons and could havetherefore contributed to the kilonova signal in AT2017gfo[90, 92, 93].With the goal of studying the scenario of the merger ofquark stars – and thus explore the possibility of SQM evap-oration in far greater detail – we have carried out the first fullygeneral-relativistic simulations of the inspiral and merger ofquark stars. Contrasting their evolution with a system of com-pact stars that have very similar properties in mass and com-pactness but are fully hadronic, we have been able to isolatethree important features of the merger of quarks stars. First,their GW spectral properties are in agreement with the quasi-universal behaviour found for hadronic stars during the inspi-ral, but differs in the post-merger phase. Second, becauseof the intrinsic self-boundness of these objects, the amountof ejected mass is considerably smaller, with binary quarkstars ejecting up to 30 times less mass than a correspondinghadronic binary. Finally, as natural to be expected for matterthat is colder and more self-bound, the ejected matter containsconsiderably smaller tails in the corresponding distributions ofvelocity and entropy. The structure of the paper is as follows. In Sec. II we dis-cuss in detail the mathematical and physical setup that wasnecessary to develop in order to carry out the numerical sim-ulations. Section Sec. III is instead dedicated to a carefultest of the validity of our approach when considering the os-cillation properties of isolated quark and hadronic stars, andtheir match with perturbative studies. Section IV is used topresent in detail the results of our simulations, including boththe overall dynamics of the merger and the outcomes in termsof GW signal and ejected matter. Finally, our conclusions andprospects for future work are presented in Sec. V. II. MATHEMATICAL AND NUMERICAL SETUPA. Equations of state
The EOS of the SQM employed in our simulations was cho-sen to be the MIT2cfl EOS [42, 94]. This EOS makes useof the MIT bag model with additional perturbative QCD cor-rections [42, 95–97], and satisfies the constraints of having amaximum mass above two solar masses as required by the ob-servations [98, 99], and a tidal deformability compatible withthe constraints from GW170817. In this EOS a colour-flavor-locked (CFL) phase is assumed to be present and consequentlythe number densities of all the flavor of quarks [up (u), down(d), and strange (s) quarks] are the same, i.e., n u = n b = n s . (1)As a result, the baryon number density defined as n B := 13 ( n u + n b + n s ) , (2)can be directly converted to the rest-mass density ρ that en-ters in the hydrodynamic equations and is evolved numeri-cally, as ρ := m B n B , where the average baryonic mass is m B := ( m u + m b + m s ) . The quark rest-masses assumedhere are m u = m d = 0 , and m s = 100 MeV , so that m B = 33 .
33 MeV . However, this conversion is not pos-sible if the CFL phase is not present. Furthermore, if therest-mass m B depends on the fractions of the quark flavorspresent – and these vary in space and time – then m B becomesa density-dependent variables and the rest-mass conservationlaw ∇ µ ( ρu µ ) = 0 , where u µ is the fluid four-velocity, isno longer valid for SQM. In this case, however, the baryonicnumber density n B will always be conserved and could playthe role of rest-mass density ρ in the GRHD equations.Possibly the most serious challenge in modelling compactstars made of SQM is that – because of their self-bound prop-erty – the surface of such objects is characterized by a sharptransition at the stellar surface, where the pressure goes to zeroat a nonzero rest-mass density. As a result, a large, intrin-sic density jump is present at the stellar surface; by contrast,hadronic stars have surfaces near which the density rapidly de-creases, but that goes to zero when the pressure is zero. Whenconsidering static solutions, as for example when construct-ing stellar models of isolated quark stars, such a density jump p ( d y n e / c m ) pp eff ρ/ρ s h crust core hχh eff − − − − p ( d y n e / c m ) pp eff − − − ρ/ρ s h crust core hχh eff FIG. 1.
Left panels: behaviour of the pressure p (top part) and of the specific enthalpy h (bottom part) of the MIT2cfl EOS shown as a functionof the rest-mass density normalised to the nominal (i.e., without a crust) value at the surface ρ/ρ s . The values presented refer to the originalpressure (orange solid lines) and to the rescaled one [cf. Eq. (4)] (green dashed lines). Note the appearance of a small jump in h for verysmall densities. Right panels: same in the right panels, but when the rescaled rest-mass density is shown in a logarithmic scale to highlight thesmoother transition. can be handled analytically by matching the stellar interiorwith the exterior vacuum [42, 94, 100, 101]. However, in thecontext of a hydrodynamic simulation, such a discontinuityrepresents the exemplary condition for the development of astrong shock that would lead to an artificial oscillation in thebest-case scenario or to a numerical failure in the most realis-tic case. Clearly, a treatment aimed at smoothing this strongdiscontinuity into a region with small but finite size is neces-sary for a numerical evolution.A naive solution at the level of the EOS may consist inthe introduction of a polytropic piece in the pressure depen-dence, i.e., p = kρ Γ , from the rest-mass density, thus effec-tively introducing a thin but nonzero “crust” in the quark star.The presence of a crust in a quark star and its implicationshave been studied in detail in a number of works [102–104].In practice, however, this polytropic treatment of the crust isnot sufficient as the extrapolated polytropic EOS cannot effec-tively decreases the specific enthalpy h around the surface. Inother words, given the small width of the crust and the physi-cal constraints on the polytropic exponent, a single-polytropeprescription is not able to reduce the specific enthalpy fromthe value it has at the base of the crust, i.e., h ∼ , tothe value it should attain at the surface of the quark star,i.e., h = 1 . As a result, a jump in the specific enthalpy ap-pears inevitable in this case. To appreciate the consequencesof this jump, it is useful to recall the jump condition across acontact discontinuity [105] (cid:74) p (cid:75) = (cid:74) h (cid:75) = 0 , (3)where (cid:74) φ (cid:75) measures the jump of the scalar function φ acrossthe contact discontinuity. Clearly, if (cid:74) h (cid:75) (cid:54) = 0 , a pressure jumpwill develop, which will blow away the surface layers of the quark star into the atmosphere . Because it is not possible toremove this jump, the quark star will rapidly diffuse away intothe atmosphere.Figure 1 reports with a solid orange line the behaviour ofthe pressure (top panels) and of the specific enthalpy (lowerpanels) as a function of the rest-mass density normalized tothe putative rest-mass density at the surface of the star (i.e., be-fore the addition of a crust), ρ s . The left panels, in particular,adopt a linear scale in ρ/ρ s , while the right panels a loga-rithmic scale to highlight different properties of the EOS. Inparticular, it is very apparent in the lower left panel that ajump in h is present near ρ/ρ s (cid:39) (see bottom-left panel).Also reported with a solid orange line in the upper panels ofFig. 1 is the behaviour of the pressure, which in the crust wemodel with a polytrope with polytropic constant k = 0 . andpolytropic index Γ = 1 . (see top-right panel). Figure 1 alsoallows one to appreciate that the specific enthalpy is very largeacross the crust, but also essentially constant, thus suggestinganother and simple solution to the problem discussed above.In practice, defining as “core” whatever is the genuine in-terior of the quark star, we operate within the core a rescalingof the rest-mass density, of the specific enthalpy, and of thespecific internal energy (cid:15) of the type ρ → ρ eff , h → h eff , (cid:15) → (cid:15) eff , (4)where ρ eff := χρ = χm B n B , (5) We recall that numerical-relativity codes require the presence of a verylow-density atmosphere outside of the compact stars for the solution of theequations of relativistic hydrodynamics (see [105] for a discussion). In oursetup, the atmosphere is set to be at a threshold density that is orders ofmagnitude below the maximum density in the star. h eff := hχ , (6) (cid:15) eff := 1 + (cid:15)χ − , (7)while the energy density and pressure are left unchanged,i.e., e eff = e , p eff = p . In essence, we are effectively intro-ducing a rescaled (larger) effective baryon mass m eff B := χm B in terms of the arbitrary constant χ (cid:29) , such that the specificenthalpy is reduced, hence reaching the goal of decreasing thedifferences between the stellar interior and the fluid in the at-mosphere. Indeed, the rescaled specific enthalpy still has ajump for ρ/ρ s (cid:39) , but this is much smaller and does not leadto a diffusion of the star.The treatment of the EOS within the crust , on the otherhand, remains the same as before, namely, with a single poly-trope joining from the stellar surface down to the values of thedensity associated with the atmosphere. In this case, the poly-tropic constant and exponent can be much larger and we haveused k = 8 . and polytropic index Γ = 1 . (see top-rightpanel).A few important aspects of our novel approach need to beremarked at this point. First, the rescaling in Eqs. (5)–(7)described above is done only at the level of the numericalsolution of the GRHD equations and to remove the difficul-ties introduced by the sharp jump at the stellar surface. How-ever, after the completion of each timestep, all the physicalquantities entering, for instance, on the right-hand-side of theEinstein equations, are evaluated with the physically correctvalue of m B and hence yield the physically contributions tothe energy-momentum tensor.Second, although artificial, the introduction of the crustdoes not provide a perceptible variation to the global prop-erties. While we will demonstrate this in Sec. III, where wewill compare the oscillation properties of an isolated star withthe perturbative expectations, it suffices to say here that therest-mass in the crust is minute, i.e., ∼ × − M (cid:12) , as isits spatial extension, which is restricted to two grid cells andtherefore has a width of (cid:39)
240 m .Third, while χ is essentially arbitrary, some values workbetter than others. A bit of experimenting has revealed thatfor our setup a value of χ = 25 . yields very accurate results(see discussion in Sec. III), thus corresponding to an effectivebaryon rest-mass of m eff B = 850 MeV .Finally, the handling of the MIT2cfl EOS for the actualmerger requires the addition of a thermal part to the EOS. Wedo this in close analogy with what done routinely for simula-tions of hadronic stars described by cold EOSs. In essence, weaccount for the additional shock heating during the merger andpost-merger phases by including thermal effects via a “hybridEOS”, that is, by adding an ideal-fluid thermal component tothe cold (underscore “c” below) EOS [105] p = p c + p th , (8) p th = ρ(cid:15) th (Γ th − , (9) (cid:15) th = (cid:15) − (cid:15) c ( ρ ) , (10)where Γ th is the thermal index. Combining Eqs. (7) with thethermal part of EOS, one can easily derive the effective quan- tities of cold part p eff , c = p c , (11) h eff , c = h c χ , (12) (cid:15) eff , c = 1 + (cid:15) c χ − . (13)Since we find it important to contrast the dynamical be-haviour of merging quark stars with the corresponding oneof hadronic stars with the same total mass, we have also con-sidered the merger of a binary system subject to a hadronicEOS for comparison. Our choice has fallen on the DD2 EOS[106], which is compatible with the present observational con-straints, and whose tidal deformability for a M = 1 . M (cid:12) star is close to that of a quark star of the same mass. Atthe same time, the radius of the quark star is significantlysmaller: R = 11 .
81 km for the MIT2cfl quark star versus R = 13 .
21 km for DD2 hadronic star.In the framework of a comparative assessment of the dy-namics of SQM and hadronic-matter binaries, and despitethe fact that the hadronic DD2 EOS is a finite-temperatureEOS employed in many GRHD or GRMHD simulations [13,15, 107–109], we have adopted a hybrid-EOS approach[cf. Eqs. (8)–(10)] also for the DD2 EOS, of which we haveretained only the zero-temperature slice. An immediate dis-advantage of this approach is that the consistent knowledgeof some thermodynamical quantities, such the temperature orthe entropy, are missing and need to be estimated in alterna-tive manners. In this case, additional assumptions – someof which are not necessarily realistic – are required to ex-tract these quantities, at least to some degree. However, sincethe same approximations are made for both EOSs, the expec-tation is that the systematic differences we can find in thisway will persist also when considering more advanced andtemperature-dependent EOSs for SQM.More specifically, in the case of ideal-fluid EOS, the tem-perature T is proportional to the average kinetic energy perparticle and we can therefore express the specific internal en-ergy (cid:15) as (cid:15) = k t Tm B , (14)where k t is a constant and m B is mass per baryon. We extendthis expression to our EOSs by rewriting it as (cid:15) = k t ( T − T c ) m B + (cid:15) c ( ρ ) , (15)where T c is the temperature of the cold part of the EOS andwhich we take to be to T c = 0 .
01 MeV . Using now Eq. (15)and recalling that for transformations at constant density d(cid:15) = T ds , we can compute the specific entropy as s = k t m B log (cid:18) (cid:15) − (cid:15) c k t T c /m B + 1 (cid:19) . (16)Finally, specifying k t = 20 for both EOSs, we can computethe entropy per baryon s B = m B s once (cid:15) and (cid:15) c are known. B. Equations of relativistic hydrodynamics
The rescaling procedure described in the previous Sectionhas a simple and direct impact on the equations of relativistichydrodynamics that we solve numerically. Adopting a flux-conservative formulation of the GRHD equations, the conser-vation of rest-mass, energy, and momentum can be cast intoan “effective” Valencia formulation ∂ t ( √ γ U eff ) = ∂ i ( √ γ F i eff ) = S eff , (17)where we have employed the rescaling in Eqs. (5)–(7) U eff := DS j E = ρ eff Wρ eff h eff W v j ρ eff h eff W − p eff , (18) F i eff := αv i D − β i DαS ij − β i S j αS i − Eβ i , (19) S eff := αS ik ∂ j γ ik / S i ∂ j β i − E∂ j ααS ij K ij − S j ∂ j α , (20)and where S ij := ρ eff h eff W v i v j + p eff γ ij . Note that thequantities γ ij , γ , α , β i , and K ij are, respectively, the spa-tial three-metric, its determinant, the lapse function, the shiftvector, and the extrinsic curvature.Note that γ ij , α , β i , and K ij are related to the spacetimemetric, which is computed from the Einstein equations via theenergy-momentum tensor that is always computed from un-rescaled quantities; hence, these spacetime variables are un-affected by the rescaling operated in the EOS. Similarly, thequantities S i , E , and S ij all depend on the product ρ eff h eff ,where the scaling cancels out (i.e., ρ eff h eff = ρh ), or on thepressure that is not the rescaled in the core (i.e., p eff = p );hence they are not modified by the rescaling. Similarly,the three-velocity v i and the corresponding Lorentz factor W := (1 − v i v i ) − / are unaffected by the rescaling in theEOS. Consequently, the whole set of GRHD equations (17)–(20) remains unchanged in the core under the rescaling. Theonly exception is represented by the first equation of the set,i.e., the one involving the conservation of rest-mass. However,since the only change is in the effective baryon mass – whichappears in the conservation equation just as a multiplicativeconstant – the change is practically trivial. C. Numerical setup: Initial data
The initial data for binary stars was generated makinguse of the publicly available
Lorene code [110], whichis a multi-domain spectral-method code computing quasi-equilibrium irrotational binary configurations of compactstars. Using this code, and as discussed above, we have com-puted initial binary configurations of both quark stars and of
TABLE I. Properties of the quark and hadron stars considered here,and distinguished in whether they refer to the isolated configurationsor to binaries. Reported are the mass M , the baryon mass M b , theradii R , and the tidal deformability Λ . In addition, in the case of thebinaries, we also report the relevant frequencies of the gravitational-wave signal ( f max , f , and f ), and the ejected matter M ej .Quark star Hadronic star EOS type MIT2cfl DD2 M [ M (cid:12) ] (isolated) 1 .
40 1 . M b [ M (cid:12) ] (isolated) 0 .
062 1 . R [km] (isolated) 11 .
90 13 .
22Λ (isolated) 657 .
78 698 . M [ M (cid:12) ] (binary) 1 .
35 1 . M b [ M (cid:12) ] (binary) 0 .
059 1 . R [km] (binary) 11 .
81 13 .
21Λ (binary) 789 .
32 857 . f max [kHz] (binary) 1 .
684 1 . f [kHz] (binary) 1 .
919 1 . f [kHz] (binary) 2 .
399 2 . M ej [10 − M (cid:12) ] (binary) 0 .
10 2 . hadronic stars. Independently of whether we have consideredthe MIT2cfl or the
DD2
EOS, the properties of the binaryhave been set to be the same: the initial separation in the bi-nary was set to be
45 km and the two stars have the same massof M = 1 . M (cid:12) .It is useful to remark that the introduction of a crust forthe quark star was important also for the calculation of theinitial data. Furthermore, the computation of the solution ofbinary quark stars with Lorene required extra care. In partic-ular, because of the steep drop in rest-mass density across thecrust of the quark star and of the presence of discontinuities inhigher-order derivatives at the crust-core interface (we recallthat the EOS is continuous but with discontinuous derivativesat crust-core interface), a single interior coordinate domainemployed to cover the quark star turned out to be insufficientand prevented the convergence to an accurate solution. For-tunately, however, the addition of a second domain at higherresolution to cover the crust was sufficient to yield the conver-gence to an accurate solution.
D. Numerical setup: evolution equations
The simulations presented here have been performed withthe publicly available GRHD code
WhiskyTHC [111–113],which is a high-order, fully general-relativistic code for thesolution of the equations of relativistic hydrodynamics andcompatible with the
Einstein Toolkit [114]. The hy-drodynamic equations were solved by using a high-resolutionshock-capture (HRSC) [105] approach having with LocalLax-Friedrichs (LLF) flux-split and Monotonicity Preserving t [ms] − . − . . . . ( ρ c ( t ) − ρ c ( )) / ρ c ( ) Quark StarQuark StarQuark Starlowmediumhigh t [ms] − . − . . . . ( ρ c ( t ) − ρ c ( )) / ρ c ( ) Hadronic StarHadronic StarHadronic Starlowmediumhigh
FIG. 2. Evolutions of the normalised variations in the stellar central rest-mass density ρ c ( t ) for either an oscillating quark star (left panel) orfor a hadronic star (right panel). Curves of different shading refer to different resolutions, with the high-resolution data being shown with thedarkest shade. McLachlan [119].In order to cover a large enough spatial domain and hencecompute accurately the information on the ejected matterwhile keeping a high resolution on the compact stars, anadaptive-mesh refinement (AMR) approach was employedand handled by the
CARPET driver [120]. The number of re-finement levels depends on the problem considered and wasof levels in the case of isolated stars and of levels inthe case of merging binaries. In this way, the correspond-ing total extents of the computational domain were set to M (cid:12) (36 km) and M (cid:12) (3051 km) , respectively. Fur-thermore, each of the scenario considered – isolated starsor binary mergers – has been evolved with simulations atdifferent resolutions. More specifically, for the finest gridwe have used three resolutions in the case of isolated stars,i.e., dx = 0 . M (cid:12) (236 m) , dx = 0 . M (cid:12) (177 m) ,and dx = 0 . M (cid:12) (118 m) , which we will indicate as“low”, “medium”, and “high” resolution in following. Sim-ilarly, for both quark and hadronic stars, we have employedtwo resolutions in the case of merging binaries, i.e., dx =0 . M (cid:12) (369 m) and dx = 0 . M (cid:12) (236 m) , which we willindicate as “low” and “high”, respectively. Note that a reso-lution of dx = 0 . M (cid:12) (236 m) , and even lower ones, (seee.g., [121, 122]) have been used routinely in simulations ofbinary neutron stars and has been shown to be sufficient tocapture accurately the dynamics of the system. III. RESULTS: ISOLATED QUARK STARS
As a first but very important test of the correct implementa-tion of the rescaling procedure presented in the previous Sec- tion, we next discuss the results from the evolution of isolatedand oscillating quark stars, comparing with the correspond-ing evolutions obtained with fully hadronic stars. The analyti-cal formulation and the numerical computation of the spectralproperties of oscillating compact stars is a very old one andhas been studied in detail in the literature [123–126], also inthe case of quark stars [86, 88]. Similarly, the investigation ofthe dynamical response of a relativistic star to a perturbationhas been studied extensively in the literature (see, e.g., [127–129] for some of the initial works).Because this serves only as a test and not to study the spec-tral properties of such stars, the oscillations we have consid-ered are simple quasi-radial oscillations of initially static stars,whose perturbative eigenfrequencies of can be obtained easilyby solving two ordinary differential equations [see Eqs. (11)–(12) in Ref. [88]] that are coupled with the equations of stellarequilibrium, i.e., the Tolmann-Oppenheimer-Volkoff (TOV)equations (see, e.g., [105]). To make the comparison mean-ingful, both the isolated MIT2cfl quark stars and the hadronicDD2 stars have the same mass of M = 1 . M (cid:12) . Solving theperturbation equations together with the TOV equations, weobtain the eigenfrequencies F i , for different mode numbers i . In practice, we limit ourselves to the fundamental model F and the first overtone F . Note that both the static modelsand the eigenfrequencies are computed after introducing thecrustal prescription in the EOS, as described in the previousSection. For the dynamical evolutions, on the other hand, weintroduce an initial purely radial perturbation in the radial ve-locity with an eigenfunction corresponding to an n = 1 , (cid:96) = 0 oscillation mode under the Cowling approximation (i.e., witha spacetime being held fixed) [130], and an amplitude of .As customary, we report in Fig. 2 the normalised variationin the stellar central density ρ c ( t ) for either the quark star (leftpanel) or for the hadronic star (right panel). Curves of dif-ferent shading refer to different resolutions, with the high- f [kHz]0 . . . . . . P o w e r Sp ec tr u m F F F Quark Star F F F Quark Star F F F Quark Starlowmediumhigh . . . . . . . f [kHz]0 . . . . . . P o w e r Sp ec tr u m F F F Hadronic Star F F F Hadronic Star F F F Hadronic Star lowmediumhigh . . . . . . FIG. 3. Power spectral densities of the timeseries in Fig. (2) when compared with the perturbative predictions for the eigenfrequencies,reported as vertical blue dashed lines. The left panel refers to a quark star, while the right one to a hadronic star; the same convention isfollowed for the different types of lines. In both cases, the insets show a magnification near the frequency of the F mode. Note that theaccuracy of the match increases as the resolution is increased. resolution data being shown with the darkest shade. Notethat, especially for the quark star, the behaviour of the oscilla-tions shows aspects of strong nonlinearity, i.e., deviation froma simple harmonic oscillation, during the first few millisec-ond, but these then disappear rather soon and the oscillationbecomes rather regular. As to be expected, these nonlinearfeatures becomes more marked as the resolution is increasedand the deviations from a perfect conservation of rest-massare (cid:46) − % .Note also the imprint of different resolutions is more ev-ident in the quark star than in the hadronic one. A phasedifference can in fact be measured in the low- and medium-resolution evolutions for the quark star, which is instead ab-sent for the hadronic star. This is again easy to interpret asdue to the much sharper density jump at the stellar surface inthe case of a quark star, which requires therefore a higher res-olution to produce a consistent behaviour. However, for bothstars, the low resolution is sufficient to capture a convergentbehaviour and the accuracy simply increases as the resolu-tion is increased. Finally, it should be noted that while thequark star exhibits oscillations that are almost symmetric withrespect to the unperturbed value, the hadronic one does not.This is due to the fact that for the latter, the oscillations areaccompanied by a significant contribution from higher-ordermodes (see discussion below).Figure 3 helps to compare and contrast the evolutions ofquark and hadron stars by taking the Fourier transforms ofthe timeseries in Fig. (2) and comparing the resulting powerspectral densities (PSDs) with the perturbative prediction forthe eigenfrequencies, reported as vertical blue dashed lines.In this way, it is apparent that the match between the dynam-ical and the perturbative results is worse for the quark star,especially at low resolutions (see the insets). However, as theresolution is increased, both evolutions nicely converge to the perturbative values and differ from them with comparable er-rors, namely, . ( . ) and . ( . ) for the quark andhadron stars at high (medium) resolution, respectively. Notealso that the low-resolution of the quark star has a consider-able amount of noise near the position of the F frequency, butthat this noise decreases with increasing resolution, leading tothe clear appearance of an overtone in the high-resolution run.Finally, note that the amount of power in the two overtonesis much larger in the case of the hadronic star, which exhibitsa clear second overtone. This behaviour in indeed visible al-ready in the timeseries reported in Fig. 2 (right panel), whichexhibits rapid variations near the maxima and minima of theoscillations.In summary, the results presented here on the dynamics ofisolated quark stars show that the approach presented in Sec.II A to handle the stellar surface is not only properly imple-mented, but that it also leads to evolutions that are able toreproduce the spectral properties of quark stars and yield nu-merically consistent evolutions. At the same time, they showthat, in general, quark stars require higher resolutions than theones that would be needed for a hadronic star with compara-ble properties in terms of mass and tidal deformability. Onceagain, this is the inevitable drawback of having to model starswith a strong self-confinement and sharper density profiles atthe surface. IV. RESULTS: BINARY QUARK STARS
In this section we discuss the dynamics of a binary quark-star merger and contrast it with the corresponding mergerobtained with a hadronic EOS. Before entering in the de-tails, and notwithstanding that these are the first fully general-relativistic simulations of the merger of SQM compact ob-
FIG. 4. Small-scale (i.e.,
110 km ) views of the rest-mass density at different but characteristic times during the merger of the two classes ofbinaries [i.e., before merger (upper-left), the merger time (upper-right),
10 ms (lower left), and
20 ms (lower right) after merger]. Eachpanel offers views of the ( x, y ) plane (bottom parts) and of the ( x, z ) plane (top parts), respectively. Note that, for each panel, the left partrefers to quark stars, while the right one to the hadronic stars. jects, it is useful to stress what are the limitations of our ap-proach. First, the treatment of the thermal part of the EOS isnecessarily approximate and modeled with a hybrid EOS ap-proach for both EOSs [cf. Eqs. (8)–(10)]. As a result of this(forced) choice, the modelling of neutrinos, e.g., via a leakage[131], or more advanced approaches [50, 132], is not possible.Second, the modelling of the quark-evaporation processes inthe ejected matter, is de-facto ignored. While both processesare expected to have a significant impact on the electromag-netic kilonova signal and on the nucleosynthesis, they are notexpected to play an important role in the gravitational-wavesignal, which can therefore be considered robust. Third, forobvious computational costs, the two resolutions used herefor the binaries are rather low, and the highest one would cor- respond to the “low” resolution in the case of the isolated starsdiscussed in the previous Section. However, as demonstratedabove, such a resolution is sufficient to provide a sufficientlyaccurate description of the stellar behaviour given that therelative error in the PSD for the low-resolution quark star is (cid:46) . . More importantly, by comparing and contrastingquark and hadronic stars with similar properties we can easilycapture – even at these low resolutions – the systematic dif-ferences between the two stellar types, which we believe aretherefore robust. FIG. 5. Same as in Fig. 4, but for the entropy per baryon.
A. Overview of the dynamics
The main features of the dynamics of the binary systems ofquark and hadronic stars are rather similar. Before the merger,the stars in the binary systems inspiral towards each other asa result of energy and angular moment loss via gravitational-wave emission. During this stage, the matter effects are domi-nated by the tidal deformability of the two stars. This stageis shown in the top-left panels of Figs. 4–5, which reportthe small-scale (i.e.,
110 km ) views of the rest-mass density(Fig. 4) or of the entropy per baryon (Fig. 5). Differentpanels refer to different but characteristic times during themerger [i.e., before merger (upper-left), the merger time(upper-right),
10 ms (lower left), and
20 ms (lower right) af-ter merger] and offer views of the ( x, y ) plane (bottom parts)and of the ( x, z ) plane (top parts), respectively. More impor- tantly for our comparison, for each panel, the left part refersto quark stars, while the right one to the hadronic stars. Whilenot straightforward to analyse, Figs. 4–5 are rich of informa-tion and help obtain a comprehensive overview of the proper-ties of the dynamics of the two binaries.It should be noted that already in the inspiral phase, thequark-star binary loses a considerably smaller amount of ma-terial from the surface. While this matter is not significantboth in terms of gravitational-wave emission and nucleosyn-thesis (the material has very small velocities, is distributedessentially isotropically, and is gravitationally bound in goodpart), it already shows a feature we will encounter again be-low.As the two stars approach each other, the amplitude of theemitted gravitational waves increases, achieving a maximumwhen the stars merge. At this point in time, intense tidal forcesand strong shocks induce a significant and dynamical ejection0of matter (see upper-right panels of Fig. 4). After merger, adifferentially rotating object is produced whose mass is largerthan the maximum mass of rotating configuration and that isnormally referred to as hypermassive neutron star (HMNS) inthe case of hadronic EOSs, but that should dubbed hypermas-sive quark star (HMQS) in the case of the MIT2cfl EOS (seelower-left and lower-right panels of Fig. 4).The dynamics of this post-merger phase is again rathersimilar between the two EOSs and follows well-known be-haviours in terms of dynamical ejecta, tidal forces and shockheating (see, e.g., [13, 15, 133]). In particular, the dynamicalejection of matter takes place mostly at low latitudes (i.e., nearthe equatorial plane, as it can be seen in the ( x, z ) sectionsof Fig. 4). On the other hand, the material at high latitudes(i.e., near the polar directions, as it can be seen in the ( x, z ) sections of Figs. 5) is considerably hotter and with a higherentropy per baryon. Furthermore, strong shocks taking placein the HMQS and in the HMNS lead to a pulsed dynamicalmass ejection, which is apparent in the ( x, y ) sections of Figs.4, 5, where it appears in terms of a striped spiral-arm structure.Essentially all of this matter is gravitationally unbound and bybeing both hot and with high entropy, it will be involved in thesubsequent nucleosynthetic processes that cannot be modeledhere.Overall, when inspecting the full set of Figs. 4–5, it is ap-parent that although no major qualitative difference emerges –as one would have expected given the very similar bulk prop-erties of the components in the two binaries – there are somequantitative differences in the dynamics of quark-star binariesand hadron-star binaries. Postponing a more detailed compar-ison to Sec. IV C, we can already appreciate a first importantresult of our comparative study: mass loss is significantly sup-pressed – up to a factor of 30 – in quark binaries when com-pared to hadronic binaries. On the other hand, the differencesin temperature and entropy are much smaller and of a factorof a few, at most.
B. Gravitational-wave emission
The emission of gravitational waves from merging binariesof compact objects is obviously one of the most interestingoutcomes of these processes and can be used to extract impor-tant information on the status of matter when these compactobjects are stars. It is therefore natural to investigate whetherthe gravitational-wave signal computed in our simulations canbe exploited to distinguish the dynamics of strange quarksstars from that of hadronic stars.As customary, we can compute the strongest component ofthe gravitational-wave signal by extracting the (cid:96) = m = 2 mode of the Weyl-curvature scalar ψ from our simulations,so that the GW amplitudes in the corresponding mode and inthe two polarization + and × can be obtained by integrating ψ twice in time h − ih × = (cid:90) t −∞ dt (cid:48) (cid:90) t (cid:48) −∞ dt (cid:48)(cid:48) ( ψ ) . (21) − − − h + × [ M p c ] quark stars (high)hadronic stars (high)quark stars (low)hadronic stars (low) − − t − t merg [ms]024 f [ k H z ] FIG. 6.
Top panel: gravitational-wave strain in the (cid:96) = 2 = m modeand in the + polarisation, h . Curves of different shading referto different resolutions, with the high-resolution data being shownwith the darkest shade. The vertical black line marks the mergertime, which corresponds to the time of maximum amplitude, and isused to align the waveforms in phase. Bottom panel: instantaneousgravitational-wave frequencies relative to the strain in the top panel.The same convention is followed for the line types.
Figure 6, reports in its top panel the gravitational-wavestrains h for both the quark-star binary (red lines) and thehadron-star binary (black). As in previous figures, differentshades refer to the two resolutions employed, with the highresolutions being indicated with a full shading. The variouswaveforms are aligned at the time of the merger, that is, thetime when the amplitude reaches its first maximum. The bot-tom part of Fig. 6, on the other hand, reports, the correspond-ing evolution of the instantaneous gravitational-wave frequen-cies using the same convention for the line types. Overall,the gravitational-wave information provided in Fig. 6 remarksthat the inspiral part of the binary dynamics is very similarbetween the two types of stars, with the frequency evolutionsthat are essentially the same in this phase and up to the merger.However, small differences do develop after the merger and,as we will comment below, they can be used to extract impor-tant signatures on the occurrence of the merger of a quark-starbinary.In order to highlight the differences that emerge after themerger we follow the methodology presented in Ref. [63]to process the post-merger GW signal and perform a spec-tral analysis whose results are shown in Fig. 7. Firstof all, we collect the GW signal in the time interval t ∈ [ − , M (cid:12) ∼ [ − . , .
70] ms , and perform aFourier transformation with a symmetric time-domain Tukeywindow function, with parameter α = 0 . . This windowfunction helps eliminating spurious oscillations of the com-puted PSD. Furthermore, since we are interested only in thepost-merger signal, we employ a fifth-order high-pass But-terworth filter for the low-frequency part of the signal with a1 f [kHz] − − − l og h h ( f ) f / i [ H z − / , M p c ] adLIGO ET CE f f quark stars (high)hadronic stars (high)quark stars (low)hadronic stars (low) FIG. 7. Effective PSDs of the quark-star (red lines) and hadron-star binaries (black lines) at different resolutions. Also reported withfilled circles are the frequencies at merger f max , as well as the sen-sitivity curves of advanced and next-generation gravitational-wavedetectors such as advanced LIGO (adLIGO), the Einstein Telescope(ET), or Cosmic Explorer (CE). Note the presence of precise spectralfeatures, f and f , whose frequencies are shown in Table I. cutoff frequency set to be f cut = f ( t − t merg = − .
38 ms) +0 . . In this way we compute the effective PSD as [63] ˜ h ( f ) := (cid:115) | ˜ h + ( f ) | + | ˜ h × ( f ) | , (22)where ˜ h + ( f ) and ˜ h × ( f ) are the Fourier transforms of the totwo gravitational-wave strains. At the same time, it is possibleto compute signal-to-noise ratio (SNR) as SNR := (cid:34)(cid:90) ∞ | h ( f ) f / | S h ( f ) dff (cid:35) / , (23)where S h ( f ) is the noise PSD of a given gravitational-wavedetector.Figure 7 reports the effective PSDs of the two binaries(red and black lines for the quark-star and the hadron-starbinaries, respectively) at the two different resolutions (darkand light shading for the high and low resolutions, respec-tively), and the frequencies at merger f max (filled circles;these frequencies are also used to match the amplitudes ofthe PSDs obtained at different resolutions). Also reportedin Fig. 7 are the sensitivity curves of advanced and next-generation gravitational-wave detectors such as advancedLIGO (adLIGO), the Einstein Telescope (ET) [134, 135],or Cosmic Explorer (CE) [136]. In this way, it is possi-ble to appreciate the presence of precise spectral frequencies,i.e., gravitational-wave spectral lines, which are named as f and f after the convention in [137], and whose values areshown in Table I. Also reported as filled solid circles on top ofthe various PSDs are the gravitational-wave frequencies at the merger (i.e., at amplitude maximum) f max . Note that a thirdfrequency, f can be obtained from the approximate relation f ≈ f + f , which models the dynamics of the two stellarcores in terms and repeated collisions and bounces while theHMNS/HMQS rotates.As first pointed out in Ref. [59], the values of f max canbe used as an important proxy of the EOS. Furthermore, thisfrequency was shown to follow a quasi-universal behaviourin terms of the tidal deformability. On the other hand, bothFig. 6 and Fig. 7 clearly indicate that the values of f max arerather similar for both quark and hadronic stars. This is be-cause, although the radii of the stars in the two binaries arerather different, i.e., R QS (cid:39) . and R HS (cid:39) . for the quark (QS) and hadron star (HS), respectively (see Ta-ble I), the corresponding tidal deformabilities are very similar,i.e., Λ QS (cid:39) and Λ HS (cid:39) . Indeed, in Ref. [64] and us-ing a large set of purely hadronic EOSs, a universal relationwas found between the f max and tidal deformability , namely log (cid:18) f max Hz (cid:19) ≈ a + a Λ / − log (cid:32) MM (cid:12) (cid:33) , (24)where a = 4 . , a = − . , and M represent the aver-age mass at infinite separation (in our simulations and for bothbinaries, M = 1 . M (cid:12) ).The quasi-universal relation (24) is shown with a greensolid line in the left panel of Fig. 8 as a function of the tidaldeformability Λ and in the right panel as a function of the av-erage stellar compactness C := M /R , where R is the averageradius. We recall that this is possible because another quasi-universal relation exists relating the tidal deformability Λ withthe stellar compactness C , namely [138, 139] C = 0 . − . . . (25)Also shown in Fig. 8 with red and black symbols are the val-ues of the measured gravitational-wave frequencies at merger f max for binary quark stars (red filled circles) and binaryhadronic stars (black filled circles), respectively. When com-paring the numerical results with the expected quasi-universalrelation with the tidal deformability (left panel of Fig. 8), itis clear that the match is very good despite the very differentnature of the two binaries and, in particular, the large differ-ence in radii. This result may appear bizarre at first as onewould expect that the gravitational-wave frequency at mergeris directly related to the stellar radius and indeed the “con-tact” frequency is normally computed as the Keplerian fre-quency when the two stars have a separation which is twicethe stellar radius f cont := C / / (2 πM ) [137, 140]. However,as recently pointed out in Ref. [141], the gravitational-wavefrequency at merger f max is closely related to the quadrupolar( (cid:96) = 2 ) F -mode oscillation ( F F ) of nonrotating and rotat-ing models; furthermore, F F is can also be expressed quasi-universally in terms of the tidal deformability, so that the re-sults shown in the left panel of Fig. 8 are essentially stating In Ref. [64], the original form of Eq. (24) was given in terms of the tidaldeformability coefficient κ T2 , which is trivially related to Λ and in the caseof equal-mass binaries as κ T2 = Λ .
500 1000 1500 2000Λ1 . . . . f [ k H z ] f max f f f max f f Bauswein+ (2019)Bauswein+ (2019) 0 .
14 0 .
16 0 . M/R . . . . f [ k H z ] f max f f f max f f Bauswein+ (2019)Bauswein+ (2019)
FIG. 8.
Left panel:
Quasi-universal relations for the frequency at merger f max and the two post-merger frequencies f and f when shown asa function of the tidal deformability Λ ([cf. Eqs. (24), (26), and (28)]). Shown with symbols are the values of the frequencies measured fromthe high-resolution simulations in the case of quark-star binaries (red symbols) or hadron-star binaries (black symbols). Right panel: same asin the left panel but when the quasi-universal relations are expressed in terms of the average stellar compactness C ([cf. Eqs. (27) and (29)]).In both panels the shading refers to the uncertainties in the relations and the orange lines report the functional fitting in Ref. [40]. that quark and hadronic stars with similar tidal deformabil-ity – and thus similar properties of their dense cores – havesimilar quadrupolar F F oscillation modes and, hence, similargravitational-wave frequency at merger f max . This is anotherimportant result of our comparative study: quark-star bina-ries have merger frequencies similar to those of hadron-starbinaries with comparable tidal deformability . In turn, this im-plies that it may be difficult to distinguish the two classes ofstars based only on the gravitational-wave signal during theinspiral.Note, however, that when expressed in terms of the stellarcompactness (see right panel of Fig. 8), the match betweenthe measured merger frequencies and the expectations fromthe quasi-universal relations is reasonably good in the case ofhadronic stars, but it is much worse for quark-star binaries.This hints to the fact that while expression (24) could be usedboth for hadronic and quark-star binaries, the correspondingexpression in terms of the stellar compactness needs to be cor-rected in the case of quark-star binaries. We will return to thispoint below.Figure 8 also shows additional quasi-universal relations forthe f and f frequencies, both as a function of the tidal de-formability (left panel) and of the average stellar compactness(right panel). For the first, low-frequency peak, the universalrelations are expressed respectively as [63, 64] f ≈ b + b C + b C + b C kHz , (26) f ≈ c + c Λ / + c Λ / + c Λ / kHz , (27)where b = − . , b = 727 . , b = − . , b =10989 . , c = 45 . , c = − . , c = 7 . , and c = − . . Similarly, for the second and most powerful feature in the post-merger spectrum, we use [64] f ≈ . − .
800 Λ / kHz . (28) f ≈ . − .
016 exp (cid:16) −√ . C + 4 . (cid:17) kHz , (29)where for the second relation (29) we have used expression(25). These frequencies are shown in Fig. 8 with green linesand shadings together with their uncertainties are estimated inRef. [64]; furthermore, they are complemented by the expres-sion for the f frequency proposed in Eq. (1) of Ref. [40],where it is dubbed f peak (orange line and shading).A comparison with numerical data clearly shows that theuniversal relations provide a rather accurate representation inthe case of hadronic stars, both when expressed in terms of thetidal deformability (left panel) and of the average stellar com-pactness (right panel). The match remains very good also interms of quark stars, but only when the universal relations aregiven in terms of the tidal deformability (in the case of the f frequencies, the relative differences are of . and . for binary quark stars and binary hadron stars, respectively).Furthermore, the match of the f frequency of the hadron-starbinary (black filled star) with the quasi-universal relation isonly marginally acceptable when using expressions (28) and(29), while it is below the value estimated by Ref. [40], whoseuncertainty bands are probably underestimated. However, thedifferences can be quite substantial when the numerical datafor binary quark stars is compared with expressions (26), (27)for the f frequency or with expressions (28), (29) for the f frequency.We can summarise these findings into another notable resultof our comparative investigation: quark-star binaries havepost-merger frequencies that obey quasi-universal relationsderived from hadron-star binaries in terms of the tidal de- formability, but not when expressed in terms of the averagestellar compactness. Hence, it may be difficult to distinguishthe two classes of stars based only on the post-merger fre-quencies and if no information on the stellar radius is avail-able . As conjectured above, this behaviour seems to suggestthat new universal relations in terms of the stellar compact-ness need to be derived when consider binaries of SQM. Thisconjecture can only be verified when simulations with otherEOSs for the SQM have been performed. C. Properties of the ejected material
As mentioned in the Introduction, besides gravitationalwaves, another important product of the merger of a binarysystem of compact stars is represented by the ejected matteras this plays an important role in subsequent nucleosynthesisand in the kilonova emission. We have also already remarkedthat the lack of a consistent treatment of the evaporation pro-cess of SQM to nucleons prevents us from establishing whatis the actual role played by the ejected matter from the binaryquark-star merger. We recall that the study of Ref. [90] hasconcluded that only a small amount of quark matter can sur-vive evaporation and thus survive to yield a very low densityof SQM in the galaxy. At the same time, the use of a hybrid-EOS approach to handle the thermal effects in the hadronicDD2 EOS – and that we have employed to achieve a consis-tent picture with the MIT2cfl EOS – prevents us from mak-ing considerations on the nucleosynthetic yields and kilonovaemission also for the binary hadron-star merger.Notwithstanding these limitations, we can nevertheless ob-tain interesting comparative measurements of the ejected mat-ter in terms of bulk properties, such as: the total amounts ofejected rest-mass for the two classes of stars and the corre-sponding distributions in velocity and entropy. To this scope,and as done routinely in this type of simulations, we place aseries of detectors in terms of spherical coordinate 2-spheresat different distances from the origin and measure the amountof matter that crosses such detectors. In this way, consid-ering gravitationally unbound the matter satisfying the so-called “geodesic” criterion, namely, matter whose covarianttime component of the four-velocity is u t < − (see [142]for a detailed discussion of various criteria for unbound mat-ter) we can measure the effective mass lost by the binaries.This is shown as a function of time from the merger in Fig. 9for a detector placed at a radius of M (cid:12) (cid:39)
750 km , wherelines of different colours follow the same convention of theprevious figures. When comparing the results of the ejectedmass from the quark-star binaries and from the hadron-starbinaries, the difference is apparent, with the former hav-ing M ej (cid:39) − M (cid:12) and the latter M ej (cid:39) × − M (cid:12) ,i.e., about 30 times more ejected mass.There are at least two explanations behind this difference.The first is that, as mentioned repeatedly, SQM is charac-terised by self-bound properties and generally harder to eject.Second, and more importantly, a considerable rest-mass dif-ference is present between SQM and hadronic matter. We re-call, in fact, that the average rest-mass per baryon of SQM is t − t merg [ms]10 − − − − − M e j [ M (cid:12) ] quark stars (high)hadronic stars (high) FIG. 9. Ejected rest-mass as a function of time as measured by adetector placed at a radius of M (cid:12) (cid:39)
750 km . Lines of differ-ent colours follow the same convention of the previous figures: redfor quark-star binaries and black for hadron-star binaries. Note thatthe data refers to the high-resolution simulations and that the ejectedmass for the quark-star binaries is significantly smaller.
33 MeV , which is significantly smaller – about a factor of 30– than the corresponding value for hadronic matter,
930 MeV .Proceeding further in our comparison, we report in Fig. 10,two distributions of the ejected matter in terms of the normof the three-velocity (left panel) and in terms of entropy perbaryon (right panel). Note that both distributions are rela-tive to the ejected matter and hence provide a measure ofthe fraction of the ejecta with given properties. Concentrat-ing first on the left panel, it is clear that the two distributionsare very similar and that differences appear only in the high-velocity tails, i.e., v (cid:38) . , which are larger by almost oneorder of magnitude in the case of binary hadron-star merg-ers. Although the actual amount of matter in these tails is tiny(i.e., (cid:46) − M (cid:12) ), they could play a role when interactingwith the interstellar medium and produce a re-brightening ofthe afterglow signal [143] (see [18] for a discussion about therole of fast ejecta).When considering instead the distribution in entropy (rightpanel of Fig. 10) – and bearing in mind that the these mea-surements are the same for the two classes of stars but equallyapproximate – it is possible to appreciate that quark-star merg-ers overall produce ejected matter with larger entropy. At thesame time, hadron-star mergers are able to eject matter withvery high entropy, i.e., s B (cid:38) k B / baryon , and about oneorder of magnitude more than quark-star mergers. These dif-ferences on the entropy distributions may have a potential ef-fects on the subsequent kilonovae observation, but the degreeto which this will happen remains unclear until the quark-evaporation mechanism is properly taken into account, moresophisticated EOSs are devised to estimate the impact of SQM4 . . . . . v [ c ]10 − − − − − M i e j / M e j quark stars (high)hadronic stars (high)
20 40 60 80 100 s B [ k B / baryon]10 − − − − − M i e j / M e j quark stars (high)hadronic stars (high) FIG. 10. Distributions of the ejected matter in terms of the norm of the three-velocity (left panel) and of the entropy per baryon (right panel).Lines of different colours follow the same convention of the previous figures: red for quark-star binaries and black for hadron-star binaries.Note that the data refers to the high-resolution simulations and that quark-star binaries have much reduced tails in both quantities. in binary mergers , and more systematical investigations arecarried out. V. CONCLUSION
Since the hypothesis that SQM is the ground state of mat-ter has been formulated more than four decades ago, a vastliterature has been developed to investigate this scenario andconsider its viability against the astronomical observations.In this way, a very large number of works have been pre-sented in which the possibility that SQM could lead to com-pact stars, i.e., quark stars, has been explored in the greatestdetail. Of this vast literature, however, only a very small frac-tion has been dedicated to the study of the dynamics of bi-nary systems of quark stars. The scarcity of studies of thisscenario, which is particularly striking after the detection ofgravitational waves and of electromagnetic counterparts fromthe merger of low-mass compact binaries such as GW170817,probably has a double origin. From an observational side,there is the expectation that a merger of binary quark starscannot be responsible of the observed kilonova emission inGW170817 and of the subsequent nucleosynthesis. From atheoretical side, on the other hand, the modelling of the inspi-ral and merger of a binary of quark stars is particularly chal-lenging as it is difficult to properly capture the large disconti-nuity that appears at the surface of these self-bound stars.The first of these difficulties has been recently addressedin a number of recent works that have invoked a quark-evaporation process into hadrons that could have taken place We note that the MIT2cfl EOS considered here and introduced in [42] doesnot include an electron fraction, making it thus impossible to account forthe evolution of the electron fraction in the ejected matter. as a result of the high temperatures reached after the merger.In this case, therefore, most of ejected SQM from the quark-star binary would have evaporated into nucleons and there-fore could have therefore contributed to the kilonova signal inAT2017gfo [90, 92, 93]. Addressing the second difficulty, onthe other hand, is the purpose of this work, where we have pre-sented a novel technique in which a very thin crust has beenadded to the EOS of a quark star and where the properties ofthe EOS have been suitably rescaled to produce a smooth andgradual change of the specific enthalpy across the crust.The introduction of this crust, whose rest-mass content isminute, i.e., ∼ × − M (cid:12) , and whose spatial extension isrestricted to two grid cells only, has been carefully tested byconsidering the oscillation properties of isolated quark stars.In this way, it was possible to show that the dynamical andsimulated response of the quark stars matches accurately theperturbative predictions and that the match becomes increas-ingly accurate as the numerical resolution is increased. Thisvalidation, which has been introduced numerous times in thepast when considering neutron stars, has levels of accuracythat are comparable with those obtained with hadronic stars.Using this novel technique we have been able to carry outthe first fully general-relativistic simulations of the merger ofbinary strange quark stars. In addition, in order to best assessthe feature of this merging process that are specific to quarkstars, we have carried out a systematic comparison of thedynamics of quark-star binaries with the corresponding be-haviour exhibited by a binary of hadron stars having the samemass and very similar tidal deformability, namely, a binary de-scribed by the DD2 EOS. In this way, it was possible to high-light several important differences between the SQM and thehadronic stars, which represent the main results of our work.First, dynamical mass loss is significantly suppressed – up toa factor of 30 – in quark binaries when compared to hadronicbinaries. Second, quark-star binaries have merger frequencies5similar to those of hadron-star binaries with comparable tidaldeformability. Hence, it may be difficult to distinguish thetwo classes of stars based only on the gravitational-wave sig-nal during the inspiral. Third, quark-star binaries have post-merger frequencies that obey quasi-universal relations derivedfrom hadron-star binaries in terms of the tidal deformability,but not when expressed in terms of the average stellar com-pactness. Hence, it may be difficult to distinguish the twoclasses of stars based only on the post-merger frequencies andif no information on the stellar radius is available. Fourth, thematter ejected from quark-star binaries has much smaller tailsin the velocity distributions, i.e., for v (cid:38) . ; this lack of fastejecta may be revealed by the corresponding lack of a radiore-brightening when the fast ejecta interact with the interstel-lar medium. Finally, while quark-star binaries eject materialthat, on average, has larger entropy per baryon, it also lacksthe important tail of very high-entropy material. Determin-ing whether these differences in the ejected will able to leavean imprint in the electromagnetic counterpart and nucleosyn-thetic yields remains unclear.The results presented here had to rely necessarily on a num-ber of simplifying assumptions and can therefore be improvedin a number of ways. First, by adopting EOSs for the SQMthat have a proper treatment of the temperature dependenceand a nonzero electron fraction. Doing so will allow us not only to accurately and self-consistently describe the thermo-dynamical evolution after the merger, but also to determinewhether the ejected material will lead to the observed chemi-cal abundances. Second, by employing a larger class of EOSsso that it will be possible to establish whether the frequencyat merger and the post-merger frequencies obey new univer-sal relations when expressed in terms of the stellar compact-ness. Finally, and most importantly, by adopting a detailedand quantitative description of the quark-evaporation mecha-nism, so that a consistent assessment can be made of the vi-ability of binary quark-star mergers as sources to the electro-magnetic counterpart in AT2017gfo. All of these improve-ments will be explored and presented in future works. ACKNOWLEDGMENTS
It is a pleasure to thank Matthias Hanauske, Jens Papen-fort, and Lukas Weih, for useful input and comments. Sup-port comes in part from “PHAROS”, COST Action CA16214;LOEWE-Program in HIC for FAIR; the ERC Synergy Grant“BlackHoleCam: Imaging the Event Horizon of Black Holes”(Grant No. 610058). The simulations were performed on theSuperMUC and SuperMUC-NG clusters at the LRZ in Garch-ing, on the LOEWE cluster in CSC in Frankfurt, and on theHazelHen cluster at the HLRS in Stuttgart. [1] The LIGO Scientific Collaboration and The Virgo Collabora-tion (LIGO Scientific Collaboration and Virgo Collaboration),Phys. Rev. Lett. , 161101 (2017), arXiv:1710.05832 [gr-qc].[2] B. P. Abbott, R. Abbott, T. D. Abbott, S. Abraham, F. Ac-ernese, K. Ackley, C. Adams, R. X. Adhikari, V. B. Adya,C. Affeldt, and et al., Astrophys. J. Lett. , L3 (2020),arXiv:2001.01761 [astro-ph.HE].[3] The LIGO Scientific Collaboration, the Virgo Collaboration,R. Abbott, T. D. Abbott, S. Abraham, F. Acernese, K. Ackley,C. Adams, R. X. Adhikari, V. B. Adya, and et al., Astrophys.J. Lett. , L44 (2020), arXiv:2006.12611 [astro-ph.HE].[4] 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, R. X. Adhikari, V. B. Adya,and et al. (LIGO Scientific Collaboration and Virgo Collabora-tion), Astrophys. J. Lett. , L12 (2017), arXiv:1710.05833[astro-ph.HE].[5] LIGO Scientific Collaboration, Virgo Collaboration, Gamma-Ray Burst Monitor, INTEGRAL, B. P. Abbott, R. Abbott,T. D. Abbott, F. Acernese, K. Ackley, C. Adams, T. Adams,P. Addesso, R. X. Adhikari, V. B. Adya, and et al. (LIGO Sci-entific Collaboration and Virgo Collaboration), Astrophys. J.Lett. , L13 (2017), arXiv:1710.05834 [astro-ph.HE].[6] R. Narayan, B. Paczynski, and T. Piran, Astrophys. J. Lett. , L83 (1992), astro-ph/9204001.[7] D. Eichler, M. Livio, T. Piran, and D. N. Schramm, Nature , 126 (1989).[8] L. Rezzolla, B. Giacomazzo, L. Baiotti, J. Granot, C. Kouve-liotou, and M. A. Aloy, Astrophys. J. Letters , L6 (2011),arXiv:1101.4298 [astro-ph.HE]. [9] E. Berger, Annual Review of Astron. and Astrophys. , 43(2014), arXiv:1311.2603 [astro-ph.HE].[10] S. Wanajo, Y. Sekiguchi, N. Nishimura, K. Kiuchi, K. Kyu-toku, and M. Shibata, Astrophys. J. , L39 (2014),arXiv:1402.7317 [astro-ph.SR].[11] A. Perego, S. Rosswog, R. M. Cabez´on, O. Korobkin,R. K¨appeli, A. Arcones, and M. Liebend¨orfer, Mon. Not.R. Astron. Soc. , 3134 (2014), arXiv:1405.6730 [astro-ph.HE].[12] O. Just, A. Bauswein, R. A. Pulpillo, S. Goriely, andH.-T. Janka, Mon. Not. R. Astron. Soc. , 541 (2015),arXiv:1406.2687 [astro-ph.SR].[13] L. Bovard, D. Martin, F. Guercilena, A. Arcones, L. Rez-zolla, and O. Korobkin, Phys. Rev. D , 124005 (2017),arXiv:1709.09630 [gr-qc].[14] F.-K. Thielemann, M. Eichler, I. V. Panov, and B. Wehmeyer,Annual Review of Nuclear and Particle Science , 253(2017), arXiv:1710.02142 [astro-ph.HE].[15] D. Radice, A. Perego, K. Hotokezaka, S. A. Fromm,S. Bernuzzi, and L. F. Roberts, Astrophys. J. , 130 (2018),arXiv:1809.11161 [astro-ph.HE].[16] B. D. Metzger, Living Reviews in Relativity , 3 (2017),arXiv:1610.09381 [astro-ph.HE].[17] S. De and D. Siegel, arXiv e-prints , arXiv:2011.07176 (2020),arXiv:2011.07176 [astro-ph.HE].[18] E. R. Most, L. J. Papenfort, S. Tootle, and L. Rezzolla, arXive-prints , arXiv:2012.03896 (2020), arXiv:2012.03896 [astro-ph.HE].[19] B. Margalit and B. D. Metzger, Astrophys. J. Lett. , L19(2017), arXiv:1710.05938 [astro-ph.HE].[20] A. Bauswein, O. Just, H.-T. Janka, and N. Stergioulas, As- trophys. J. Lett. , L34 (2017), arXiv:1710.06843 [astro-ph.HE].[21] L. Rezzolla, E. R. Most, and L. R. Weih, Astrophys. J. Lett. , L25 (2018), arXiv:1711.00314 [astro-ph.HE].[22] M. Ruiz, S. L. Shapiro, and A. Tsokaros, Phys. Rev. D ,021501 (2018), arXiv:1711.00473 [astro-ph.HE].[23] E. Annala, T. Gorda, A. Kurkela, and A. Vuorinen, Phys. Rev.Lett. , 172703 (2018), arXiv:1711.02644 [astro-ph.HE].[24] D. Radice, A. Perego, F. Zappa, and S. Bernuzzi, Astrophys.J. Lett. , L29 (2018), arXiv:1711.03647 [astro-ph.HE].[25] E. R. Most, L. R. Weih, L. Rezzolla, and J. Schaffner-Bielich,Phys. Rev. Lett. , 261103 (2018), arXiv:1803.00549 [gr-qc].[26] I. Tews, J. Carlson, S. Gandolfi, and S. Reddy, Astrophys. J. , 149 (2018), arXiv:1801.01923 [nucl-th].[27] S. De, D. Finstad, J. M. Lattimer, D. A. Brown, E. Berger, andC. M. Biwer, Physical Review Letters , 091102 (2018),arXiv:1804.08583 [astro-ph.HE].[28] B. P. Abbott, R. Abbott, T. D. Abbott, F. Acernese, K. Ack-ley, C. Adams, T. Adams, P. Addesso, R. X. Adhikari, V. B.Adya, and et al. (LIGO Scientific Collaboration and VirgoCollaboration), Physical Review Letters , 161101 (2018),arXiv:1805.11581 [gr-qc].[29] M. Shibata, E. Zhou, K. Kiuchi, and S. Fujibayashi, Phys. Rev.D , 023015 (2019), arXiv:1905.03656 [astro-ph.HE].[30] S. Koeppel, L. Bovard, and L. Rezzolla, Astrophys. J. Lett. , L16 (2019), arXiv:1901.09977 [gr-qc].[31] A. Nathanail, E. R. Most, and L. Rezzolla, Astrophys. J. Lett,in press , arXiv:2101.01735 (2021), arXiv:2101.01735 [astro-ph.HE].[32] F. J. Fattoyev, J. Piekarewicz, and C. J. Horowitz, PhysicalReview Letters , 172702 (2018), arXiv:1711.06615 [nucl-th].[33] V. Paschalidis, K. Yagi, D. Alvarez-Castillo, D. B. Blaschke,and A. Sedrakian, Phys. Rev. D , 084038 (2018),arXiv:1712.00451 [astro-ph.HE].[34] G. F. Burgio, A. Drago, G. Pagliara, H.-J. Schulze, and J.-B.Wei, Astrophys. J. , 139 (2018).[35] G. Monta˜na, L. Tol´os, M. Hanauske, and L. Rezzolla, Phys.Rev. D , 103009 (2019), arXiv:1811.10929 [astro-ph.HE].[36] R. O. Gomes, P. Char, and S. Schramm, Astrophys. J. , 139(2019), arXiv:1806.04763 [nucl-th].[37] C.-M. Li, Y. Yan, J.-J. Geng, Y.-F. Huang, and H.-S. Zong,Phys. Rev. D , 083013 (2018), arXiv:1808.02601 [nucl-th].[38] J. J. Li, A. Sedrakian, and M. Alford, Phys. Rev. D ,063022 (2020).[39] E. R. Most, L. J. Papenfort, V. Dexheimer, M. Hanauske,S. Schramm, H. St¨ocker, and L. Rezzolla, Physical ReviewLetters , 061101 (2019), arXiv:1807.03684 [astro-ph.HE].[40] A. Bauswein, N.-U. F. Bastian, D. B. Blaschke, K. Chatziioan-nou, J. A. Clark, T. Fischer, and M. Oertel, Physical ReviewLetters , 061102 (2019), arXiv:1809.01116 [astro-ph.HE].[41] L. R. Weih, M. Hanauske, and L. Rezzolla, Phys. Rev. Lett. , 171103 (2020), arXiv:1912.09340 [gr-qc].[42] E.-P. Zhou, X. Zhou, and A. Li, Phys. Rev. D , 083015(2018), arXiv:1711.04312 [astro-ph.HE].[43] A. Drago and G. Pagliara, Astrophys. J. Lett. , L32 (2018),arXiv:1710.02003 [astro-ph.HE].[44] M. Shibata and K. Ury¯u, Phys. Rev. D , 064001 (2000), gr-qc/9911058.[45] M. D. Duez, Y. T. Liu, S. L. Shapiro, and B. C. Stephens, Phys.Rev. D , 104030 (2004), astro-ph/0402502.[46] M. Anderson, E. W. Hirschmann, L. Lehner, S. L. Liebling,P. M. Motl, D. Neilsen, C. Palenzuela, and J. E. Tohline, Phys. Rev. D , 024006 (2008), arXiv:0708.2720 [gr-qc].[47] L. Baiotti, B. Giacomazzo, and L. Rezzolla, Phys. Rev. D ,084033 (2008), arXiv:0804.0594 [gr-qc].[48] A. Bauswein and H.-T. Janka, Phys. Rev. Lett. , 011101(2012), arXiv:1106.1616 [astro-ph.SR].[49] S. Bernuzzi, T. Dietrich, W. Tichy, and B. Br¨ugmann, Phys.Rev. D , 104021 (2014), arXiv:1311.4443 [gr-qc].[50] D. Radice, F. Galeazzi, J. Lippuner, L. F. Roberts, C. D. Ott,and L. Rezzolla, Mon. Not. R. Astron. Soc. , 3255 (2016),arXiv:1601.02426 [astro-ph.HE].[51] Y. Sekiguchi, K. Kiuchi, K. Kyutoku, M. Shibata,and K. Taniguchi, Phys. Rev. D , 124046 (2016),arXiv:1603.01918 [astro-ph.HE].[52] L. Lehner, S. L. Liebling, C. Palenzuela, O. L. Caballero,E. O’Connor, M. Anderson, and D. Neilsen, Classical andQuantum Gravity , 184002 (2016), arXiv:1603.00501 [gr-qc].[53] M. Ruiz, V. Paschalidis, A. Tsokaros, and S. L. Shapiro, Phys.Rev. D , 124077 (2020), arXiv:2011.08863 [astro-ph.HE].[54] L. J. Papenfort, R. Gold, and L. Rezzolla, Phys. Rev. D ,104028 (2018), arXiv:1807.03795 [gr-qc].[55] V. Nedora, S. Bernuzzi, D. Radice, B. Daszuta, A. En-drizzi, A. Perego, A. Prakash, M. Safarzadeh, F. Schianchi,and D. Logoteta, arXiv e-prints , arXiv:2008.04333 (2020),arXiv:2008.04333 [astro-ph.HE].[56] L. Baiotti and L. Rezzolla, Rept. Prog. Phys. , 096901(2017), arXiv:1607.03540 [gr-qc].[57] V. Paschalidis, Classical and Quantum Gravity , 084002(2017), arXiv:1611.01519 [astro-ph.HE].[58] A. Bauswein and H.-T. Janka, Phys. Rev. Lett. , 011101(2012), arXiv:1106.1616 [astro-ph.SR].[59] J. S. Read, L. Baiotti, J. D. E. Creighton, J. L. Friedman, B. Gi-acomazzo, K. Kyutoku, C. Markakis, L. Rezzolla, M. Shi-bata, and K. Taniguchi, Phys. Rev. D , 044042 (2013),arXiv:1306.4065 [gr-qc].[60] A. Bauswein, N. Stergioulas, and H.-T. Janka, Phys. Rev. D , 023002 (2014), arXiv:1403.5301 [astro-ph.SR].[61] K. Takami, L. Rezzolla, and L. Baiotti, Phys. Rev. Lett. ,091104 (2014), arXiv:1403.5672 [gr-qc].[62] S. Bernuzzi, A. Nagar, S. Balmelli, T. Dietrich, and M. Ujevic,Phys. Rev. Lett. , 201101 (2014), arXiv:1402.6244 [gr-qc].[63] K. Takami, L. Rezzolla, and L. Baiotti, Phys. Rev. D ,064001 (2015), arXiv:1412.3240 [gr-qc].[64] L. Rezzolla and K. Takami, Phys. Rev. D , 124051 (2016),arXiv:1604.00246 [gr-qc].[65] S. Bose, K. Chakravarti, L. Rezzolla, B. S. Sathyaprakash,and K. Takami, Phys. Rev. Lett. , 031102 (2018),arXiv:1705.10850 [gr-qc].[66] A. Bauswein, S. Goriely, and H.-T. Janka, Astrophys. J. ,78 (2013), arXiv:1302.6530 [astro-ph.SR].[67] C. Palenzuela, S. L. Liebling, D. Neilsen, L. Lehner, O. L.Caballero, E. O’Connor, and M. Anderson, Phys. Rev. D ,044045 (2015), arXiv:1505.01607 [gr-qc].[68] T. Dietrich, S. Bernuzzi, M. Ujevic, and W. Tichy, Phys. Rev.D , 044045 (2017), arXiv:1611.07367 [gr-qc].[69] T. Dietrich, M. Ujevic, W. Tichy, S. Bernuzzi, andB. Br¨ugmann, Phys. Rev. D , 024029 (2017),arXiv:1607.06636 [gr-qc].[70] K. Kyutoku, K. Kiuchi, Y. Sekiguchi, M. Shibata,and K. Taniguchi, Phys. Rev. D , 023009 (2018),arXiv:1710.00827 [astro-ph.HE].[71] S. V. Chaurasia, T. Dietrich, M. Ujevic, K. Hendriks, R. Dudi,F. M. Fabbri, W. Tichy, and B. Br¨ugmann, arXiv e-prints ,arXiv:2003.11901 (2020), arXiv:2003.11901 [gr-qc]. [72] R. Gill, A. Nathanail, and L. Rezzolla, Astrophys. J. , 139(2019), arXiv:1901.04138 [astro-ph.HE].[73] A. Bauswein, H.-T. Janka, R. Oechslin, G. Pagliara, I. Sagert,J. Schaffner-Bielich, M. M. Hohle, and R. Neuh¨auser, Phys.Rev. Lett. , 011101 (2009), arXiv:0812.4248.[74] A. Bauswein, R. Oechslin, and H.-T. Janka, Phys. Rev. D ,024012 (2010), arXiv:0910.5169 [astro-ph.SR].[75] E. Witten, Phys. Rev. D , 272 (1984).[76] C. Alcock, E. Farhi, and A. Olinto, Astrophys. J. , 261(1986).[77] P. Haensel, J. L. Zdunik, and R. Schaefer, Astron. Astrophys. , 121 (1986).[78] N. Itoh, Progress of Theoretical Physics , 291 (1970).[79] R. X. Xu, S. I. Bastrukov, F. Weber, J. W. Yu, andI. V. Molodtsova, Phys. Rev. D , 023008 (2012),arXiv:1110.1226 [astro-ph.HE].[80] A. Drago, A. Lavagno, and I. Parenti, Astrophys. J. , 1519(2007), arXiv:astro-ph/0512652.[81] J. W. Yu and R. X. Xu, Mon. Not. R. Astron. Soc. , 489(2011), arXiv:1005.1792 [astro-ph.HE].[82] A. Drago, A. Lavagno, and G. Pagliara, Phys. Rev. D ,043014 (2014), arXiv:1309.7263 [nucl-th].[83] A. Drago and G. Pagliara, Phys. Rev. C , 045801 (2015),arXiv:1506.08337 [nucl-th].[84] A. Drago, A. Lavagno, G. Pagliara, and D. Pigato, EuropeanPhysical Journal A , 40 (2016), arXiv:1509.02131 [astro-ph.SR].[85] A. Drago and G. Pagliara, European Physical Journal A , 41(2016), arXiv:1509.02134 [astro-ph.SR].[86] G. Panotopoulos and I. Lopes, Phys. Rev. D , 083013(2017), arXiv:1709.06643 [gr-qc].[87] G. Wiktorowicz, A. Drago, G. Pagliara, and S. B. Popov, As-trophys. J. , 163 (2017), arXiv:1707.01586 [astro-ph.HE].[88] J. Bora and U. D. Goswami, arXiv e-prints , arXiv:2007.06553(2020), arXiv:2007.06553 [gr-qc].[89] X.-Y. Lai, Y.-W. Yu, E.-P. Zhou, Y.-Y. Li, and R.-X. Xu,Research in Astronomy and Astrophysics , 024 (2018),arXiv:1710.04964 [astro-ph.HE].[90] N. Bucciantini, A. Drago, G. Pagliara, and S. Traversi, arXive-prints , arXiv:1908.02501 (2019), arXiv:1908.02501 [astro-ph.HE].[91] X. Y. Lai, C. J. Xia, Y. W. Yu, and R. X. Xu, arXiv e-prints ,arXiv:2009.06165 (2020), arXiv:2009.06165 [astro-ph.HE].[92] R. De Pietri, A. Drago, A. Feo, G. Pagliara, M. Pasquali,S. Traversi, and G. Wiktorowicz, Astrophys. J. , 122(2019), arXiv:1904.01545 [astro-ph.HE].[93] J. E. Horvath, O. G. Benvenuto, E. Bauer, L. Paulucci,A. Bernardo, and H. R. Viturro, Universe , 144 (2019).[94] Z. Zhu, A. Li, and L. Rezzolla, Phys. Rev. D , 084058(2020), arXiv:2005.02677 [astro-ph.HE].[95] E. S. Fraga, R. D. Pisarski, and J. Schaffner-Bielich, Phys.Rev. D , 121702 (2001), arXiv:hep-ph/0101143 [hep-ph].[96] M. Alford, M. Braby, M. Paris, and S. Reddy, Astrophys. J. , 969 (2005), nucl-th/0411016.[97] A. Li, Z.-Y. Zhu, and X. Zhou, Astrophys. J. , 41 (2017).[98] P. B. Demorest, T. Pennucci, S. M. Ransom, M. S. E.Roberts, and J. W. T. Hessels, Nature , 1081 (2010),arXiv:1010.5788 [astro-ph.HE].[99] J. Antoniadis, P. C. C. Freire, N. Wex, T. M. Tauris, R. S.Lynch, and et al., Science , 448 (2013), arXiv:1304.6875[astro-ph.HE].[100] T. Damour and A. Nagar, Phys. Rev. D , 084035 (2009),arXiv:0906.0096 [gr-qc].[101] S. Postnikov, M. Prakash, and J. M. Lattimer, Phys. Rev. D , 024016 (2010), arXiv:1004.5098 [astro-ph.SR].[102] C. Kettner, F. Weber, M. K. Weigel, and N. K. Glendenning,Phys. Rev. D , 1440 (1995).[103] Y. F. Huang and T. Lu, Astron. Astrophys. , 189 (1997).[104] X. H. Wu, S. Du, and R. X. Xu, Mon. Not. R. Astron. Soc. , 4526 (2020), arXiv:2006.11514 [astro-ph.HE].[105] L. Rezzolla and O. Zanotti, Relativistic Hydrodynamics (Ox-ford University Press, Oxford, UK, 2013).[106] S. Typel, G. R¨opke, T. Kl¨ahn, D. Blaschke, and H. H. Wolter,Phys. Rev. C , 015803 (2010), arXiv:0908.2344 [nucl-th].[107] D. Radice, S. Bernuzzi, W. Del Pozzo, L. F. Roberts, and C. D.Ott, Astrophys. J. Lett. , L10 (2017), arXiv:1612.06429[astro-ph.HE].[108] M. Hanauske, J. Steinheimer, A. Motornenko, V. Vovchenko,L. Bovard, E. R. Most, L. J. Papenfort, S. Schramm, andH. St¨ocker, Particles , 44 (2019).[109] E. R. Most, L. J. Papenfort, and L. Rezzolla, Mon. Not. R. As-tron. Soc. , 3588 (2019), arXiv:1907.10328 [astro-ph.HE].[110] E. Gourgoulhon, P. Grandcl´ement, K. Taniguchi, J.-A. Marck,and S. Bonazzola, Phys. Rev. D , 064029 (2001), gr-qc/0007028.[111] D. Radice and L. Rezzolla, Astron. Astrophys. , A26(2012), arXiv:1206.6502 [astro-ph.IM].[112] D. Radice, L. Rezzolla, and F. Galeazzi, Mon. Not. R. Astron.Soc. L. , L46 (2014), arXiv:1306.6052 [gr-qc].[113] D. Radice, L. Rezzolla, and F. Galeazzi, Class. Quantum Grav. , 075012 (2014), arXiv:1312.5004 [gr-qc].[114] R. Haas, S. R. Brandt, W. E. Gabella, M. Gracia-Linares,B. Karakas¸, R. Matur, M. Alcubierre, D. Alic, G. Allen,M. Ansorg, M. Babiuc-Hamilton, L. Baiotti, W. Benger,E. Bentivegna, S. Bernuzzi, T. Bode, B. Bruegmann, M. Cam-panelli, F. Cipolletta, G. Corvino, S. Cupp, R. D. Pietri,P. Diener, H. Dimmelmeier, R. Dooley, N. Dorband, M. El-ley, Y. E. Khamra, Z. Etienne, J. Faber, T. Font, J. Frieben,B. Giacomazzo, T. Goodale, C. Gundlach, I. Hawke, S. Haw-ley, I. Hinder, S. Husa, S. Iyer, T. Kellermann, A. Knapp,M. Koppitz, P. Laguna, G. Lanferman, F. L¨offler, J. Masso,L. Menger, A. Merzky, J. M. Miller, M. Miller, P. Moesta,P. Montero, B. Mundim, A. Nerozzi, S. C. Noble, C. Ott,R. Paruchuri, D. Pollney, D. Radice, T. Radke, C. Reiss-wig, L. Rezzolla, D. Rideout, M. Ripeanu, L. Sala, J. A.Schewtschenko, E. Schnetter, B. Schutz, E. Seidel, E. Seidel,J. Shalf, K. Sible, U. Sperhake, N. Stergioulas, W.-M. Suen,B. Szilagyi, R. Takahashi, M. Thomas, J. Thornburg, M. To-bias, A. Tonita, P. Walker, M.-B. Wan, B. Wardell, H. Witek,M. Zilh˜ao, B. Zink, and Y. Zlochower, The einstein toolkit(2020), to find out more, visit http://einsteintoolkit.org.[115] A. Suresh and H. T. Huynh, Journal of Computational Physics , 83 (1997).[116] A. Mignone, P. Tzeferacos, and G. Bodo, Journal of Com-putational Physics , 5896 (2010), arXiv:1001.2832 [astro-ph.HE].[117] D. Alic, C. Bona-Casas, C. Bona, L. Rezzolla, and C. Palen-zuela, Phys. Rev. D , 064040 (2012), arXiv:1106.2254 [gr-qc].[118] D. Alic, W. Kastaun, and L. Rezzolla, Phys. Rev. D , 064049(2013), arXiv:1307.7391 [gr-qc].[119] D. Brown, P. Diener, O. Sarbach, E. Schnetter, and M. Tiglio,Phys. Rev. D , 044023 (2009), arXiv:0809.3533 [gr-qc].[120] E. Schnetter, S. H. Hawley, and I. Hawke, Classical and Quan-tum Gravity , 1465 (2004), arXiv:gr-qc/0310042 [gr-qc].[121] R. De Pietri, A. Feo, F. Maione, and F. L¨offler, Phys. Rev. D , 064047 (2016), arXiv:1509.08804 [gr-qc].[122] R. De Pietri, A. Feo, J. A. Font, F. L¨offler, F. Maione, M. Pasquali, and N. Stergioulas, Phys. Rev. Lett. , 221101(2018), arXiv:1802.03288 [gr-qc].[123] S. Chandrasekhar, Astrophys. J. , 417 (1964).[124] K. S. Thorne and A. Campolattaro, Astrop. J. , 591 (1967).[125] P. N. McDermott, H. M. van Horn, and C. J. Hansen, Astro-phys. J. , 725 (1988).[126] L. Rezzolla, in Summer School on Astroparticle Physics and Cosmology,ICTP Lecture Series, arXiv:gr-qc/0302025, Vol. 3 (2003).[127] J. A. Font, N. Stergioulas, and K. D. Kokkotas, Mon. Not. R.Astron. Soc. , 678 (2000), gr-qc/9908010.[128] J. A. Font, T. Goodale, S. Iyer, M. Miller, L. Rezzolla, E. Sei-del, N. Stergioulas, W.-M. Suen, and M. Tobias, Phys. Rev. D , 084024 (2002), gr-qc/0110047.[129] L. Baiotti, I. Hawke, P. J. Montero, F. L¨offler, L. Rezzolla,N. Stergioulas, J. A. Font, and E. Seidel, Phys. Rev. D ,024035 (2005), gr-qc/0403029.[130] T. G. Cowling, Mon. Not. R. Astron. Soc. , 367 (1941).[131] F. Galeazzi, W. Kastaun, L. Rezzolla, and J. A. Font, Phys.Rev. D , 064009 (2013), arXiv:1306.4953 [gr-qc].[132] L. R. Weih, H. Olivares, and L. Rezzolla, Mon. Not. R. Astron.Soc. , 2285 (2020), arXiv:2003.13580 [gr-qc].[133] M. Shibata and K. Hotokezaka, Annual Reviewof Nuclear and Particle Science , 41 (2019), https://doi.org/10.1146/annurev-nucl-101918-023625.[134] M. Punturo et al., Class. Quantum Grav. , 084007 (2010).[135] M. Punturo et al., Class. Quantum Grav. , 194002 (2010).[136] B. P. Abbott and et al., Classical and Quantum Gravity ,044001 (2017), arXiv:1607.08697 [astro-ph.IM].[137] K. Takami, L. Rezzolla, and L. Baiotti, Phys. Rev. Lett. ,091104 (2014), arXiv:1403.5672 [gr-qc].[138] K. Yagi and N. Yunes, Phys. Rep. , 1 (2017),arXiv:1608.02582 [gr-qc].[139] C. Raithel, F. ¨Ozel, and D. Psaltis, Astrophys. J. , L23(2018), arXiv:1803.07687 [astro-ph.HE].[140] T. Damour, A. Nagar, and L. Villain, Phys. Rev. D , 123007(2012), arXiv:1203.4352 [gr-qc].[141] H. Ho-Yin Ng, P. Chi-Kit Cheong, L.-M. Lin, and T. G. F. Li,arXiv e-prints , arXiv:2012.08263 (2020), arXiv:2012.08263[astro-ph.HE].[142] L. Bovard and L. Rezzolla, Classical and Quantum Gravity ,215005 (2017), arXiv:1705.07882.[143] K. Hotokezaka, K. Kiuchi, M. Shibata, E. Nakar, and T. Pi-ran, Astrophys. J.867