Simulating Binary Neutron Stars with Hybrid Equation of States: Gravitational Waves, Electromagnetic Signatures, and Challenges for Numerical Relativity

TL;DR

This paper explores how strange quark matter cores in neutron stars influence gravitational wave signals and electromagnetic signatures during mergers, highlighting numerical challenges and implications for interpreting GW170817 data.

Contribution

It introduces numerical relativity simulations of neutron star mergers with hybrid equations of state including strange quark matter, analyzing their gravitational and electromagnetic signatures.

Findings

Strong phase transitions affect numerical convergence order.

Existing waveform models are adequate within uncertainties but require higher resolution.

Electromagnetic signatures in the studied setup are too faint to explain GW170817's kilonova.

Abstract

The gravitational wave and electromagnetic signatures connected to the merger of two neutron stars allow us to test the nature of matter at supranuclear densities. Since the Equation of State governing the interior of neutron stars is only loosely constrained, there is even the possibility that strange quark matter exists inside the core of neutron stars. We investigate how strange quark matter cores affect the binary neutron star coalescence by performing numerical relativity simulations. Interestingly, the strong phase transition can cause a reduction of the convergence order of the numerical schemes to first order if the numerical resolution is not high enough. Therefore, an additional challenge is added in producing high-quality gravitational wave templates for Equation of States with a strong phase transition. Focusing on one particular configuration of an equal mass configuration consistent with GW170817, we compute and discuss the associated gravitational wave signal and some of the electromagnetic counterparts connected to the merger of the two stars. We find that existing waveform approximants employed for the analysis of GW170817 allow describing this kind of systems within the numerical uncertainties, which, however, are several times larger than for pure hadronic Equation of States, which means that even higher resolutions have been employed for an accurate gravitational wave model comparison. We also show that for the chosen Equation of State, quasi-universal relations describing the gravitational wave emission after the moment of merger seem to hold and that the electromagnetic signatures connected to our chosen setup would not be bright enough to explain the kilonova associated to GW170817.

Tables

Table 1. Table 1: Properties of the employed configuration. The columns refer to the individual masses in isolation, the baryonic mass of the individual stars, the stars’ compactnesses, the radii of the stars, the tidal deformability of the individual stars Λ A , B superscript Λ 𝐴 𝐵 {\Lambda}^{A,B} , the ADM-mass and angular momentum of the system, and the estimated residual eccentricity e 𝑒 e .

$M^{A,B}\left[M_{\odot }\right]$	$M_b^{A,B}\left[M_{\odot }\right]$	$C^{A,B}$	$R^{A,B}\left[\mathrm{km}\right]$	${\Lambda }^{A,B}$	$M_{\mathrm{ADM}}\left[M_{\odot }\right]$	$J_{\mathrm{ADM}}\left[M_{\odot }^2\right]$	$e$
1.365	1.527	0.18	11.0	378	2.705	8.324	$5\times {10}^{-3}$

Table 2. Table 2: Grid configurations employed in the simulations for the different resolutions. The columns refer to the simulation’s name, the number of moving boxes L mv superscript 𝐿 mv L^{\rm mv} , the number of grid points per direction in the static levels n 𝑛 n , the number of grid points per direction in the moving levels n mv superscript 𝑛 mv n^{\rm mv} , the grid spacing in the innermost level h 6 subscript ℎ 6 h_{6} (corresponding to l = 6 𝑙 6 l=6 ), and the grid spacing of the coarsest level h 0 subscript ℎ 0 h_{0} .

Name	$L^{\mathrm{mv}}$	$n$	$n^{\mathrm{mv}}$	$h_6$	$h_0$
R1	$4$	$192$	$64$	$0.240$	$15.36$
R2	$4$	$288$	$96$	$0.160$	$10.24$
R3	$4$	$384$	$128$	$0.120$	$7.68$
R4	$4$	$480$	$160$	$0.096$	$6.14$

Table 3. Table 3: Piecewise polytrope of the HyS EoS used in the simulations. The units are chosen so that the rest-mass density ρ 𝜌 \rho is in g / cm 3 g superscript cm 3 \rm g/cm^{3} , Γ Γ \Gamma is dimensionless and K 𝐾 K is such that the pressure p 𝑝 p is in g / cm 3 g superscript cm 3 \rm g/cm^{3} . The i 𝑖 i -th row refers to the polytrope p = K i ​ ρ Γ i 𝑝 subscript 𝐾 𝑖 superscript 𝜌 subscript Γ 𝑖 p=K_{i}\rho^{\Gamma_{i}} , for ρ i ≤ ρ ≤ ρ i + 1 subscript 𝜌 𝑖 𝜌 subscript 𝜌 𝑖 1 \rho_{i}\leq\rho\leq\rho_{i+1} .

$\rho$	$K$	$\Gamma$
$0.00000$	$3.99873\times {10}^{-8}$	$1.357$
$1.46231\times {10}^{14}$	$1.80658\times {10}^{-31}$	$3.005$
$2.60400\times {10}^{14}$	$3.76582\times {10}^{12}$	$0.000$
$4.66533\times {10}^{14}$	$1.00162\times {10}^{-167}$	$12.242$
$5.06058\times {10}^{14}$	$1.29766\times {10}^{-108}$	$8.222$
$6.07892\times {10}^{14}$	$2.18919\times {10}^{-43}$	$3.810$
$7.22396\times {10}^{14}$	$7.07726\times {10}^{-28}$	$2.766$
$8.50306\times {10}^{14}$	$7.86344\times {10}^{-21}$	$2.294$
$9.92358\times {10}^{14}$	$7.59597\times {10}^{-17}$	$2.028$
$1.14929\times {10}^{15}$	$2.58041\times {10}^{-14}$	$1.860$
$1.32183\times {10}^{15}$	$1.38093\times {10}^{-12}$	$1.746$
$1.51072\times {10}^{15}$	$2.40271\times {10}^{-11}$	$1.664$
$1.71669\times {10}^{15}$	$2.01349\times {10}^{-10}$	$1.604$
$1.94048\times {10}^{15}$	$1.02422\times {10}^{-9}$	$1.557$
$2.18282\times {10}^{15}$	$4.06158\times {10}^{-9}$	$1.518$

Equations

$$\Omega =-\frac{3}{4{\pi }^2}{a}_4{\mu }^4+\frac{3}{4{\pi }^2}{a}_2{\mu }^2+B+{\Omega }_e,$$

$$\frac{2}{3}{n}_u-\frac{1}{3}{n}_d-\frac{1}{3}{n}_s=n_e,$$

$${\mu }_d={\mu }_s={\mu }_u+{\mu }_e,$$

$$p$$

$$n_b$$

$$ϵ$$

$$p(ϵ)=\frac{1}{3}(ϵ-4B)-\frac{a_2^2}{12{\pi }^2{a}_4}\left(1+\sqrt{1+\frac{16{\pi }^2{a}_4}{a_2^2}(ϵ-B)}\right),$$

$$g_q(p)=g_h(p),$$

$$g=\frac{ϵ+p}{n_b}=3\mu .$$

$$p(g)=\frac{a_4}{108{\pi }^2}{g}^4-\frac{9a_2}{108{\pi }^2}{g}^2-B.$$

$$B>\frac{g^2}{108{\pi }^2}(a_4{g}^2-9a_2).$$

$$B_{\min }=\frac{g_{\min }^2}{108{\pi }^2}(a_4{g}_{\min }^2-9a_2).$$

$$g_{2f}>g_h,$$

$$B$$

$${\bar{B}}_{\min }$$

$$g_q(p)={\left\{\frac{9a_2}{2a_4}\left[1+\sqrt{1+\frac{(4\pi {)}^2{a}_4}{3a_2^2}(p+B)}\right]\right\}}^{1/2}.$$

$$g=\frac{ϵ+p}{n_b},$$

$${ϵ}_q=\frac{a_4}{36{\pi }^2}{g}_T^4-\frac{a_2}{12{\pi }^2}{g}_T^2+B,$$

$$p(\rho )=K_i{\rho }^{{\Gamma }_i},d\left(\frac{ϵ}{\rho }\right)=-pd\left(\frac{1}{\rho }\right),{\rho }_{i-1}<\rho <{\rho }_i$$

$$u=t-r_{\mathrm{ext}}-2M\ln (r_{\mathrm{ext}}/2M-1),$$

$$f_2=(3.27\pm 0.1)\mathrm{kHz},$$

Full text

\Title

Simulating Binary Neutron Stars with Hybrid Equation of States: Gravitational Waves, Electromagnetic Signatures, and Challenges for Numerical Relativity \AuthorHenrique Gieg1,∗, Tim Dietrich2, Maximiliano Ujevic1

\corresCorrespondence: henrique@gieg.com.br

1 Introduction

The observation of gravitational waves (GWs), GW170817, and electromagnetic (EM) signatures, AT2017gfo and GRB170817A, originating from the same astrophysical source have initiated a new era of multi-messenger astronomy Abbott et al. (2017a, b, c).

Because of their large compactness, neutron stars (NSs) contain matter at supranuclear densities, which makes them a perfect laboratory to determine matter under extreme conditions governing the NS interior. Today’s Equation of State (EoS) constraints inferred from GW170817 arise from the analysis of the GW signal, e.g., Abbott et al. (2017a); Dai et al. (2018); De et al. (2018); Abbott et al. (2019, 2018a, 2018b) or from a combined GW and EM analysis, e.g, Radice et al. (2018); Bauswein et al. (2017); Coughlin et al. (2018); Radice and Dai (2018); Coughlin et al. (2018). In addition, one can also perform statistical analyses of a large set of possible EoSs incorporating information from nuclear physics theory, e.g., Annala et al. (2018); Most et al. (2018).

The current analysis of GW170817 disfavors NSs with too large radii and tidal deformabilities, the exact state of the supranuclear matter inside the NS core is still unknown. While recent first-principle studies disfavor the presence of strange quark matter (SQM) inside NSs with masses $\lesssim 1.4M_{\odot}$, Ref. Annala et al. (2019), it is not clear if more massive NSs could have SQM cores or if the merger could undergo a strong phase transition Most et al. (2019); Bauswein et al. (2019).

In this article, we study possible effects of SQM during the binary coalescence, for this purpose, we solve the equations of general relativity combined with the equations of general relativistic hydrodynamics (GRHD) with the help of full 3+1 -dimensional (3 dimensions in space and 1 dimension in time) numerical relativity (NR) simulations. Such simulations allow a characterization of the GW and EM radiation. Over the last years, the NR community has made significant progress in simulating neutron star spacetimes on many fronts, e.g., extending the simulation’s accuracy, incorporating microphysical models in the simulations, and exploring different regions of the binary neutron star (BNS) parameter space Radice et al. (2014); Hotokezaka et al. (2015); Dietrich et al. (2017); Kiuchi et al. (2017); Dietrich et al. (2018); Rezzolla et al. (2011); Neilsen et al. (2014); Sekiguchi et al. (2015); Palenzuela et al. (2015); Foucart (2018); Palenzuela et al. (2015); Ruiz et al. (2018); Ciolfi et al. (2017); Kiuchi et al. (2018); Radice (2017); Shibata et al. (2017). We will show that standard techniques are capable of dealing with NSs containing hadronic and SQM. Up to our knowledge, most previous simulations with EoSs allowing strong phase transitions have studied the imprint of a SQM phase after the merger of the two stars Most et al. (2019); Bauswein et al. (2019), in addition, also the phenomenological implications of the two-families scenario (in which hadronic stars and pure strange quark stars coexist) has been recently analysed Pietri et al. (2019), but no numerical relativity inspiral simulation with stars containing SQM has been performed. Thus, the effect of employing a hybrid EoS, in which hadrons convert into strange quarks through a thermodynamical phase transition, on the inspiral, has to be investigated in more detail. We find that a strong phase transition impairs the NR inspiral evolution and that the convergence order of the simulations can drop to first order if the employed resolutions are not sufficient. While it might be possible to deal with phase transitions incorporating more sophisticated numerical techniques or higher resolutions, our work is (up to our knowledge) the first to point out this arising complication.

The structure of the article is as follows. First, we discuss how we construct a hybrid star EoS using the tdBag model, we continue by describing the employed setup and numerical techniques. Then, we present a first qualitative description of the binary coalescence and continue by discussing the GW emission during the inspiral and postmerger phase, as well as, the related ejecta outflow and expected EM counterparts. Unless otherwise stated, we employ geometric units, i.e., $G=c=1$. Physical units are sometimes given for better intuition.

2 A new hybrid Star Equation of State

2.1 tdBag model for strange quark matter

tdBag Farhi and Jaffe (1984) is a bag-like model inspired by the standard MIT bag model Chodos et al. (1974); Detar (1979)111For convenience we will use natural units in this section setting $\hbar=c=1$.. It incorporates features expected for cold strange quark matter composed of up, down, and strange quarks in free Fermi seas resulting from the deconfinement, i.e., the breaking of hadronic matter into its constituents once baryon densities compatible to the nuclear density $n_{0}\sim 0.16$ fm*-3* are reached. tdBag accounts for perturbative effects of strong interaction between quarks Alford et al. (2005) and the effects of color superconductivity at low temperatures Alford et al. (2008). There are many possible color superconducting phases, in this work we assume that quarks are in the color-flavor-locked phase (CFL) Lugones and Horvath (2002); Alford et al. (2008), in which the three quark flavors are symmetric (equal baryon density of up, down and strange) and exhibit superfluidity.

In the massless up and down quarks approximation $m_{u}=m_{d}=0$, the tdBag model is defined by the grand potential per volume ($\Omega$) as

[TABLE]

where $a_{4}$, $a_{2}$, and $B$ are free parameters and $\mu\equiv(\mu_{u}+\mu_{d}+\mu_{s})/{3}$ is the quark chemical potential, $\mu_{i}$ refers to the chemical potential of the $i$-th flavor, and $\Omega_{e}$ is the electrons’ contribution. Before proceeding to the construction of hybrid EoS containing a SQM core, we are going to briefly discuss the meaning of the parameters presented in Eq. (1):

1. The quartic coefficient $a_{4}$ can be interpreted as a correction to the pressure of the free Fermi sea arising from quantum chromodynamics (QCD). Fraga et al. Fraga et al. (2001) show that, to second order in the strong coupling constant, a reasonable value for this parameter is $a_{4}\sim 0.7$. The values of $a_{4}=0$ and $a_{4}=1$ correspond to maximum or none QCD interaction, respectively.

2. The quadratic coefficient $a_{2}$ accounts for the contribution of the strange quark mass and color superconductivity. The value of this parameter is given by $a_{2}=m_{s}^{2}-4\Delta^{2}$, where $m_{s}\sim 100$ MeV is the strange quark mass and $\Delta\sim 0-100$ MeV is the energy gap between free and Cooper pairs of quarks Alford et al. (2005); Lugones and Horvath (2002); Klahn and Fischer (2015).

3. The bag constant $B$, introduced phenomenologically into the theory to describe the deconfinement, is the energy density required to produce a “vacuumless” volume in which quarks exist.

4. The contribution of the electron potential ($\Omega_{e}$) for our purposes is negligible; cf. Ref. Pereira et al. (2018) and references therein. Nevertheless, if the SQM is not in the CFL phase Lugones and Horvath (2002) there must be electrons within the quark phase to ensure (i) local charge neutrality and (ii) chemical equilibrium. These two conditions may be expressed respectively as

[TABLE]

where $n_{i}$ denotes the number density. For SQM in the CFL phase, the quark matter is symmetric, thus, local charge neutrality is satisfied without the presence of electrons. In both cases, CFL ($\Omega_{e}=0$) or non-CFL phase (neglecting $\Omega_{e}$), Eq. (1) leads to an analytical EoS for SQM.

2.2 Thermodynamics and parameter analysis

In this section we constrain the parameter space $(a_{4},a_{2},B)$ presented in Eq. (1) assuming the absolutely stable strange quark matter hypothesis Bodmer (1971); Witten (1984) and thermodynamic properties. With the help of the grand potential $\Omega$, we found that the pressure $p$, the baryon number density $n_{b}$, and the energy density $\epsilon$ can be written as

[TABLE]

It is possible to obtain $p(\epsilon)$ from Eq. (4), (5) and (6), and the result is

[TABLE]

which is the analytical EoS for SQM parametrized by $a_{4}$, $a_{2}$, and $B$.

As long as we are dealing with a phenomenological model, it is of fundamental importance to constrain the parameter space within a region in which it can produce EoSs for stellar matter. To justify our next choices, we make a digression on the absolutely stable strange quark matter hypothesis (i.e. no strange matter will decay into hadronic matter under a thermodynamic process), proposed by Bodmer Bodmer (1971) and Witten Witten (1984). They consider the ground state of matter to be three Fermi seas of deconfined quarks in equal proportions of up, down and strange flavors. The reason for this is that the energy per baryon ($E/A$) in the ground state is smaller than that of the most stable pressureless nuclei, 56Fe with $E/A\sim 930$ MeV. Until now, no experimental observations could refute this hypothesis. Indeed, as deconfinement is expected to take place under extreme conditions such as the ones in the core of compact stars, many efforts have been made to uncover which observables in astrophysical phenomena could rule out or confirm the existence of SQM in stars.

In the case of a Hybrid Star (HyS), i.e., a star containing a SQM core surrounded by hadronic matter (HM), the deconfinement of HM into SQM must take place at a strictly positive pressure. Assuming that the deconfinement is a first order phase transition, i.e., the interface between the hadron and quark phases is characterized by a sharp discontinuity (no mixed phase is present). At the interface the Gibbs condition must be satisfied

[TABLE]

where $g_{q}$, $g_{h}$ are the Gibbs free energy per baryon of the quarks and hadrons, respectively, with the generic $g$ given by

[TABLE]

Eq. (8) represents chemical and mechanical equilibrium at the interface. Using Eqs. (9) and (4), one can write $p$ as a function of the Gibbs free energy

[TABLE]

To enforce the phase transition pressure to be strictly positive ($p>0$), the following relation must hold

[TABLE]

According to the absolutely stable SQM hypothesis the minimum of $g$ is $g_{\rm min}(p=0)\approx 930$ MeV. This value of $g$ sets $B_{\rm min}$, for which hadronic matter deconfines in three quark flavors, such that

[TABLE]

This equation defines the region on the $B$, $a_{2}$, $a_{4}$ space above which a transition from HM to SQM happens. On the other hand, the hadronic part of the HyS must not be in metastable equilibrium, i.e., one shall guarantee that it would not decay in 2-flavor quark matter. This condition can be ensured if

[TABLE]

where $g_{2f}$ is the Gibbs free energy per baryon for 2-flavor quark matter. A second similar condition for $B$ can be found Pereira et al. (2018) adapting Eq. (1) to the $2f$ case in the massless approximation ($m_{u}=m_{d}=0$), the result is

[TABLE]

Therefore, conditions (12) and (14) constrain the parameters space $(a_{4},a_{2},B)$ in three regions (see Fig. 1 for an example),

• (i)

$B<\bar{B}_{\rm min}$: not allowed region due to the fact that hadronic matter is more stable than quark matter, then we cannot use a SQM EoS.

• (ii)

$\bar{B}_{\rm min}<B<B_{\rm min}$: absolutely stable strange stars (SS). There is no phase transition HM $\rightarrow$ SQM and the star is entirely made of 3-flavor quarks.

• (iii)

$B>B_{\rm min}$: HyS. There exists a phase transition between HM and SQM for some energy density $\epsilon$.

The EoS for SS is given directly by Eq. (7). To obtain the EoS for HyS we must solve Gibbs condition (8) with the help of a hadronic EoS to find the deconfinement pressure $p_{T}$ and the corresponding energy densities $\epsilon_{h}$ and $\epsilon_{q}$ for hadronic and quark matter, respectively. Assuming a first order phase transition, the HyS EoS is then constructed as a piecewise function with the prescription

• (i)

$\epsilon\leqslant\epsilon_{h}$: HM EoS,

• (ii)

$\epsilon_{h}<\epsilon<\epsilon_{q}$: $p=p_{T}$,

• (iii)

$\epsilon\geqslant\epsilon_{q}$: SQM EoS.

If deconfinement produces a sharp discontinuity between the phases or an extended region in which hadrons and quarks are mixed, is an open question that relies primarily on how global electric charge neutrality is achieved by the star. This could happen in two scenarios: (i) assuming local charge neutrality, i.e., that the electric charge density is zero everywhere in the stellar interior or (ii) assuming that electric charge density is non-zero and distributed within the star such that the total charge is zero. Despite the latter being a relaxation of the first condition, it has profound implications for the star’s structure, e.g. the rise of Coulomb interactions between quarks in different geometrical shapes (drops, rods, and slabs) and the surrounding hadronic matter may lead to a mixed phase of quarks and hadrons extending along the stellar interior Glendenning (2000). On the other hand, if local charge neutrality is considered, quarks and hadrons separate into two phases with a sharp discontinuity interface between them. The existence of a mixed phase is favored if the Coulomb interactions were to lower the system’s energy Glendenning (2000) in comparison to maintain the interface by surface tension. Recently, it had been shown Lugones et al. (2013) that it is energetically favorable for SQM to sharp discontinuity interface if the quark matter is described by the Nambu-Jona-Lasinio model Nambu and Jona-Lasinio (1961).

2.3 Building hybrid star EoS

The recent observation of GWs from GW170817 established constraints on the NS EoS. The analysis of the data favored soft EoS instead of stiff EoS. For our analysis, we use for the HM phase the SLy EoS Douchin and Haensel (2001), which is currently in agreement with EoS constraints. Recalling that we assumed the deconfinement as a first order phase transition (at the sharp interface between quarks and hadrons) the Gibbs conditions (8) must hold for both phases.

We start writing $g$ as a function of $p$ for both phases. This is straightforward for the SQM, solving Eq. (10) for $g$, one obtains

[TABLE]

For the chosen hadronic EoS, we construct the Gibbs free energy per baryon with the definition

[TABLE]

from which we obtain an expression for $g_{h}\equiv g_{h}(p)$ by interpolating the hadronic EoS. The Gibbs condition (8) can be satisfied numerically by finding the root of $f(p)=g_{q}(p)-g_{h}(p)$. The result is the pressure $p_{T}$ at which the two curves $g(p)$ meet (see Fig. 2) and corresponds to the constant pressure along the phase transition. The value $g_{T}=g(p_{T})$ is the common Gibbs free energy per baryon. The next step is to obtain the energy density $\epsilon$ associated with $p_{T}$ (or $g_{T}$) for both EoS. In the case of SQM, from Eqs. (6) and (9), we obtain

[TABLE]

which sets the energy density at the phase transition for the SQM phase. Similar calculations in the hadronic phase must be carried out to find $\epsilon_{h}$. The final result is a HyS EoS produced by connecting the hadronic EoS to the SQM EoS Eq. (7) with $p=p_{T}$ for $\epsilon_{h}<\epsilon<\epsilon_{q}$, as shown in Fig. 2.

2.4 HyS EoS as piecewise polytropes

For the dynamical evolution of the binary HyS with the BAM code Bruegmann et al. (2008); Thierfelder et al. (2011); Dietrich et al. (2015); Bernuzzi and Dietrich (2016) and to calculate the initial condition using the SGRID code Tichy (2009), we express the EoS as a piecewise polytrope (PwP):

[TABLE]

with $K_{i}$ and $\Gamma_{i}$ being defined in the individual intervals $\rho_{0}<...<\rho_{i}<...<\rho_{n}$ such that the continuity of $p(\rho)$ is guaranteed.

The SLy part of our hybrid EoS uses the piecewise polytrope representation given in Ref. Read et al. (2009). For the SQM phase, Eq. (7) has to be interpolated as polytropes in rest mass density. This is done producing a tabulated version of the SQM EoS and dividing this part of the EoS into fifteen pieces from which we built the SQM piecewise polytrope. The complete HyS piecewise polytrope can be found in Appendix A. In Fig. 3 we compare the results of the mass-radius curve for an isolated TOV star obtained using the tabulated EoS and the piecewise polytrope EoS. For the gravitational mass used in this article, we find a difference less than $0.5\%$ in the star radius.

3 Configurations details and numerical methods

3.1 Hybrid star configuration

For our primary study of the binary HyS coalescence, we chose a setup which is close to the estimated binary properties of GW170817 (details are given in Tab. 1). Our configuration consists of two non-spinning stars with individual masses of $M=1.365M_{\odot}$. The total mass of $M=M^{A}+M^{B}=2.73M_{\odot}$ is consistent with the estimated total mass of GW170817 ($2.73_{-0.01}^{+0.04}[M_{\odot}]$, Ref. Abbott et al. (2019)). The tidal deformability of the individual stars is $\Lambda^{A,B}=378.4$ and the radii of the stars are $\sim 11.0\rm km$, which also agrees with current multi-messenger constraints arising from the analysis of the GW, the kilonova, and the short GRB Coughlin et al. (2018) of $R\in[10.9,13.6]\rm km$.

3.2 Numerical methods

We compute the initial configuration solving the conformal thin-sandwich equations Wilson and Mathews (1995); Wilson et al. (1996); York (1999) with the SGRID code Tichy (2009, 2012); Dietrich et al. (2015). SGRID uses pseudospectral methods to accurately compute spatial derivatives. We performed an iterative procedure to reduce the orbital eccentricity varying the system’s initial radial velocity and eccentricity parameter, as discussed in Moldenhauer et al. (2014); Kyutoku et al. (2014); Dietrich et al. (2015), until eccentricity is $\sim$ $10^{-3}$.

Dynamical evolutions are performed with the BAM code Bruegmann et al. (2008); Thierfelder et al. (2011); Dietrich et al. (2015); Bernuzzi and Dietrich (2016), which uses the method of lines employing finite differences for the spatial discretization of the metric variables and a high-resolution shock-capturing (HRSC) scheme for the computation of the matter fluxes of the general relativistic hydrodynamics equations. The spacetime is evolved with the Z4c evolution system Bernuzzi and Hilditch (2010); Weyhausen et al. (2012); Hilditch et al. (2013). The BAM code employs an adaptive mesh refinement strategy with a hierarchy of nested Cartesian grids. The grid spacing of each level is half the grid spacing of its coarser parent level. A number of inner levels can be moved dynamically to cover the NS’s during the evolution. In this article, we use $7$ refinement levels labeled $l=0,...,6$  from the coarsest to the finest. The detailed grid specifications for the employed resolutions are given in Tab. 2. The GRHD numerical fluxes are computed with a local Lax-Friedrich’s (LLF) scheme and the fluid’s characteristic variables are reconstructed using a fifth order weighted essentially non-oscillatory (WENOZ) algorithm Borges et al. (2008); Bernuzzi and Dietrich (2016). The same GRHD settings have been employed in Bernuzzi and Dietrich (2016); Dietrich et al. (2017, 2018) leading to a 2nd order convergence of the GW phase independent of the configuration details.

4 Quantitative Merger Dynamics

We start discussing our binary HyS mergers with a qualitative description of the coalescence of the two stars as depicted in Fig. 4. The first panel shows the initial setup for our evolution as computed from SGRID. The two stars have an initial coordinate separation of $\sim 61.5\rm km$, which corresponds to a proper distance of $\sim 76.4\rm km$. As visible in the figure, we use a gray scale for the SQM, while HM is shown on a blue to red scale. The second panel shows the system after approximately 6 orbits, i.e., 12 GW cycles. At this time one sees that the star’s surface gets ‘smeared out’ and is not perfectly preserved. Thus, lower density material surrounds the two stars. Similar observations are typical for full NR simulations and have been discussed in e.g. Faber and Rasio (2012) and references therein. The main reason for this, is the sharp transition of the density and the usage of an artificial atmosphere for the GRHD simulation, see e.g. Bernuzzi and Dietrich (2016). A similar observation is true about three to four orbit before the moment of merger (third panel). The fourth panel shows the system at the merger, i.e., when the GW amplitude peaks. One sees that the two SQM cores of the individual stars come into contact at this time. After the moment of merger, the central density inside the SQM core increases. The formed remnant (fifth panel of Fig. 4) survives for about $\sim 5$ milliseconds and undergoes two oscillations (three oscillations for lower resolutions) before it collapses to a black hole (BH). The final system (sixth panel) consists of a BH surrounded by an accretion disk composed of HM.

5 Gravitational Wave Signal

5.1 The Inspiral

We present the full gravitational waveform obtained from our simulations in Fig. 5 (top panel). Quantities related to the GW signal are plotted against the retarded coordinate time $u$, defined as

[TABLE]

where we choose $r_{\rm ext}=1000M_{\odot}$. The inspiral is characterized by an increasing frequency and amplitude of the GW, the so-called chirp-signal. We find that the merger time (marked by a vertical line) shifts by several milliseconds between the different resolutions. It is important to point out that for similar resolutions, cf. Dietrich et al. (2017a, b), the dephasing and consequently the shift in the merger time between similar resolutions have been significantly smaller for EoSs without phase transition. The bottom panel of Fig. 5 shows the phase differences between individual resolutions. One finds that for the three lowest resolutions, the phase is only converging with first order with respect to the grid spacing. Rescaling the phase differences assuming first order convergence leads to very good agreement with the phase difference of lower resolutions (dashed yellow line). This is of particular importance since it has been pointed out that the combination of the high-order flux reconstruction for the characteristic variables with the 5th order WENOZ scheme leads to a clear second-order convergence of the GW phase Bernuzzi and Dietrich (2016); Dietrich et al. (2017, 2018) for a variety of hadronic EoSs employing the same resolution and code settings. Increasing the resolution, we find that second-order convergence is recovered (dash-dotted green line). Ref. Bernuzzi and Dietrich (2016) pointed out that the NS’s surface causes this second-order convergence even though the numerical scheme is 5th order convergent for smooth hydrodynamical problems. Here, we find that in the case of a first order phase transition inside the NSs, the convergence order can even drop further to first order, but that with higher resolutions the convergence order increases. The reduced convergence order explains the large dephasing between the individual resolutions and displays the additional challenges in NR simulations of NS binaries employing a hybrid EoS.

Further work employing different hybrid EoSs and numerical schemes (varying the numerical flux scheme and limiters) is needed to understand if the observed feature is universal or might be caused by a particular choice of our parameters.

Finally, assuming second order convergence, we perform a Richardson extrapolation to obtain a better prediction for the GW phase, denoted by $\rm Re_{3,4}$. We compare $\rm Re_{3,4}$ with the GW approximant IMRPhenomPv2_NRTidal Dietrich et al. (2019, 2017); Hannam et al. (2014), which has been employed by the Ligo and Virgo Collaborations for the analysis of GW170817 in Abbott et al. (2019, 2018a, 2018b, 2018). This approximant (IMRPhenomPv2_NRTidal) is a phenomenological waveform model which has been developed for the description of BNS systems. During its construction only NR simulations with hadronic EOSs were used Dietrich et al. (2017). Therefore, we use this model to compare our results of hybrid star binaries with predictions for purely hadronic EOSs. The exact waveform IMRPhenomPv2_NRTidal, which we compare to, uses the same masses and tidal deformabilities as our NR simulation. Finally, the IMRPhenomPv2_NRTidal and our hybrid star waveform are aligned within a time window $u\in[5,10]\rm ms$ by minimizing the phase difference within this interval to allow a reasonable comparison and to compute the phase difference. The phase difference between resolution $\rm Re_{3,4}$ and IMRPhenomPv2_NRTidal is shown as a dash-dotted purple line in Fig. 5 (bottom panel). We find that until the moment of merger, the phase difference $\phi({\rm Re_{3,4}})-\phi({\rm IMRPhenomPv2\_NRTidal})$ is below the estimated error $\Delta\phi(\rm Re_{3,4},R_{4})$ 222Note that the larger dephasing at the beginning of the simulation can be explained by junk radiation., which is the first validation of a waveform model against a NR dataset including a strong phase transition inside the star. That the difference between the NR data and the waveform approximant is below the uncertainty of the data strengthens the reliability of models used for GW data analysis even if these have been derived for purely hadronic EOSs. However, due to the reduced convergence order and the large phase error of our simulations, further evolutions with higher resolutions are required for a more quantitative and stronger test.

5.2 The Postmerger

After the merger of the two stars, the rotating merger remnant emits GWs due to its asymmetric shape (see the fifth panel of Fig. 4). The typical emission during the postmerger happens on frequencies ranging from $2-4\rm kHz$. Which frequencies are excited and the strength of each individual frequency depends on the properties of the merger remnant and, thus, on the mass, spin, and EoS of the binary. We refer to existing studies available in the literature, e.g., Takami et al. (2014, 2015); Rezzolla and Takami (2016); Bauswein and Janka (2012); Stergioulas et al. (2011); Bauswein et al. (2012, 2014); Clark et al. (2016); Bauswein and Stergioulas (2015); Bernuzzi et al. (2015); Dietrich et al. (2017b, a, 2015), for further details and discussions about the mechanisms causing GW emission during the postmerger evolution. However, we want to understand if the presence of a SQM core changes the standard picture of the characteristic postmerger spectrum estimated for hadronic matter, see also Bauswein et al. (2019); Most et al. (2019).

For a better understanding of the frequency evolution, we present the spectrogram of the dominant (2,2)-mode of the GW signal in Fig. 6; we refer to Ref. Chaurasia et al. (2018) for a detailed discussion about the computation of the spectrogram. For our setup, the merger remnant collapses within $\sim 5\ \rm ms$ to a BH, due to this very short time the postmerger GW emission is not very strong. Nevertheless, one can extract the main emission frequency $f_{2}$ and obtains an emission of

[TABLE]

where the uncertainty is determined as the difference between the two highest available resolutions and also includes the uncertainty caused by the finite width of the peak in the frequency domain spectrum. We compare this $f_{2}$ frequency with the quasi-universal relation expressed in Ref. Bernuzzi et al. (2015) derived from pure hadronic EoSs. For a combined tidal deformability of $\tilde{\Lambda}=378$ 333The combined tidal deformability $\tilde{\Lambda}$ equals the individual tidal deformabilities $\Lambda^{A,B}$ presented in Tab. 1 for equal-mass binary stars. Ref. Bernuzzi et al. (2015) predicts $f_{2}=3.55\rm kHz$. We visualize this estimate as a horizontal dotted line in Fig. 6.

Overall, the difference of $0.3\rm kHz$ between our measurement and the quasi-universal relation is consistent with the uncertainty of the phenomenological fit Bernuzzi et al. (2015) combined with the uncertainty of our NR data. Thus, no clear imprint of the SQM core of the merger remnant on the postmerger evolution is present. This is opposite to the work of Bauswein et al. (2019) in which the merger remnant undergoes a first order phase transition which changes significantly the postmerger spectrum leading to a main emission frequency about $450\\ rmHz$ higher than the one predicted for a pure hadronic EoS. However, since the SQM imprint on the postmerger dynamics will depend on the exact EoS, our results are not in contradiction or tension to Bauswein et al. (2019). Similar to our results Ref. Most et al. (2019) finds a shorter merger remnant lifetime in the presence of a phase transition, but a similar postmerger frequency evolution as for pure hadronic matter.

6 Ejecta and electromagnetic counterparts

In addition to the emission of GWs, NS mergers are associated with a variety of EM signals such as short gamma-ray bursts (GRB), a synchrotron afterglow, and a thermal kilonova.

The thermal kilonova originates from neutron-rich outflows and the radioactive decay of $r$-process elements in the ejected matter, see e.g. Metzger (2017) and references therein. Kilonovae produce an almost isotropic EM emission and are visible in the ultraviolet/optical/infrared band. Current analyses of the kilonova AT2017gfo associated with GW170817 estimate the ejecta mass to be $\sim 5\times 10^{-2}M_{\odot}$. In our work, we will assume $\sim 5.4\times 10^{-2}M_{\odot}$ as found in Coughlin et al. (2018).

We present the amount of ejected matter present over the course of our simulation in Fig. 7 as a blue dashed line. We find that after the merger, the maximum density of the NS remnant increases and 3 oscillations are found. During this time most of the ejection happens and at the end of our simulation the ejecta mass is about $\sim 0.8\times 10^{-2}M_{\odot}$. This estimate agrees to about $30\%$ with resolution R3, which shows the overall consistency of our dataset. For comparison, the ejecta mass obtained in equal-mass BNS simulations employing SLy EoS with similar grid configurations and masses (Tab. III of Ref. Dietrich et al. (2017a), with $M^{A,B}=1.37M_{\odot}$ and Tab. V of Ref. Dietrich et al. (2015), with $M^{A,B}=1.35M_{\odot}$) are, respectively, $1.6\times 10^{-2}M_{\odot}$ and $1.22\times 10^{-2}M_{\odot}$. Consequently, our binary HyS configuration ejects only half of the amount of dynamical ejecta mass when compared to simulations employing the purely hadronic EOS SLy Dietrich et al. (2017a, 2015) 444We note that the phenomenological fit for hydonic EOSs for the dynamical ejecta mass proposed in Coughlin et al. (2018) (which is an updated version of the work presented in Dietrich and Ujevic (2017)) predicts a dynamical ejecta mass of $0.28\times 10^{-2}M_{\odot}$ , which lies even below our results. The discrepancy between the fit and the SLy simulations is caused by the large shock driven ejecta for the SLy EOS, which seems not to be fully captured in Coughlin et al. (2018). In addition, this comparison shows that further work with a larger set of HyS EOSs is required to obtain further confidence in the ejecta properties of NR simulations..

Due to the short evolution time after the merger, we are unable to estimate the total amount of ejecta due to the missing disk wind driven outflows which happen on longer timescales, see e.g. Kohri et al. (2005); Surman et al. (2006); Metzger et al. (2008); Dessart et al. (2009); Fernández and Metzger (2013); Perego et al. (2014); Siegel et al. (2014); Just et al. (2015); Rezzolla and Kumar (2015); Ciolfi and Siegel (2015); Siegel and Metzger (2017). However, we find that the disk mass drops to about $\sim 3.4\times 10^{-2}M_{\odot}$ towards the end of our simulation. Compared to the SLy simulations in Refs. Dietrich et al. (2017a, 2015) this is only about $\sim 15$ to $40\%$ of the SLy simulations. This means that if even half of the disk would be ejected (see Fig. 7), the kilonova associated with our studied HyS configuration would not be as bright as it was observed for GW170817.

Overall, the comparison of the ejecta and disk mass of our binary HyS and an SLy BNS suggests that HyS systems retain more mass within the stars along the time evolution and create dimmer EM counterparts. Further studies are necessary in order to understand whether this behavior is universal or restricted to our SQM EoS parameters choice.

To estimate the absolute magnitude of the kilonova connected to the HyS merger, we use the kilonova model of Kasen et al. (2017); Coughlin et al. (2018). The model of Kasen et al. (2017) employs a multi-dimensional Monte Carlo code to solve the multi-wavelength radiation transport equation for a relativistically expanding medium. To obtain the final predicted lightcurves, we interpolate between the existing Monte Carlo simulations using the Gaussian Process Regression techniques discussed in Coughlin et al. (2018). For our estimates, we use the dynamical ejecta component $M_{\rm ej,dyn}=8\times 10^{-3}M_{\odot}$ which leaves the system with an average velocity of $v_{\rm ej,dyn}=0.2$. Furthermore, we assume a lanthanide fraction $X_{\rm ej,dyn}=10^{-1}$. For the disk wind ejecta, we assume a velocity of $v_{\rm ej,wind}=0.1c$, a lanthanide fraction of $10^{-4}$, and an ejecta mass of $M_{\rm ej,wind}=0.017M_{\odot}$. The largest absolute magnitudes are obtained in the near-infrared, e.g. i- and K-band, while in the optical and ultraviolet the signal is fainter, see Fig. 8.

In contrast to the kilonova, short GRBs are created by highly relativistic outflows, powered by the merger remnant accretion disk, e.g., Paczynski (1986); Eichler et al. (1989). Based on the measured accretion disk mass of $\sim 3\times 10^{-2}M_{\odot}$ at the end of our simulation, we can obtain an order of magnitude estimate for the energy which could be released by a potential short GRB. Assuming that the GRB energy is proportional to the disk mass and that the proportionality constant is of about $10^{-2}$ as estimated in Coughlin et al. (2018), we find a GRB energy of $3\sim 10^{-4}M_{\odot}$, i.e., $6\times 10^{50}$ erg, which is consistent with current observational constraints for GRB170817A, e.g., van Eerten et al. (2018).

7 Conclusion

In this work, we presented and discussed the general procedure to construct HyS EoSs combining the tdBag model for SQM and a HM EoS for low-density material. Our procedure can be applied to study a variety of parameters and HM EoS choice and enables future studies of NSs binary containing SQM cores. The GW signal obtained from our NR simulations of a non-spinning, equal mass binary HyS merger exhibits a larger dephasing caused by numerical discretization than similar simulations for purely hadronic BNS configurations. We point out that the presence of a strong phase transition inside the stars can lead to a drop of the convergence to 1st order for insufficient resolutions even though higher order methods are applied. This is evidence that HyS mergers add additional complexity to the NR simulations. Further simulations must be done in order to understand if these features are universal or caused by the particular choice of EoS in connection with the employed numerical methods. Furthermore, we show that the postmerger GW main emission frequency for our constructed EoS is consistent with quasi-universal relations derived from hadronic EoSs within the uncertainties of the relation. Performing an analysis of the ejecta properties and EM signatures, the kilonova and the short GRB, we are able to draw the full multi-messenger picture for our HyS merger simulation. Our work is the first step towards a better coverage of the BNS parameter space and enables future investigations of a wider range of EoS with full NR simulations.

Appendix A HyS Piecewise Polytrope

In the following, we present the EOS parameters defining the employed HyS EOS; cf. Eq. (17) in section (2.4).