Binary Neutron Star Mergers: Mass Ejection, Electromagnetic Counterparts and Nucleosynthesis

Binary Neutron Star Mergers:
Mass Ejection, Electromagnetic Counterparts and Nucleosynthesis

David Radice11affiliation: Institute for Advanced Study, 1 Einstein Drive, Princeton, NJ 08540, USA 22affiliation: Department of Astrophysical Sciences, Princeton University, 4 Ivy Lane, Princeton, NJ 08544, USA , Albino Perego33affiliation: Istituto Nazionale di Fisica Nucleare, Sezione Milano Bicocca, gruppo collegato di Parma, I-43124 Parma, Italy 44affiliation: Dipartimento di Fisica, Università degli Studi di Milano Bicocca, Piazza della Scienza 3, 20126 Milano, Italy , Kenta Hotokezaka22affiliation: Department of Astrophysical Sciences, Princeton University, 4 Ivy Lane, Princeton, NJ 08544, USA , Steven A. Fromm55affiliation: NSCL/FRIB and Department of Physics & Astronomy, Michigan State University, 640 S Shaw Lane East Lansing, MI 48824, USA ,
Sebastiano Bernuzzi66affiliation: Theoretisch-Physikalisches Institut, Friedrich-Schiller-Universität Jena, 07743, Jena, Germany 33affiliation: Istituto Nazionale di Fisica Nucleare, Sezione Milano Bicocca, gruppo collegato di Parma, I-43124 Parma, Italy , and Luke F. Roberts55affiliation: NSCL/FRIB and Department of Physics & Astronomy, Michigan State University, 640 S Shaw Lane East Lansing, MI 48824, USA

We present a systematic numerical relativity study of the mass ejection and the associated electromagnetic transients and nucleosynthesis from binary neutron star (NS) mergers. We find that a few 10^{-3}\,M_{\odot} of material are ejected dynamically during the mergers. The amount and the properties of these outflow depend on binary parameters and on the NS equation of state (EOS). A small fraction of these ejecta, typically {\sim}10^{-6}\,M_{\odot}, is accelerated by shocks formed shortly after merger to velocities larger than 0.6\,{\rm c} and produces bright radio flares on timescales of weeks, months, or years after merger. Their observation could constrain the strength with which the NSs bounce after merger and, consequently, the EOS of matter at extreme densities. The dynamical ejecta robustly produce second and third r-process peak nuclei with relative isotopic abundances close to solar. The production of light r-process elements is instead sensitive to the binary mass ratio and the neutrino radiation treatment. Accretion disks of up to {\sim}0.2\,M_{\odot} are formed after merger, depending on the lifetime of the remnant. In most cases, neutrino- and viscously-driven winds from these disks dominate the overall outflow. Finally, we generate synthetic kilonova light curves and find that kilonovae depend on the merger outcome and could be used to constrain the NS EOS.

Subject headings:
Stars: neutron – Nuclear reactions, nucleosynthesis, abundances – Gamma-ray burst: general

1. Introduction

Merging neutron stars are loud gravitational wave (GW) sources and power bright electromagnetic (EM) transients, as recently demonstrated by GW170817 (Abbott et al., 2017d, 2018b, e, c; Arcavi et al., 2017; Coulter et al., 2017; Drout et al., 2017; Evans et al., 2017; Hallinan et al., 2017; Kasliwal et al., 2017; Nicholl et al., 2017; Smartt et al., 2017; Soares-Santos et al., 2017; Tanvir et al., 2017; Troja et al., 2017; Mooley et al., 2017; Ruan et al., 2018; Lyman et al., 2018). AT2017gfo, the EM counterpart to GW170817, had a UV/optical/infrared component with blackbody-like spectrum that has been interpreted as the result of the radioactive decay of about 0.05\ M_{\odot} of material ejected during and shortly after the merger (Chornock et al., 2017; Cowperthwaite et al., 2017; Drout et al., 2017; Nicholl et al., 2017; Tanaka et al., 2017; Tanvir et al., 2017; Perego et al., 2017a; Villar et al., 2017; Waxman et al., 2017; Metzger et al., 2018; Kawaguchi et al., 2018). Non thermal emissions from the radio to the \gamma-ray bands was also observed. The origin of this EM component are less clear, but recent observations point to the presence of synchrotron radiation generated by a relativistic jetted outflow interacting with the interstellar medium (ISM; Abbott et al. 2017c; Kasliwal et al. 2017; Margutti et al. 2017; Mooley et al. 2017; Lazzati et al. 2018; Margutti et al. 2018; Nakar & Piran 2018; Hotokezaka et al. 2018b; Alexander et al. 2018; Barkov et al. 2018; Mooley et al. 2018; Ghirlanda et al. 2018).

These groundbreaking observations are being used to constrain general relativity (GR; Abbott et al. 2017c; Baker et al. 2017; Pardo et al. 2018), cosmological parameters (Abbott et al., 2017a; Hotokezaka et al., 2018c), the equation of state (EOS) of NSs (Margalit & Metzger, 2017; Shibata et al., 2017a; Rezzolla et al., 2018; Ruiz et al., 2018; Bauswein et al., 2017; Annala et al., 2018; Radice et al., 2018c; Fattoyev et al., 2018; Paschalidis et al., 2018; Drago et al., 2018; Most et al., 2018b; De et al., 2018; Tews et al., 2018; Malik et al., 2018; Abbott et al., 2018a; Tsang et al., 2018; Hinderer et al., 2018; Wei et al., 2018; Nandi et al., 2018), and the origin of short \gamma-ray bursts (SGRBs; Abbott et al. 2017c; Mooley et al. 2017; Lazzati et al. 2018; Finstad et al. 2018; Beniamini et al. 2018), among others.

NS mergers also generate neutron-rich outflows that are thought to be a likely site for the production of r-process nuclei (Lattimer & Schramm, 1974; Symbalisty & Schramm, 1982; Meyer, 1989; Eichler et al., 1989; Freiburghaus et al., 1999; Goriely et al., 2011; Korobkin et al., 2012; Wanajo et al., 2014; Just et al., 2015; Thielemann et al., 2017). This is supported by the UV/optical/infrared follow up observations to GW170817 (e.g., Kasen et al., 2017; Rosswog et al., 2018; Hotokezaka et al., 2018a). The observed signals are thought to have been produced by the decay of radioactive isotopes produced during the r-process. However, uncertainty persists on the exact rates and nucleosynthetic yields from mergers. These needs to be addressed by future observations and theoretical studies.

Numerical simulations are the cornerstone to the modeling of multimessenger signatures and nucleosynthetic yields from NS mergers. Binary NS mergers are essentially multi-physics problems, involving strong-field gravity, strongly interacting matter, relativistic (magneto-)hydrodynamics, and neutrino transport. Due to the complexity of the problem, most previous binary NS merger simulations either employed an approximate treatment for the gravitational field of the NSs (Ruffert et al., 1996; Rosswog et al., 1999; Rosswog & Davies, 2003; Rosswog & Liebendoerfer, 2003; Rosswog et al., 2003; Oechslin et al., 2006; Rosswog et al., 2013; Korobkin et al., 2012; Bauswein et al., 2013), or included general relativistic (GR) effects, but compromised on the treatment of NS matter (Shibata & Uryu, 2000; Shibata et al., 2003, 2005; Baiotti et al., 2008; Kiuchi et al., 2010; Rezzolla et al., 2010, 2011; Hotokezaka et al., 2011, 2013b; Palenzuela et al., 2013; Hotokezaka et al., 2013a; Ruiz et al., 2016; Dietrich et al., 2017b, a). In the last few years, however, full-GR simulations with a microphysical treatment of NS matter and with different levels of approximation for neutrino-radiation effects have also become available (Sekiguchi et al., 2011, 2015; Palenzuela et al., 2015; Radice et al., 2016b; Lehner et al., 2016; Sekiguchi et al., 2016; Foucart et al., 2016b; Bovard et al., 2017).

Merger simulations have clarified the mechanisms driving the mass ejection and highlighted the importance of neutrino effects. Complementary studies on the long-term evolution of merger remnants also revealed the importance of the magnetically-, neutrino-, and viscously-driven secular outflows (Dessart et al., 2009; Metzger et al., 2008, 2009; Lee et al., 2009; Fernández & Metzger, 2013; Siegel et al., 2014; Just et al., 2015; Metzger & Fernández, 2014; Perego et al., 2014; Martin et al., 2015; Wu et al., 2016; Siegel & Metzger, 2017; Lippuner et al., 2017; Fujibayashi et al., 2017, 2018; Siegel & Metzger, 2018; Metzger et al., 2018; Fernández et al., 2018). These might in some cases dominate the nucleosynthetic yields and produce the bulk of the thermal radiation from NS mergers. However, quantitative questions remain concerning the dependency of the multimessenger emissions and of the r-process production on the binary parameters and the EOS. Indeed, most previous studies considered only several binary configurations, with the exception of Hotokezaka et al. (2013a) and Dietrich et al. (2017b, a), who however used an ideal-gas prescription to approximate thermal effects in the NS matter, and Bauswein et al. (2013) who used temperature dependent microphysical EOS, but adopted an approximate treatment of gravity and neglected weak reactions.

In this work, we present a systematic study of the mass ejection, nucleosynthetic yields, and EM counterparts of NS merger based on 59 high-resolution numerical relativity simulations. We employ four microphysical temperature-dependent nuclear EOSs and include the impact of neutrino losses. We consider total binary masses between 2.4\,M_{\odot} and 3.4\,M_{\odot} and mass ratios between 0.85 and 1. Our datasets also includes 13 simulations with an effective treatment of neutrino reabsorption and 6 simulations with viscosity. We quantify the dependency of dynamical ejecta mass and intrinsic properties on the binary parameters and NS EOS. We compute nucleosynthetic yields for the dynamical ejecta and show that second and third r-process peaks are robustly produced with relative isotopic abundances close to solar, while the production of light r-process elements is sensitive to the binary mass ratio and the treatment of neutrino radiation. We estimate kilonova light curves and show that their properties depend on the lifetime of the merger remnant. We also compute the expected radio signal from the interaction between the ejecta and the ISM. This signal could be used to probe the strength of the shocks generated after merger and, thus, indirectly the EOS of matter at extreme densities and temperatures. Our simulations also reveal a new outflow mechanism operating in unequal mass binaries and enabled by viscosity. These “viscous-dynamical” ejecta are launched due to the thermalization of mass exchange streams between the secondary and the primary NS shortly before merger and is discussed in more detail in a companion paper (Radice et al., 2018b).

The rest of this paper is organized as follows. Section 2 summarises the numerical methods and the initial data employed for the simulations. Section 3 discusses dynamical and secular mass ejection. Section 4 reports on our r-process nucleosynthesis calculations and yields. Section 5 is dedicated to the discussion of the EM counterparts from binary NS mergers, focusing on the kilonova signal and on the radio remnant powered by the interaction of the ejecta with the ISM.

Section 6 is dedicated to discussion and conclusions. Finally, Appendices A and B present a discussion of finite resolution effects and implementation details of the viscosity treatment employed by the simulations.

2. Methods

2.1. Initial Data

We construct initial data in quasi-circular orbit using the Lorene pseudo-spectral code (Gourgoulhon et al., 1999). The initial separation between the NSs is set to 40\,{\rm km}, corresponding to {\sim}2{-}3 orbits before merger. The EOSs used for the initial data are constructed from the minimum temperature slice of the EOS table used for the evolution assuming neutrino-less beta-equilibrium. To create the initial data tables we also subtract from the pressure the contribution of photon radiation, which dominates at the lowest densities due to the assumption of constant temperature. When reading the initial data in the evolution code we set the electron fraction according to the beta equilibrium condition, and we reset the specific internal energy according to the minimum temperature slice in the EOS table used for the evolution. Small errors in the initial data and in the mapping from the zero to the finite temperature EOSs induce small oscillations in the NSs prior to merger. We quantify these in terms of the relative change of the central density of the stars which we find typically to be \lesssim 2{-}3\%.

2.2. General-Relativistic Hydrodynamics

We evolve the initial data using the WhiskyTHC code (Radice & Rezzolla, 2012; Radice et al., 2014a, b, 2015).

WhiskyTHC separately evolves the proton and neutron number densities

\nabla_{\mu}(n_{p,n}u^{\mu})=R_{p,n}\,, (1)

where n_{p}=Y_{e}n is the proton number density, n_{n} is the neutron number density, n=n_{p}+n_{n} is the baryon number density, u^{\mu} the fluid four-velocity, and Y_{e} is the electron fraction of the material. R_{p}=-R_{n} is the net lepton number deposition rate due to the absorption and emission of neutrinos and anti-neutrinos (see Section 2.3).

We model NS matter as a perfect fluid with stress energy tensor

T_{\mu\nu}=(e+p)u_{\mu}u_{\nu}+pg_{\mu\nu}\,, (2)

where e is the energy density and p the pressure. We solve the Euler equations for the balance of energy and momentum

\nabla_{\nu}T^{\mu\nu}=Qu^{\mu}\,, (3)

where Q is the net energy deposition rate due to the absorption and emission of neutrinos (see Section 2.3).

WhiskyTHC discretizes Eqs. (1) and (3) using high-resolution shock-capturing (HRSC) schemes. For the simulations presented in this work we use a central Kurganov-Tadmor type scheme (Kurganov & Tadmor, 2000) employing the HLLE flux formula (Einfeldt, 1988) and non-oscillatory reconstruction of the primitive variables with the MP5 scheme of Suresh (1997). For numerical reasons we embed the NSs in a low density medium, with \rho_{0}=m_{b}\,n\simeq 6\times 10^{4}\,{\rm g}\,{\rm cm}^{-3}, where m_{b} is the atomic mass unit. We use the best available numerical schemes for the treatment of low density regions and for the advection of the fluid composition. We employ the positivity-preserving limiter from Radice et al. (2014b) to ensure rest-mass conservation even in the presence of the artificial density floor, and we use the consistent multi-fluid advection method of Plewa & Müller (1999) to ensure separate local conservation of the proton and neutron number densities. Furthermore, we extract the outflows properties when the density is still several orders of magnitude higher than that of the artificial atmosphere.

The spacetime is evolved using the Z4c formulation of Einstein’s equations (Bernuzzi & Hilditch, 2010; Hilditch et al., 2013) as implemented in the CTGamma code (Pollney et al., 2011; Reisswig et al., 2013b), which is part of the Einstein Toolkit (Löffler et al., 2012). CTGamma implements fourth-order finite-differencing of the equations and we use fifth-order Kreiss-Oliger dissipation to ensure the nonlinear stability of the evolution (Kreiss & Oliger, 1973). The coupling between the hydrodynamics and the spacetime evolution is handled using the method of lines (MoL). We adopt the optimal strongly-stability preserving third-order Runge-Kutta scheme (Gottlieb et al., 2008) as time integrator. The timestep is set according to the speed-of-light Courant-Friedrichs-Lewy (CFL) condition with CFL factor 0.15.

Our simulations domain covers a cube of 3,024 km in diameter whose center is at the center of mass of the binary. Our code uses Berger-Oliger conservative adaptive mesh refinement (AMR; Berger & Oliger 1984) with sub-cycling in time and refluxing (Berger & Colella, 1989; Reisswig et al., 2013a) as provided by the Carpet module of the Einstein Toolkit (Schnetter et al., 2004). We setup an adaptive mesh-refinement (AMR) grid structure with 7 refinement levels. The finest refinement level covers both NSs during the inspiral and the remnant after the merger and has a typical resolution of h\simeq 185\,{\rm m}. For selected binaries, we also perform simulations with finest grid resolution of 123\,{\rm m} and 246\,{\rm m}. See Table 2 for a summary of our models.

2.3. Neutrino Leakage Scheme

We treat compositional and energy changes in the material due to weak reactions using a leakage scheme (Galeazzi et al. 2013; Radice et al. 2016b; see also van Riper & Lattimer 1981; Ruffert et al. 1996; Rosswog & Liebendoerfer 2003; O’Connor & Ott 2010; Sekiguchi 2010; Neilsen et al. 2014; Perego et al. 2016; Ardevol-Pulpillo et al. 2018 for other implementations). Our scheme tracks reactions involving electron \nu_{e} and anti-electron type \bar{\nu}_{e} neutrinos separately. Heavy-lepton neutrinos are lumped together in a single effective specie labeled \nu_{x}. We account for the reactions listed in Table 1. We use the formulas from the references listed in the tables to compute the neutrino production rates R_{\nu}, \nu\in\{\nu_{e},\bar{\nu}_{e},\nu_{x}\}, the associated energy release Q_{\nu}, and neutrino absorption \kappa_{\nu,a} and scattering \kappa_{\nu,s} opacities. In doing so, we assume that the neutrinos follow Fermi-Dirac distributions with chemical potentials obtained assuming beta-equilibrium with thermalized neutrinos as in Rosswog & Liebendoerfer (2003). Following Ruffert et al. (1996), we also distinguish between the number density weighted opacities \kappa_{\nu,a}^{0} and \kappa_{\nu,s}^{0} that determine the rate at which neutrinos diffuse out of the material, and the energy density weighted opacities \kappa_{\nu,a}^{1} and \kappa_{\nu,s}^{1} that determine the rate at which energy diffuses out of the material due to the loss of neutrinos.

Table 1Weak reaction rates and references for their implementation. We use the following notation \nu\in\{\nu_{e},\bar{\nu}_{e},\nu_{x}\} denotes a neutrino, \nu_{x} denote any heavy-lepton neutrino, N\in\{n,p\} denotes a nucleon, and A denotes a nucleus.
Reaction Reference
\nu_{e}+n\leftrightarrow p+e^{-} Bruenn (1985)
\bar{\nu}_{e}+p\leftrightarrow n+e^{+} Bruenn (1985)
e^{+}+e^{-}\rightarrow\nu+\bar{\nu} Ruffert et al. (1996)
\gamma+\gamma\rightarrow\nu+\bar{\nu} Ruffert et al. (1996)
N+N\rightarrow\nu+\bar{\nu}+N+N Burrows et al. (2006)
\nu+N\rightarrow\nu+N Ruffert et al. (1996)
\nu+A\rightarrow\nu+A Shapiro & Teukolsky (1983)

The total neutrino opacities \kappa_{\nu,a}^{\alpha}+\kappa_{\nu,s}^{\alpha}, with \alpha\in\{0,1\}, are used to compute an estimate to the optical depth \tau^{\alpha}_{\nu} following the scheme of Neilsen et al. (2014), which is well-suited to the complex geometries present in NS mergers. We use the optical depth to define effective emission rates as Ruffert et al. (1996):

R_{\nu}^{\rm eff}=\frac{R_{\nu}}{1+t_{\rm diff}^{0}\,(t_{\rm loss}^{0})^{-1}}\,, (4)

where we have introduced the effective diffusion time t_{\rm diff}

t_{\rm diff}^{0}=\mathcal{D}\frac{(\tau^{0}_{\nu})^{2}}{\kappa_{\nu,a}^{0}+% \kappa_{\nu,s}^{0}}\,, (5)

the neutrino emission timescale

t_{\rm loss}^{0}=\frac{R_{\nu}}{n_{\nu}}\,, (6)

and n_{\nu} is the neutrino number density computed assuming beta-equilibrium with neutrinos. The constant \mathcal{D} is a tuning parameter that we set to 6. The effective energy emission rates Q_{\nu}^{\rm eff} are computed along the same lines, but using \tau^{1}_{\nu} and \kappa^{1}_{\nu,a}, \kappa^{1}_{\nu,s} instead of \tau^{0}_{\nu}, \kappa^{0}_{\nu,a}, and \kappa^{0}_{\nu,s}, respectively.

Neutrinos are split into a trapped component n_{\nu}^{\rm trap} and a free-streaming component n_{\nu}^{\rm fs}. Free-streaming neutrinos are emitted according to the effective rate R_{\nu}^{\rm eff} of Eq. (4) and with average energy Q_{\nu}^{\rm eff}/R_{\nu}^{\rm eff} and then evolved according to the M0 scheme we introduced in Radice et al. (2016b) and that is briefly summarized below. In our implementation the pressure due to the trapped neutrino component is neglected, since it is found to be unimportant in the conditions relevant for NS mergers (Galeazzi et al., 2013).

The M0 scheme evolves the number density of the free-streaming neutrinos assuming that they move along radial null rays with four vector k^{\alpha} normalized so that k^{\alpha}u_{\alpha}=-1. Under these assumptions it is possible to show that that the free-streaming fluid rest frame neutrino number density n_{\nu}^{\rm fs} satisfies (Radice et al., 2016b)

\nabla_{\alpha}[n_{\nu}^{\rm fs}k^{\alpha}]=R_{\nu}^{\mathrm{eff}}-\kappa_{\nu% ,a}^{\rm eff}n_{\nu}^{\rm fs}\,, (7)

where R_{\nu}^{\rm eff} is the effective luminosity from Eq. (4) and the effective absorption rates are defined as

\kappa_{\nu,a}^{\rm eff}=e^{-\tau_{\nu}^{0}}\,\left(\frac{E_{\nu}^{\rm fs}}{E_% {\nu}^{\beta}}\right)^{2}\,\kappa_{\nu,a}^{0}\,. (8)

We have also introduced the average energy of free-streaming neutrinos E_{\nu}^{\rm fs} and the average neutrino energy in beta-equilibrium E_{\nu}^{\beta}, both defined in the fluid rest frame. Note that in the simulations we reported in Radice et al. (2016b) the term in parenthesis on the right hand side of Eq. (8) was neglected. We report a comparison of results obtained with and without the inclusion of this term in Sec. 4.1.

Our scheme estimates the free-streaming neutrino energy under the additional assumption, only approximately satisfied in our simulations, of stationarity of the metric. Accordingly, \partial_{t} is assumed to be a Killing vector so that p^{\alpha}_{\nu}(\partial_{t})_{\alpha}, p^{\alpha}_{\nu} being the neutrinos four-momentum, is a conserved quantity. Under this assumption it is possible to show that the average free-streaming neutrino energy satisfies (Radice et al., 2016b)

k^{t}\partial_{t}(E_{\nu}^{\rm fs}\chi)+k^{r}\partial_{r}(E_{\nu}^{\rm fs}\chi% )=\frac{\chi}{n_{\nu}^{\rm fs}}\left(Q_{\nu}^{\rm eff}-E_{\nu}^{\rm fs}R_{\nu}% ^{\rm eff}\right)\,, (9)

where \chi=-k^{\alpha}(\partial_{t})_{\alpha}.

In the simulations in which it is employed the M0 scheme is switched on shortly before the two NSs collide, when neutrino matter interactions start to become dynamically important. We solve Eqs. (7) and (9) on a uniform spherical grid extending to 512\,G/c^{2}\,M_{\odot}\simeq 756\,{\rm km} and having n_{r}\times n_{\theta}\times n_{\phi}=3096\times 32\times 64 grid points.

The coupling between neutrinos and matter is treated using an operator split approach discussed in detail in Radice et al. (2016b). The source terms for Eqs. (1) and (3) are, respectively,

R_{p}=({\kappa}_{{\nu}_{e},a}^{\rm eff}n_{\nu_{e}}^{\rm fs}-{\kappa}_{\bar{\nu% }_{e},a}^{\rm eff}n_{\bar{\nu}_{e}}^{\rm fs})-(R^{\mathrm{eff}}_{\nu_{e}}-R^{% \mathrm{eff}}_{\bar{\nu}_{e}})\,, (10)


\begin{split}\displaystyle Q=({\kappa}_{\nu_{e},a}^{\rm eff}n_{\nu_{e}}^{\rm fs% }E_{\nu_{e}}&\displaystyle+{\kappa}_{\bar{\nu}_{e},a}^{\rm eff}n_{\bar{\nu}_{e% }}^{\rm fs}E_{\bar{\nu}_{e}})\\ &\displaystyle-(Q^{\mathrm{eff}}_{\nu_{e}}+Q^{\mathrm{eff}}_{\bar{\nu}_{e}}+Q^% {\mathrm{eff}}_{\nu_{x}})\,.\end{split} (11)

The M0 scheme is more approximate than the frequency-integrated M1 scheme adopted by, e.g., Sekiguchi et al. (2015) and Foucart et al. (2015). However it has the advantage of computational efficiency, it includes gravitational and Doppler effects, albeit in an approximate way, and does not suffer from unphysical radiation shocks in the important region above the merger remnant that instead plagues M1 schemes (Foucart et al., 2018a).

2.4. Viscosity

For a subset of binaries we use the general-relativistic large eddy simulations method (GRLES; Radice 2017) to explore the impact of subgrid-scale turbulent angular momentum transport. Accordingly, we decompose the stress energy tensor of the fluid T^{\mu\nu} as

\mathchoice{T_{\kern 0.0pt\hbox{$\scriptstyle\mu\nu$}}\kern-10.49984pt\kern 0.% 0pt\kern 0.0pt{}^{\hbox{$\scriptstyle{}$}\kern 0.0pt\hbox{$\kern 0.0pt% \scriptstyle\hbox{}\hbox{}$}}}{T_{\kern 0.0pt\hbox{$\scriptstyle\mu\nu$}}\kern% -10.49984pt\kern 0.0pt\kern 0.0pt{}^{\raise 0.43pt\hbox{$\scriptstyle{}$}\kern 0% .0pt\raise 0.43pt\hbox{$\kern 0.0pt\scriptstyle\hbox{}\hbox{}$}}}{T_{\kern 0.0% pt\hbox{$\scriptscriptstyle\mu\nu$}}\kern-7.499886pt\kern 0.0pt\kern 0.0pt{}^{% \raise 0.43pt\hbox{$\scriptscriptstyle{}$}\kern 0.0pt\raise 0.43pt\hbox{$\kern 0% .0pt\scriptscriptstyle\hbox{}\hbox{}$}}}{T_{\kern 0.0pt\hbox{$% \scriptscriptstyle\mu\nu$}}\kern-7.499886pt\kern 0.0pt\kern 0.0pt{}^{\raise 0.% 43pt\hbox{$\scriptscriptstyle{}$}\kern 0.0pt\raise 0.43pt\hbox{$\kern 0.0pt% \scriptscriptstyle\hbox{}\hbox{}$}}}=En_{\mu}n_{\nu}+S_{\mu}n_{\nu}+S_{\nu}n_{% \mu}+S_{\mu\nu}\,, (12)


\displaystyle E=T_{\mu\nu}n^{\mu}n^{\nu}=\rho hW^{2}-p\,, (13)
\displaystyle S_{\mu}=-\gamma_{\mu\alpha}n_{\beta}T^{\alpha\beta}=\rho hW^{2}v% _{\mu}\,, (14)
\displaystyle S_{\mu\nu}=\gamma_{\mu\alpha}\gamma_{\mu\beta}T^{\alpha\beta}=S_% {\mu}v_{\nu}+p\gamma_{\mu\nu}+\tau_{\mu\nu}\,, (15)

and n^{\mu} is the normal to the space-like slice hyper-surface, while \gamma_{\mu\nu}, v^{\mu}, and W are, respectively, the spatial metric, the three-velocity, and the Lorentz factor. \tau_{\mu\nu} is a purely spatial tensor representing the effect of subgrid scale turbulence. As in Radice (2017), we use the following ansatz for \tau_{ij}:

\tau_{ij}=-2\nu_{{}_{T}}(\rho+p)W^{2}\left[\frac{1}{2}\big{(}D_{i}v_{j}+D_{j}v% _{i}\big{)}-\frac{1}{3}D_{k}v^{k}\gamma_{ij}\right]\,, (16)

where \nu_{{}_{T}}=\ell_{\rm mix}c_{s} is the turbulent viscosity, c_{s} is the sound speed, D_{i} denotes the covariant derivative compatible with the spatial metric, and \ell_{\rm mix} is a free parameter we vary to study the sensitivity of our results to turbulence. In the context of accretion disk theory turbulent viscosity is typically parametrized in terms of a dimensionless constant \alpha linked to \ell_{\rm mix} through the relation \ell_{\rm mix}=\alpha\,c_{s}\,\Omega^{-1}, where \Omega is the angular velocity of the fluid (Shakura & Sunyaev, 1973). Recently, Kiuchi et al. (2018) performed very high resolution general-relativistic magnetohydrodynamics (GRMHD) simulations of a NS merger with sufficiently high seed magnetic fields (10^{15}\,{\rm G}) to be able to resolve the magnetorotational instability (MRI) in the merger remnant and reported averaged \alpha values for different rest-mass density shells. Combining their estimate of \alpha with values of c_{s} and \Omega from our simulations we find values of \ell_{\rm mix}=0{-}30\,{\rm m}. Here, we conservatively vary \ell_{\rm mix} between 0 (default; no subgrid model) and 50\,{\rm m} (very efficient angular momentum transport).

WhiskyTHC consistently includes the contributions of \tau_{\mu\nu} to the fluid stress energy tensor in the calculation of the right hand side of the metric and fluid equations. Flux terms resulting from the inclusion of \tau_{\mu\nu} need to be treated with care, because a naive application of Godunov-type methods would result in a scheme suffering from an odd-even decoupling instability (e.g., Lowrie & Morel, 2001). WhiskyTHC uses a proper combination of left and right biased finite-differencing operators to discretize terms arising from the derivatives of \tau_{\mu\nu} in a flux-conservative fashion. The details are given in Appendix B.

2.5. Models

Figure 1.— NS masses and radii predicted by the EOSs considered in this study. The BHB\Lambda\phi and DD2 EOS predict the same radii for NSs less massive than {\sim}1.5\,M_{\odot}. BHB\Lambda\phi predicts smaller radii for more massive NSs and a smaller maximum mass. The LS220 and SFHo EOSs have similar maximum masses and similar compactness close to the maximum mass. However, LS220 predicts radii almost 1\,{\rm km} larger than SFHo for a canonical 1.4\,M_{\odot} NS.

For our survey we consider four different nuclear EOS: the DD2 EOS (Typel et al., 2010; Hempel & Schaffner-Bielich, 2010), the BHB\Lambda\phi EOS (Banik et al., 2014), the LS220 EOS (Lattimer & Swesty, 1991), and the SFHo EOS (Steiner et al., 2013). These EOSs predict NS maximum masses and radii within the range allowed by current astrophysical constraints, including the recent LIGO/Virgo constraint on tidal deformability (Abbott et al., 2017d, 2018b; De et al., 2018; Abbott et al., 2018a). The LS220 EOS is based on a liquid droplet Skyrme model while the other three EOSs are based on nuclear statistical equilibrium with a finite volume correction coupled to a relativistic mean field theory for treating high-density nuclear matter. DD2 and SFHo use different parameterizations and values for modeling the mean-field nuclear interactions. The BHB\Lambda\phi EOS uses the same nucleon interactions as the DD2 EOS, but also includes interacting \Lambda hyperons that can be produced at high density and soften the EOS.

Although these EOSs differ in many aspects, including their finite temperature properties and their dependence on the neutron richness of the system, we can generally characterize them by the TOV solutions they predict, which we show in Fig. 1. SFHo, LS220, DD2, and BHB\Lambda\phi support 2.06, 2.06, 2.42, and 2.11 M_{\odot} cold, non-rotating maximum NS masses and have R_{1.4} of 11.9, 12.7, 13.2, and 13.2 km, respectively. Since NS radii correlate with the pressure at roughly twice saturation density (Lattimer, 2012), we refer to EOSs having smaller R_{1.4} as being “softer” and to EOSs having larger R_{1.4} as being stiffer. Although all four models have saturation density symmetry energies within experimental bounds, LS220 has a significantly steeper density dependence of its symmetry energy than the other models. In all models, the finite temperature behavior of the EOS is mainly determined by the nucleon effective mass, with smaller effective masses leading to higher temperatures for constant entropy. The LS220 EOS assumes that the nucleon mass is the bare nucleon mass at all densities, while SFHo has m_{N}^{*}/m_{N}=0.76 at saturation density, and both DD2 and BHB\Lambda\phi have m_{N}^{*}/m_{N}=0.56, where m_{N}^{*} is the effective nucleon mass and m_{N} is the bare nucleon mass.

We considered 35 distinct binaries with total masses M\in[2.4,3.4]\,M_{\odot} and mass ratios in the range q\in[0.85,1]. All binaries have been simulated using the leakage scheme discussed above, selected binaries have also been simulated with the M0 scheme, at different resolutions, and/or with viscosity, for a total of 59 simulations. A summary of all simulations and key results is given in Table 2. Each run is labelled after the employed EOS, the masses of the two NSs at infinite separation, the simulated physics, the value of the mixing length \ell_{\rm mix} in meters, if larger than zero, and, in the case of binaries run at multiple resolutions, the resolution (LR: low resolution, HR: high resolution). For example, LS220_M140120_M0_L25 is a binary with NSs masses 1.4\,M_{\odot} and 1.2\,M_{\odot} that was simulated with the LS220 EOS, employing the M0 scheme, with \ell_{\rm mix}=25\,{\rm m}, and run at our standard resolution h=185\,{\rm m}. We will make all of our GW waveforms publicly available as part of the CoRe catalog (Dietrich et al., 2018).

We simulate all binaries for at least 20\,{\rm ms} after the merger, or until few milliseconds after black hole (BH) formation if this occurs earlier. With the exceptions of the runs including viscosity, which we discuss in more details in a companion paper (Radice et al., 2018b), the outflow rate, precisely defined in Sec 3, has dropped to zero at the end of our simulations. Our models include binaries spanning all range of possible outcomes predicted for NS binary mergers (e.g., Shibata, 2016). Some of our binaries produce BHs promptly at the time of the merger, others produce hypermassive neutron stars (HMNSs) that collapse on time scales of several milliseconds (Baumgarte et al., 2000; Shibata & Taniguchi, 2006; Baiotti et al., 2008; Sekiguchi et al., 2011; Bernuzzi et al., 2016), or long-lived massive NS remnants, expected to be stable on secular timescales or, in some cases, indefinitely (Giacomazzo & Perna, 2013; Foucart et al., 2016a; Radice et al., 2018a). Table 2 reports the time to BH formation in milliseconds from the merger t_{\rm BH}, or a lower limit if the central object does not collapse within the simulation time.

3. Mass Ejection

Tidal interactions and shocks exerted on the NSs close to the time of merger trigger the ejection of material on a dynamical time scale, the so called dynamical ejecta (e.g., Hotokezaka et al., 2013b; Bauswein et al., 2013; Radice et al., 2016b). Dynamical ejecta mass and average properties are reported in Table 2. In simulations that do not include viscosity the outflow rate drops to zero few milliseconds after the merger. However, it is expected that more material will become unbound from the remnant on longer timescales, due to magnetic effects and/or nuclear recombination (e.g., Fernández & Metzger, 2013; Perego et al., 2014; Siegel & Metzger, 2017; Fujibayashi et al., 2018; Fernández et al., 2018), which we cannot presently study with our simulations. We refer to this latter component as the secular ejecta to distinguish it from the dynamical ejecta defined above. The simulations performed with viscosity show the early development of viscously-driven outflows. However, due to the high computational costs, we do not follow the postmerger remnant for sufficiently long times to study the secular ejecta.

3.1. Dynamical Ejecta

Figure 2.— Volume rendering of the electron fraction of the ejecta for the simulation SFHo_M135135_M0. The ray-casting opacity is linear in the logarithm of the rest-mass density. From the top-left in clockwise direction, the transparency minimum – maximum in the opacity scale are (10^{11}-10^{14}) {\rm g}\ {\rm cm}^{-3}, (10^{8}-10^{11}) {\rm g}\ {\rm cm}^{-3}, (10^{8}-10^{11}) {\rm g}\ {\rm cm}^{-3}, and (10^{7}-10^{11}) {\rm g}\ {\rm cm}^{-3}. The last panel of this figure should be compared with Fig. 14 where we plot a cut of the data in the xz-plane.

We discuss the qualitative properties of the dynamical ejecta in the case of the SFHo_M135135_M0 binary, which we take as fiducial. Figure 2 summarizes its most salient features. The bulk of the dynamical outflow is contained within a wide {\sim}60^{\circ} angle from the orbital plane. In the simulations that do not account for neutrino absorption the outflow is neutron rich, with average electron fraction \langle Y_{e}\rangle<0.2. When neutrino absorption is included in the simulations, as is the case for the SFHo_M135135_M0 binary, the ejecta is reprocessed to higher values of Y_{e}, but remains neutron rich \langle Y_{e}\rangle<0.25.

As shown in Fig. 2, the tidal component of the dynamical ejecta is emitted first close to the orbital plane, followed by a more isotropic shock heated ejecta component. However, the shock heated component generally has higher velocity, so it rapidly overtakes the tidal component. The two components interact and, as a consequence, the tidal tail is partially reprocessed by weak reactions to slightly higher values of Y_{e}\simeq 0.1. Later, we observe the emergence of a neutrino driven wind component of the outflow, concentrated at high latitudes. Finally, for this binary, mass ejection shuts off shortly after the formation of a BH.

Figure 3.— Maximum rest-mass density and outflow rate for the SFHo_M135135_M0 binary. The outflow rate is computed at a radius of R=443\,{\rm km} and shifted in time by R\,(0.5c)^{-1}, 0.5c being roughly the 99 percentile of the velocity of the ejecta at the radius R. We show both the total outflow rate as well as the outflow rate of only the material with asymptotic velocity larger than 0.6\,c. We find that the bulk of the ejecta is launched when the HMNS first bounces back after the merger.

As the two NSs merge, their inner cores are violently compressed against each other. For some of our models this compression is sufficient to trigger a runaway collapse of the remnant and a BH forms within a single dynamical timescale. However, for most of our binaries the angular momentum of the remnant is sufficiently large to prevent the collapse and the remnant’s core undergoes a centrifugal bounce. The bounce starts a cycle of large scale oscillations of the remnant. As also discussed in detail by Bauswein et al. (2013), for most binaries the most abundant component of the outflow is triggered during the first expansion of the remnant. This is demonstrated in Fig. 3, where we show the maximum density and the outflow rates measured for our fiducial SFHo_M135135_M0 binary. The former is extracted from the flux of unbound material leaving a coordinate sphere with radius R=300\,G/c^{2}\,M_{\odot}\simeq 443\,{\rm km} and is retarded according to the velocity of the tip of the ejecta at the detection sphere. As evidenced in the figure, and confirmed by a careful inspection of the multidimensional data, the bulk of the outflow is triggered at the time when the merger remnant rebounds. However, after the first bounce, the subsequent oscillations do not produce significant mass ejection. The outflow rate remains positive for a few milliseconds as slower material, that has also started expanding during the first ejection episode, reaches the detector.

Figure 4.— Angular distribution and composition of the ejecta for the 1.35\ M_{\odot} vs 1.35\ M_{\odot} binary with BHB\Lambda\phi EOS. Right panel: neutrino cooling only. Left panel: simulation with neutrino absorption. Neutrino irradiation is necessary to generate high-Y_{e}, polar outflows.

The inclusion of neutrino re-absorption results in the formation of a new outflow component. This additional outflow stream is composed of material that is ablated from the surface of the hypermassive neutron star (HMNS) and its accretion disk. This material is reprocessed to Y_{e}>0.25 and channeled along the polar direction \theta<45^{\circ}. We remark that this outflow component is absent if neutrino heating is not included (Radice et al. 2016b; Perego et al. 2017a; see also Fig. 4 for another representative case). Furthermore, because in our simulation the tidal streams interact with the shocked component of the ejecta, the former also see their Y_{e} slightly increased with the inclusion of neutrino absorption, as can be seen from the shift in Y_{e} of the material close to the orbital plane (see region \theta\simeq 90^{\circ} in Fig. 4).

3.1.1 Outflow Properties

Table 2Dynamical ejecta and remnant disks for all simulations. We report the model name, the grid resolution h, the NS masses at infinite separation M_{a} and M_{b}, the time of BH formation t_{\rm BH}, the remnant disk masses M_{\rm disk}, the total ejected mass M_{\mathrm{ej}}, and the amount of fast moving ejecta M_{\rm ej}^{v\geq 0.6c}. For the dynamical ejecta we report the mass-averaged electron fraction \langle Y_{e}\rangle, specific entropy per baryon \langle s\rangle, and asymptotic velocity v_{\rm ej}. We also give the rms opening angle of the outflow streams from the orbital plane \theta_{\rm ej}=\sqrt{\langle\theta^{2}\rangle}, and the total kinetic energy of the ejecta T_{\rm ej}. All ejecta properties are measured on a sphere with coordinate radius of 300\,G/c^{2}\,M_{\odot}\simeq 443\ \mathrm{km}.
Model h M_{a} M_{b} t_{\rm BH} M_{\rm disk} M_{\rm ej} M_{\rm ej}^{v\geq 0.6c} \langle Y_{e}\rangle \langle s\rangle v_{\rm ej} \theta_{\rm ej} T_{\rm ej}
[{\rm m}] [M_{\odot}] [M_{\odot}] [{\rm ms}] [10^{-2}\ M_{\odot}] [10^{-5}\ M_{\odot}] [k_{\rm B}] [c] [10^{50}\ {\rm erg}]
BHBlp_M125125_LK 185 1.25 1.25 >\ 21.9 NA 0.13 0.000 0.13 15 0.18 22 0.46
BHBlp_M1365125_LK 185 1.365 1.25 >\ 24.0 18.73 0.06 1.367 0.14 17 0.21 25 0.35
BHBlp_M130130_LK 185 1.3 1.3 >\ 26.9 NA 0.07 0.000 0.16 22 0.12 33 0.11
BHBlp_M135135_LK 185 1.35 1.35 >\ 21.3 14.45 0.07 0.746 0.15 20 0.17 28 0.26
BHBlp_M135135_LK_HR 123 1.35 1.35 >\ 23.3 14.26 0.05 0.015 0.16 20 0.18 24 0.19
BHBlp_M135135_M0 185 1.35 1.35 >\ 37.0 12.75 0.14 0.283 0.26 24 0.14 38 0.40
BHBlp_M140120_LK 185 1.4 1.2 >\ 23.7 20.74 0.11 0.229 0.11 13 0.16 22 0.36
BHBlp_M140120_M0 185 1.4 1.2 >\ 28.2 22.56 0.16 0.994 0.19 17 0.17 30 0.60
BHBlp_M140140_LK 185 1.4 1.4 12.0 7.05 0.09 0.232 0.15 18 0.17 28 0.34
BHBlp_M140140_LK_HR 123 1.4 1.4 10.3 5.38 0.10 0.793 0.14 17 0.20 26 0.50
BHBlp_M144139_LK 185 1.44 1.39 10.4 8.28 0.06 0.511 0.18 22 0.20 30 0.30
BHBlp_M150150_LK 185 1.5 1.5 02.3 1.93 0.05 0.727 0.17 20 0.23 28 0.33
BHBlp_M160160_LK 185 1.6 1.6 01.0 0.09 0.00 0.000 - - - - -
DD2_M120120_LK 185 1.2 1.2 >\ 24.7 NA 0.08 0.000 0.15 21 0.14 31 0.16
DD2_M125125_LK 185 1.25 1.25 >\ 32.4 NA 0.04 0.000 0.18 27 0.15 33 0.10
DD2_M1365125_LK 185 1.365 1.25 >\ 24.2 20.83 0.04 0.443 0.15 21 0.20 25 0.20
DD2_M130130_LK 185 1.3 1.3 >\ 22.9 NA 0.12 0.005 0.13 15 0.18 21 0.45
DD2_M135135_LK 185 1.35 1.35 >\ 24.4 15.69 0.03 0.024 0.18 27 0.18 31 0.12
DD2_M135135_LK_HR 123 1.35 1.35 >\ 23.4 15.05 0.02 0.001 0.18 28 0.19 30 0.09
DD2_M135135_M0 185 1.35 1.35 >\ 20.4 16.16 0.14 0.168 0.23 21 0.17 31 0.49
DD2_M140120_LK 185 1.4 1.2 >\ 23.6 19.26 0.09 0.307 0.12 15 0.18 23 0.36
DD2_M140120_M0 185 1.4 1.2 >\ 26.9 19.48 0.16 0.570 0.21 18 0.16 34 0.54
DD2_M140140_LK 185 1.4 1.4 >\ 24.5 12.36 0.04 0.529 0.17 22 0.22 28 0.26
DD2_M140140_LK_HR 123 1.4 1.4 >\ 24.6 16.85 0.09 0.431 0.14 17 0.18 28 0.39
DD2_M144139_LK 185 1.44 1.39 >\ 23.5 14.40 0.05 0.158 0.17 22 0.20 28 0.26
DD2_M150150_LK 185 1.5 1.5 >\ 23.1 16.70 0.07 0.462 0.20 23 0.17 32 0.26
DD2_M160160_LK 185 1.6 1.6 02.3 1.96 0.12 2.759 0.14 13 0.24 20 0.83
LS220_M120120_LK 185 1.2 1.2 >\ 23.2 17.43 0.14 0.000 0.12 15 0.15 25 0.33
LS220_M1365125_LK 185 1.365 1.25 >\ 26.7 16.86 0.11 0.013 0.10 13 0.16 23 0.32
LS220_M135135_LK 185 1.35 1.35 20.3 7.25 0.06 0.000 0.11 16 0.16 23 0.17
LS220_M135135_LK_LR 246 1.35 1.35 >\ 27.6 NA 0.11 0.000 0.13 16 0.16 26 0.32
LS220_M135135_M0 185 1.35 1.35 >\ 22.6 9.06 0.19 0.009 0.25 19 0.15 35 0.50
LS220_M135135_M0_L05 185 1.35 1.35 >\ 24.5 14.07 0.27 0.061 0.27 20 0.13 38 0.58
LS220_M135135_M0_L25 185 1.35 1.35 18.1 4.65 0.20 0.395 0.26 21 0.14 42 0.51
LS220_M135135_M0_L50 185 1.35 1.35 >\ 32.0 8.59 0.20 0.419 0.26 24 0.16 44 0.66
LS220_M135135_M0_LTE 185 1.35 1.35 >\ 21.1 14.24 0.11 0.000 0.20 17 0.13 30 0.24
LS220_M140120_LK 185 1.4 1.2 >\ 23.5 22.82 0.19 0.000 0.09 11 0.15 20 0.47
LS220_M140120_M0 185 1.4 1.2 >\ 24.8 23.38 0.24 0.001 0.18 14 0.15 29 0.63
LS220_M140120_M0_L05 185 1.4 1.2 >\ 28.2 20.33 0.28 0.017 0.18 14 0.15 29 0.76
LS220_M140120_M0_L25 185 1.4 1.2 >\ 28.5 12.07 0.41 0.625 0.16 13 0.20 29 1.97
LS220_M140120_M0_L50 185 1.4 1.2 >\ 31.2 13.97 0.70 8.289 0.18 13 0.22 28 3.97
LS220_M140140_LK 185 1.4 1.4 09.9 4.58 0.14 0.087 0.14 16 0.17 29 0.48
LS220_M140140_LK_HR 123 1.4 1.4 09.4 2.68 0.07 0.168 0.15 17 0.17 33 0.26
LS220_M140140_LK_LR 246 1.4 1.4 08.6 3.12 0.17 0.019 0.16 15 0.19 28 0.68
LS220_M144139_LK 185 1.44 1.39 07.2 3.91 0.19 0.014 0.14 15 0.16 30 0.54
LS220_M145145_LK 185 1.45 1.45 02.3 2.05 0.16 1.230 0.14 14 0.21 22 0.82
LS220_M150150_LK 185 1.5 1.5 00.9 0.16 0.03 0.001 0.08 11 0.19 13 0.13
LS220_M160160_LK 185 1.6 1.6 00.6 0.07 0.03 0.000 0.07 8 0.21 8 0.14
LS220_M171171_LK 185 1.71 1.71 00.5 0.06 0.03 0.000 0.08 9 0.22 8 0.14
SFHo_M1365125_LK 185 1.365 1.25 >\ 26.4 8.81 0.15 1.888 0.14 14 0.23 24 0.92
SFHo_M135135_LK 185 1.35 1.35 12.0 6.23 0.35 1.924 0.17 14 0.24 28 2.24
SFHo_M135135_LK_HR 123 1.35 1.35 06.9 1.78 0.23 0.982 0.17 14 0.21 29 1.23
SFHo_M135135_LK_LR 246 1.35 1.35 03.4 2.79 0.36 1.877 0.18 14 0.24 27 2.35
SFHo_M135135_M0 185 1.35 1.35 07.7 1.23 0.42 3.255 0.22 17 0.22 33 2.40
SFHo_M140120_LK 185 1.4 1.2 >\ 24.3 11.73 0.12 1.302 0.14 14 0.20 27 0.59
SFHo_M140120_M0 185 1.4 1.2 >\ 32.7 15.65 0.30 2.368 0.22 17 0.16 34 1.05
SFHo_M140140_LK 185 1.4 1.4 01.1 0.01 0.04 3.853 0.19 37 0.35 24 0.60
SFHo_M144139_LK 185 1.44 1.39 00.9 0.09 0.04 3.056 0.18 16 0.33 20 0.53
SFHo_M146146_LK 185 1.46 1.46 00.7 0.02 0.00 0.000 - - - - -

Table 2 reports ejecta masses, average electron fraction \langle Y_{e}\rangle, entropy \langle s\rangle, asymptotic velocity v_{\rm ej}, and kinetic energy T_{\rm ej} of the ejecta for all models. We also report the rms opening angle of the outflows about the equatorial plane \theta_{\rm rms}. Note that the definition of \theta_{\rm rms} implies that at least 75\% of the ejecta is confined within 2\theta_{\rm rms} of the orbital plane.

The dynamical ejecta masses are not fully converged in our simulations. From the comparison of results obtained at different resolutions we estimate the relative errors to be of the order of {\sim}50\%. In the following, we will assume the uncertainty in the ejecta masses to be

\Delta M_{\rm ej}=0.5\,M_{\rm ej}+(5\times 10^{-5})\,M_{\odot}\,. (17)

On the other hand, other ejecta properties, such as the average asymptotic velocities, the entropy, and the composition, appear to be converged (see Appendix A for a more detailed discussion).

Table 3Dynamical ejecta masses reported in the literature for the (1.35+1.35)\,M_{\odot} binary simulated with either the DD2, or the SFHo EOS. For each reported value we indicate the reference and whether the simulation included heating and compositional changes due to neutrino emission (\nu cool) and/or absorption (\nu abs). There are differences of factors of a few among published results.
EOS \nu cool. \nu abs. M_{\rm ejecta} Ref.
DD2 \checkmark - 0.3 This work
DD2 \checkmark \checkmark 1.4 This work
DD2 - - 3.1 Bauswein et al. (2013)
DD2 \checkmark - 0.6 Bovard et al. (2017)
DD2 \checkmark - 0.4 Lehner et al. (2016)
DD2 \checkmark - 0.9 Sekiguchi et al. (2015)
DD2 \checkmark \checkmark 2.1 Sekiguchi et al. (2015)
SFHo \checkmark - 3.5 This work
SFHo \checkmark \checkmark 4.2 This work
SFHo - - 4.8 Bauswein et al. (2013)
SFHo \checkmark - 3.5 Bovard et al. (2017)
SFHo \checkmark - 3.4 Lehner et al. (2016)
SFHo \checkmark - 10.0 Sekiguchi et al. (2015)
SFHo \checkmark \checkmark 11.0 Sekiguchi et al. (2015)

Despite the large inferred numerical uncertainty, we find overall good qualitative agreement between the ejecta masses estimated from our simulations and others presented in the literature. However, quantitative differences are present. In particular, we compare dynamical ejecta masses estimated for the (1.35+1.35)\,M_{\odot} binaries with the DD2 and the SFHo EOSs, which have been considerer in several previous works and by us. The results are reported in Table 3. Sekiguchi et al. (2015) and Sekiguchi et al. (2016) considered these binaries with and without the inclusion of neutrino absorption. In the former case they reported ejecta masses about {\sim}3 times larger than what we find. The disagreement is somewhat less severe for the simulations with neutrino absorption included. Lehner et al. (2016) and Bovard et al. (2017) performed simulations only with neutrino cooling and report ejecta masses in good agreement with those of our DD2_M135135_LK and SFHo_M135135_LK binaries. Note, however, that the latter study also employed the WhiskyTHC code, albeit with a different grid setup and lower resolution, so it does not represent a completely independent confirmation of our results. Finally, Bauswein et al. (2013) performed simulations using a smoothed particle hydrodynamics (SPH) code with an approximate treatment of general-relativistic (GR) and neglected the effect of neutrinos. They report ejecta mass of 0.48\times 10^{-2}\,M_{\odot} for the SFHo binary. This is in good agreement with our value of 0.35\times 10^{-2}\,M_{\odot}, especially when taking into account that neglecting neutrino cooling results in an overestimate of the ejecta mass (Radice et al., 2016b). On the other hand, the ejecta mass they report for the DD2 EOS is a factor {\sim}10 larger than what we infer from our simulations.

Figure 5.— Dynamical ejecta masses (left panel) and velocities (right panel) versus their fit according to Eqs. (18) and (22). The data points show the results from our fiducial subset of simulations, which have also been used to re-calibrate the fitting coefficients. For the velocity fit, we only include/plot models with total ejecta mass larger than 5\times 10^{-5} M_{\odot}. There is a correlation between M_{\rm ej} and v_{\rm ej} as extracted from our simulations and their predictors M_{\rm ej;fit} and v_{\rm ej;fit}.

Dietrich & Ujevic (2017) derived empirical formulas predicting the ejecta mass and velocity from binary NS mergers extending previous work on the ejecta from BHNS mergers (Foucart, 2012; Kawaguchi et al., 2016)111See Foucart et al. (2018b) for updated fits to BHNS merger data.. They calibrated their fitting formulas using data from Hotokezaka et al. (2013b), Bauswein et al. (2013), Dietrich et al. (2015), Sekiguchi et al. (2016), Lehner et al. (2016), and Dietrich et al. (2017b). Note that most of these data completely neglected neutrino effects (both emission and absorption). Moreover, as we discuss above, some of these works predict rather different ejecta masses for the same binaries. The resulting fits have been used in Abbott et al. (2017b) to estimate the contribution of the dynamical ejecta to the kilonova that followed GW170817 and, more recently, by Coughlin et al. (2018) to constrain the tidal deformability of the binary progenitor to GW170817.

We recalibrate the model of Dietrich & Ujevic (2017) using our data. We only consider the simulations performed without neutrino heating and at our reference resolution h=185\,{\rm m}, so as to have a homogeneous dataset. We refer to this subset of the runs as being our fiducial subset of simulations. The omission of neutrino re-absorption could result in a systematic underestimate of the ejecta mass by up to factors of a few (see Table 2). However, we expect the qualitative trends to be robust.

Following Dietrich & Ujevic (2017) we fit the ejecta mass as

\begin{split}&\displaystyle\frac{M_{\rm ej;fit}}{10^{-3}M_{\odot}}=\Bigg{[}% \alpha\left(\frac{M_{b}}{M_{a}}\right)^{1/3}\left(\frac{1-2C_{a}}{C_{a}}\right% )+\\ &\displaystyle\qquad\beta\left(\frac{M_{b}}{M_{a}}\right)^{n}+\gamma\left(1-% \frac{M_{a}}{M_{a}^{\ast}}\right)\Bigg{]}M_{a}^{\ast}+(a\leftrightarrow b)+% \delta\,,\end{split} (18)

where M_{a}, M_{b} and M_{a}^{\ast}, M_{b}^{\ast} are the NS gravitational and baryonic masses, respectively, and

C_{a,b}=\frac{GM_{a,b}}{c^{2}R_{a,b}} (19)

are the stars’ compactnesses. We use Eq. (17) to estimate the numerical uncertainty in the ejecta mass and we use a standard least square fit to determine the coefficients \alpha, \beta, \gamma, \delta, and n in Eq. (18). We find

\displaystyle\alpha=-0.657\,, \displaystyle\beta=4.254\,, \displaystyle\gamma=-32.61\,, (20)
\displaystyle\delta=5.205\,, \displaystyle n=-0.773\,. (21)

We find differences of {\sim}50{-}100\% compared to the values reported by Dietrich & Ujevic (2017). In particular, our fit systematically predicts lower dynamical ejecta masses compared to Dietrich & Ujevic (2017).

The right panel of Fig. 5 shows the ejecta masses from our fiducial subset of simulations plotted against the predicted ejecta masses from Eq. (18). Overall, we find that M_{\rm ej} correlates with M_{\rm ej;fit}, suggesting that the model of Dietrich & Ujevic (2017) captures some of the physics behind the mass ejection. However, there are differences between the values measured in simulations and those predicted by the fit as large as {\sim}200\%. This is a factor of a few in excess of the finite-resolution uncertainties we estimate for our simulations. We note that Dietrich & Ujevic (2017) also reported similarly large fit residuals. More worrisome is the fact that the model seems to be systematically missing trends in the data. For example it over predicts the ejecta for most of the massive DD2 binaries. This suggests that Eq. (18) does not capture all of the physical effects relevant for the mass ejection.

We fit the mass-weighted average asymptotic velocity v_{\rm ej} of the ejecta using an expression similar to the one used by Dietrich & Ujevic (2017):

v_{\rm ej;fit}=\left[\alpha\left(\frac{M_{a}}{M_{b}}\right)(1+\gamma C_{a})% \right]+(a\leftrightarrow b)+\beta\,. (22)

On the basis of the comparison between results obtained at different resolution, we use a constant value \Delta v_{\rm ej}=0.02\,c as fiducial uncertainty for the ejecta velocity and perform a least square fit of our fiducial subset of simulations using Eq. (22). When doing so, we exclude simulations with ejecta masses smaller than 10^{-5}\,M_{\odot}, since the ejecta properties cannot be reliably measured in those cases. Note that changing the value of \Delta v_{\rm ej} has no effect on the results of the fitting procedure. Using a standard least square method we find

\displaystyle\alpha=-0.287\,, \displaystyle\beta=0.494\,, \displaystyle\gamma=-3.000\,. (23)

These coefficients are in good agreement with those reported by Dietrich & Ujevic (2017) for the ejecta velocity in the orbital plane. This suggests that the ejecta velocity predicted by the simulations are robust. We also find good qualitative agreement between the ejecta velocities from simulations and those predicted by the model, as shown in the right panel of Fig. 5. However, there is a significant spread with residuals as large as 0.1\,c.

Figure 6.— Dynamical ejecta masses M_{\rm ej} and asymptotic velocities v_{\rm ej} as a function of the tidal parameter \tilde{\Lambda}. The data points show the results from our fiducial subset of simulations. No correlation is found between the ejecta mass and \tilde{\Lambda}. There is, however, a weak tentative inverse correlation between the asymptotic velocity of the ejecta and \tilde{\Lambda}.

Bauswein et al. (2013) compared dynamical ejecta mass obtained for the (1.35+1.35)\,M_{\odot} binary using different EOSs and found that it was inversely proportional to the NS radii suggesting that large ejecta masses could be produced in the case of compact NSs. This has been used in Nicholl et al. (2017) to suggest that the large inferred ejecta mass for GW170817 implies a small NS radius R_{1.35}\lesssim 12\,{\rm km}. The strongest support for the correlation found by Bauswein et al. (2013) comes from binaries that were simulated with an approximate treatment of thermal effects. Our results (Table 2) also indicate that softer EOSs might produce more ejecta. However, we do not find evidences for a clear correlation that would support the conclusions of Nicholl et al. (2017).

Figure 6 shows ejecta masses and velocities as a function of the tidal deformability of the binary \tilde{\Lambda} (e.g., Flanagan & Hinderer, 2008; Favata, 2014). By comparing results from simulations with different EOSs we confirm that there is some dependency of the ejecta mass on the stiffness of the EOS, with softer EOSs such as SFHo and LS220 producing more ejecta than BHB\Lambda\phi and DD2. We also find indications that prompt BH formation results in smaller than average dynamical ejecta masses. However, we are not able to identify a relation between the ejecta masses and properties of the EOS such as the radius of a reference NS or the tidal deformability. This is not too surprising given that most of the dynamical ejecta in our simulations is the result of the centrifugal bounce of the merger remnant. Its characteristics are likely to depend on thermal effects and details of the EOS at high densities that are not captured by \tilde{\Lambda}. Our data exclude the possibility of constraining \tilde{\Lambda} from the measurement of the properties of the dynamical ejecta alone. On the other hand, it might be possible to use kilonova observations to constrain other properties of the binary, for instance using an improved version of Eq. (18), and thus indirectly improve the bounds on \tilde{\Lambda} by restricting the priors used in the GW data analysis. However, these applications would necessarily require more precise theoretical predictions than those currently available.

Figure 7.— Mass of shock heated ejecta with s>10\,k_{\rm B} plotted against the mass of the cold ejecta. The center-left part of the figure is tentatively populated by stiff EOS. Tidally dominated ejecta dominates over shock-heated ejecta only for binaries with prompt BH formation (lower-middle part of the figure).

Other properties of the outflow, such as proton fraction and entropy, also show some dependency on the EOS. The reason is that different EOSs result in different relative amount of shocked and tidal ejecta. This is quantified in Fig. 7, where we show the mass of the tidal and shocked ejecta for our fiducial subset of simulations, i.e., those with only neutrino cooling and h=185\,{\rm m}. For this analysis we tentatively identify as shocked ejecta all of the unbound material crossing a coordinate sphere with radius \simeq 443\,{\rm km} and having specific entropy per baryon larger than 10\,k_{\rm B}. Material ejected with smaller entropies is assumed to be originating from tidal interactions. We find that stiff EOSs, such as BHB\Lambda\phi and DD2, typically have smaller M_{\rm tidal}^{\rm ej} than softer EOSs, such as LS220 and SFHo. Softer EOSs also eject more mass overall. The shocked component of the dynamical ejecta dominates the overall mass ejection for most of the binaries we have considered. The only exceptions are cases where BH formation occurs shortly after merger (\lesssim 1\,{\rm ms}) so that there is no centrifugal bounce of the remnant.

The balance between shocked and tidal ejecta is also dependent on the binary mass ratio (see Sec. 4.1). On the one hand, binaries with larger mass asymmetry produce more massive tidal outflow streams. On the other hand, asymmetric binaries result in the partial disruption of the secondary star during merger, less violent mergers, and smaller amount of shocked ejecta (Table 2; Hotokezaka et al. 2013b; Bauswein et al. 2013; Lehner et al. 2016; Sekiguchi et al. 2016; Dietrich et al. 2017b; Bovard et al. 2017). We conclude that the balance between tidal ejecta and shocked ejecta is set by a complex interplay between microphysical effects and binary properties.

Figure 8.— Average electron fraction \langle Y_{e}\rangle vs. rms angular spread \theta_{\rm ej} of the ejecta. The data points show the results from our fiducial subset of simulations. We only include models with total ejecta mass larger than 5\times 10^{-5} M_{\odot}. The shock-heated component of the ejecta is absent in the cases with prompt-BH formation.

Tidally-driven outflows have low average Y_{e} and are preferentially concentrated close to the orbital plane. Conversely, the shocked ejecta are spread over a larger angle of the orbital plane and have larger \langle Y_{e}\rangle, so a correlation between \theta_{\rm ej} and \langle Y_{e}\rangle is expected. This is shown in Fig. 8. For clarity, the figure only shows our fiducial subset of simulations, for which numerical and microphysical parameters have been set in a consistent way. The binaries in the lower left corner of the figure are LS220 binaries resulting in prompt BH formation. These are the same binaries appearing in the lower part of Fig. 7. The outflows from these binaries are dominated by the tidal ejection of material. The two outliers with \theta_{\rm ej}=20{-}25 degrees are the SFHo_M140140_LK and SFHo_M144139_LK, which also undergo prompt collapse, but have outflows mostly driven by shocks. With the possible exception of the binaries resulting in prompt BH formation, all others show a clear correlation between \langle Y_{e}\rangle and \theta_{\rm ej}. This suggests that a constrain on the opening angle of the dynamical ejecta, perhaps obtained by combining the observation of multiple systems with different orientations, could constrain the strength of the bounce of the massive NS after merger and the composition of the dynamical ejecta.

Overall, we find good qualitative agreement between our results for the dynamical ejecta and those reported in previous studies that adopted a more idealized treatment for the EOS of NSs and/or approximate GR (Hotokezaka et al., 2013a; Bauswein et al., 2013; Dietrich et al., 2017b, a). At the same time, there are substantial quantitative differences in the dynamical ejecta mass and properties from those studies, as well as with other works that considered a limited number of binary configurations, but included full-GR and neutrino effects (Sekiguchi et al., 2015, 2016). The total ejecta mass is also not fully converged in our simulations (Table 2). While there appear to be robust trends and correlations between ejecta properties, binary parameters, and the EOS of NSs, more comprehensive and better resolved studies would be needed to create reliable quantitative models of the dynamical ejecta.

3.1.2 Fast Moving Ejecta

Figure 9.— Cumulative distribution of ejecta velocities for the (1.35+1.35)\,M_{\odot} binary with four EOSs and at the reference resolution. Neutrino re-absorption has not been included in these simulations.

Metzger et al. (2015) re-analyzed data from Bauswein et al. (2013) and identified 106 SPH particles, corresponding to {\sim}10^{-4}\,M_{\odot} of material, that were dynamically ejected with velocities in excess of 0.6\,c. If indeed present, these fast moving ejecta would expand sufficiently rapidly to still contain a significant fraction of free neutrons at freeze out. The decay of the neutrons in the outermost part of the ejecta would then produce a bright UV/optical counterpart to the merger on a timescale of several minutes to an hour (Kulkarni, 2005; Metzger et al., 2015). On the other hand, because of the small number of SPH particles, it cannot be excluded that this fast component of the ejecta is due to numerical noise, as also recognized by Metzger et al. (2015).

More recently, a fast moving component of the dynamical ejecta was proposed as a possible origin for the synchrotron radiation detected from GW170817 in the first {\sim}100 days (Mooley et al., 2017; Hotokezaka et al., 2018b). This interpretation is currently disfavored on the light of more recent observations showing an abrupt decline in the source luminosity in all bands (Alexander et al., 2018) and VLBI observations showing apparent superluminal motion of the radio source indicative of collimation of the outflow (Mooley et al., 2018; Ghirlanda et al., 2018). Nevertheless, such fast moving component of the outflows was identified in Hotokezaka et al. (2018b) using data from the simulations of Kiuchi et al. (2017). The latter simulations, however, employed piecewise polytropic fits to the cold NS EOS augmented with an ideal-gas component to describe thermal effects. These approximations affect the thermodynamical properties of the shocks responsible for the ejection of this material (Bauswein et al., 2013), so they need to be independently verified.

Our previous simulations (Radice et al., 2016b) did not show evidences for the presence of a fast moving component of the ejecta. However, in Radice et al. (2016b) we only considered one equal mass configuration with M_{a}=M_{b}=1.39\ M_{\odot} simulated using the LS220 EOS. While, Metzger et al. (2015) considered a large sample of EOSs and binary masses. In our new simulations we find evidence for a fast moving component of the outflow with asymptotic velocities in excess of 0.6\,c (Table 2). The amount of the fast moving ejecta strongly depends on the EOS and other binary parameters. For instance, for total binary masses up to 2.8\,M_{\odot}, with the exceptions of the simulations performed with viscosity, discussed in a companion paper (Radice et al., 2018b), the LS220 EOS does not seem to predict appreciable amount of fast outflows, in agreement with Radice et al. (2016b). On the other hand, binaries simulated with the BHB\Lambda\phi, DD2, or SFHo EOSs, as well as higher-mass LS220 binaries, typically eject {\sim}10^{-6}{-}10^{-5}\,M_{\odot} of fast-moving material. This is still one or two orders of magnitude less than reported in Metzger et al. (2015), but suggest that the ejection of at least a small amount of fast material does indeed take place during NS mergers. That said, given the small overall mass involved, we cannot completely exclude its origin as being numerical. Indeed, we find large variations in the amount of the fast moving ejecta with resolution and much better resolved simulations would be needed to fully quantify the mass of these fast outflows. On the other hand, the smooth distribution of the ejecta in velocity space, shown in Fig. 9, seem to suggests that the properties of the outflows are reasonably well captured even down to these masses.

We discuss the possible observational consequences of the fast moving component of the ejecta in Sec. 5. Here, we focus on their origin. Metzger et al. (2015) suggested that shocks generated at the time of the collision between the NSs could be accelerated by the steep density gradient at the surface of the remnant and drive outflows with trans-relativistic velocities. This is a scenario first considered by Kyutoku et al. (2014) who studied its possible high-energy signatures. The same scenario has also been recently revisited in detail by Ishii et al. (2018). They used high-resolution 1D Lagrangian simulations to resolve the shock acceleration and predict the amount of ejected material. They found that this mechanism can indeed produce free neutrons. However, they found the amount of free neutrons that can be produced in this way to be {\sim}10^{-7}{-}10^{-6}\,M_{\odot}. This might be insufficient to explain the results of Metzger et al. (2015). On the other hand, the values found by Ishii et al. (2018) are not inconsistent with our data.

[1em] SFHo_M135135_M0ppppppppppp

Figure 10.— Angular distribution of the flux of ejecta with asymptotic velocities larger than 0.6\,c for the SFHo_M135135_M0 binary. The data is extracted on a coordinate sphere with radius 443\ \mathrm{km}. For clarity, the time is retarded in the same way as in Fig. 3. We find two distinct episodes resulting in the production of fast ejecta: one associated with the merger and one associated with the first bounce of the remnant.

Comparable mass binaries simulated with the SFHo EOS, the softest in our set, result in the most violent mergers. As the NSs collide material is squeezed out from the collisional interface at large velocities, as in the simulation discussed by Metzger et al. (2015). This is evident from Fig. 3, where we show the outflow rate of the fast-moving component of the ejecta with asymptotic velocities v_{\rm ej}>0.6\,c. Because of its peculiar origin, this outflow component is preferentially channeled to high latitudes and in the directions not obstructed by the NSs. This is shown, in the case of the SFHo_M135135_M0 binary, in the upper panel of Fig. 10. However, this first clump of material amounts only to about a quarter of all of the fast ejecta. Most of the fast moving ejecta, instead, appears to originate when the shock launched from the first bounce of the remnant breaks out of the forming ejecta cloud, as indicated in Fig. 3. This second component of the fast ejecta is also highly anisotropic, but preferentially concentrated close to the orbital plane, as shown in the lower panel of Fig. 10 for the SFHo_M135135_M0 binary. This is possibly because of the oblate shape of the merger debris cloud that favors the acceleration of the material close to the orbital plane.

Figure 11.— Angular distribution and asymptotic velocity of the ejecta for the (1.35+1.35)\,M_{\odot} and (1.4+1.2)\,M_{\odot} binaries with the SFHo EOS simulated without neutrino re-absorption. \theta=0 corresponds to the polar axis, while \theta=90^{\circ} is the orbital plane. The data is extracted on a coordinate sphere with radius 443\ \mathrm{km}. The fast-moving tail of the ejecta is confined to a region close to the orbital plane for the unequal mass SFHo_M140120_LK binary (right panel), but it is somewhat more isotropically distributed for the equal mass SFHo_M135135_LK binary (left panel).

Smaller mass ratio binaries, or binaries simulated with other EOSs, do not however show evidences for an early fast-ejecta component from the collisional interface between the NSs. For these binaries the fast-moving component of the ejecta is entirely constituted of she material accelerated after the bounce of the merger remnant. This material is preferentially confined to the region close to the orbital plane, as shown in Fig. 11. We stress that this angular distribution is inconsistent with the ejecta originating at the collisional interface between the NSs. The reason for the different behavior of more unequal mass binaries and/or binaries with stiffer EOSs is that the former result in less violent collisions than those predicted for comparable mass NSs with the SFHo EOS, either because of the larger NS radii, or because of the partial disruption of the secondary star shortly before merger. We note that the inclusion of neutrino heating only results in a modest quantitative correction to the distribution of the fast-moving component of the ejecta, but the overall qualitative picture is unchanged.


[1em] BHBlp_M135135_LKpppppppppppp

[1em] SFHo_M135135_LKppppppppppp

Figure 12.— Angular distribution of the ejecta with asymptotic velocities larger than 0.6\,c for the (1.35+1.35)\ M_{\odot} binaries simulated without neutrino re-absorption, and at the fiducial resolution. The data is extracted on a coordinate sphere with radius 443\ \mathrm{km}. The LS220_M135135_LK binary is not shown, since it does not produce an appreciable amount of fast moving ejecta. The fast-moving component of the outflow is anisotropic and covers a relatively small portion of the sky around the binary, especially in the case of the BHB\Lambda\phi and DD2 binaries.

We find the distribution of the fast-component of the ejecta to be not only anisotropic in latitude, but also in the azimuthal direction, as shown in Fig. 12. The fast-moving material forms narrow streams that cover only a small fraction of area of a sphere centered at the location of the binary merger. In contrast, the overall ejecta distribution is almost uniformly spread in the azimuthal direction (e.g., Bovard et al., 2017). We remark that non-spinning equal mass binaries have a discrete rotational symmetry of 180^{\circ} around the orbital angular momentum axis. However, this symmetry is typically broken at the time of the merger, when turbulence operating in the shear layer between the two NSs can exponentially amplify small initial asymmetries, such as those due to floating point truncation in the simulations (Paschalidis et al., 2015; Radice et al., 2016a). This is reflected in the asymmetry of the fast-moving tail of the dynamical ejecta and is another indication of the fact that the bulk of the fast-moving ejecta is launched with some delay from the merger. Note that the early fast ejecta for the equal mass binaries with the SFHo EOS, which instead originates at the time of the merger, is symmetric (see Fig. 10).

Figure 13.— Fast moving component of the ejecta plotted against the binary tidal deformability \tilde{\Lambda}. The data points show the results from our fiducial subset of simulations. Note that several binaries produce no appreaciable amount of fast moving ejecta and are not seen in the figure. There is only a tentative correlation between the mass of the fast moving ejecta and the tidal deformability of the binary.

We report the total mass of the ejecta with v_{\rm ej}>0.6\,c in Fig. 13. The error bars are estimated following Eq. (17). The motivation for this choice is that, while we find relatively large numerical uncertainties in the total ejecta mass, the relative distribution of the ejecta as a function of velocity appears to be robust (see Sec. A). Consequently, the numerical uncertainty in the total ejecta mass should be well correlated with that of the fast component. Each simulation is labelled by the tidal deformability of the binary \tilde{\Lambda}. The tidal parameter \tilde{\Lambda} is related to the compactness of the binary and should correlate with the strength with which different NSs binaries collide at the time of merger. We might expect that \tilde{\Lambda} should also correlate with the amount of fast ejecta produced during the collision between the two stars. Indeed, we find that there is a weak correlation between the two in Fig. 13. However, there are systematic effects that are not captured by \tilde{\Lambda}. For instance, the BHB\Lambda\phi and the DD2 binaries have very close compactnesses, since the two EOSs are identical for NSs with gravitational mass M\lesssim 1.5\,M_{\odot}. However, the BHB\Lambda\phi binaries typically produce a significantly larger amount of fast moving ejecta compared to the DD2 binaries. The reason is likely the more violent bounce of the merger remnant due to the formation of \Lambda hyperons in the aftermath of the BHB\Lambda\phi binary mergers (Radice et al., 2017). Perhaps more importantly, there are several binaries, spanning a wide range of \tilde{\Lambda}, that do not produce any appreciable amount of fast moving ejecta.

3.2. Secular Ejecta

Part of the tidal tails from the NSs remains bound to form a rotationally supported disk around a central remnant more precisely defined below. In the cases in which the former survives for more than about one millisecond, we observe the formation of hot {\sim}10{-}20\,{\rm MeV} streams of material expelled from what is originally the interface region between the NSs. This material assembles into an excretion disk. Consequently, there is a correlation between the life time of the remnant and the remnant disk mass (Radice et al., 2018c), see also Table 2.

Figure 14.— Electron fraction and density in the xz-plane for the SFHo_M135135_M0 simulation. The thin-black contours enclose regions with density larger than 10^{6}, 10^{7}, 10^{8}, 10^{9}, 10^{10}, and 10^{11} {\rm g}\ {\rm cm}^{-3}. The thick magenta light encloses the region with density larger than 10^{13} {\rm g}\ {\rm cm}^{-3}. This figure should be compared with the last panel of Fig. 2 where a volume rendering of the same data is shown.

The disks are geometrically thick and moderately neutron rich (see Fig. 14; Siegel & Metzger 2018). We observe the propagation of m=2 spiral density waves originated as the streams from the massive NS remnant impact the disk. After the first 10{-}20 ms from the merger, these streams subside, and m=1 spiral density waves, induced by the one-armed spiral instability of the remnant NS (Paschalidis et al., 2015; East et al., 2016; Radice et al., 2016a), become dominant.

If the merger does not result in the formation of a BH within few dynamical timescales, the remnant is composed of a relatively slowly rotating inner core surrounded by a rotationally supported envelope (Shibata et al., 2005; Kastaun et al., 2016; Hanauske et al., 2017; Ciolfi et al., 2017). In particular, regions with rest-mass densities below \simeq 10^{13}\,{\rm g}\,{\rm cm}^{-3} are mostly rotationally supported (Hanauske et al., 2017). Note that the rotational structure of the remnant might be strongly affected by the effective shear viscosity arising due to the turbulent fluid motion in the remnant (Radice, 2017; Shibata et al., 2017b; Kiuchi et al., 2018; Fujibayashi et al., 2018; Radice et al., 2018a). In our analysis we define as the central part of the remnant the region with \rho_{0}\geq 10^{13}\,{\rm g}\ {\rm cm}^{-3}, and we estimate the remnant disk mass M_{\rm disk} from the integral of the rest-mass density of the region with \rho_{0}<10^{13}\ {\rm g}\ {\rm cm}^{-3}. We remark that the same criterion has also been adopted by Shibata et al. (2017a). Furthermore, as shown in Fig. 14, the density threshold 10^{13}\,{\rm g}\,{\rm cm}^{-3} does indeed approximately corresponds to the boundary of the centrally condensed remnant. Our results are given in Table 2. Unfortunately, the 3D output necessary to estimate the disk mass in post processing has been accidentally deleted for six of our simulations. For those cases the disk masses are not given.

Figure 14 shows the electron fraction and the density contours shortly after merger, during the formation of the accretion disk. At this time the disk has not yet reached its maximum extent and is expanding behind the cloud of the dynamical ejecta. The former can be recognized for their higher electron fraction and are located at radii \gtrsim 75\,{\rm km}. The neutron rich outflow visible at radii \gtrsim 100 km is part of the tidal tail, while the higher Y_{e} material between the forming accretion disk and the tidal tail is part of the outflow generated during the first bounce of the merger remnant.

Figure 15.— Remnant disk masses as a function of the tidal parameter \tilde{\Lambda}. The data points show the results from our fiducial subset of simulations. Disk formation is suppressed in the case of prompt BH formation, corresponding to small \tilde{\Lambda}’s. The final disk masses saturate for large values of \tilde{\Lambda}.

We find that the remnant disk masses tightly correlated with the tidal deformability of the binary \tilde{\Lambda}. This is shown in Fig. 15 where we plot the disk masses for our fiducial models. The error bars are estimated from the comparison of simulations performed at different resolutions (see Table 2). In particular, we estimate the uncertainty on the disk mass due to numerical errors to be {\sim}30\%. To be conservative, we estimate the uncertainties on the disk masses as

\Delta M_{\rm disk}=0.5\,M_{\rm disk}+(5\times 10^{-4})\,M_{\odot} (24)

We remark that the error bars only account for finite-resolution uncertainties and we cannot exclude that missing physics, or more extreme binaries with larger mass asymmetries and/or extreme spins, could deviate from the trend shown in Fig. 15. Moreover, because a significant fraction ({\sim}30{-}90\%) of the disk is accreted promptly after BH formation, the transition between low and high mass disks visible in Fig. 15 would likely become sharper if we could evolve all of the hypermassive remnants to collapse.

Notwithstanding these caveats, we find that our data is reasonably well fitted with the simple expression

\frac{M_{\rm disk}}{M_{\odot}}=\max\left\{10^{-3},\alpha+\beta\tanh\left(\frac% {\tilde{\Lambda}-\gamma}{\delta}\right)\right\} (25)

with fitting coefficients \alpha=0.084, \beta=0.127, \gamma=567.1, and \delta=405.14. In the case of GW170817, the amount of ejecta estimated from the modeling of the kilonova signal {\sim}0.05\ M_{\odot} is about an order of magnitude too large to be explained by the dynamical ejecta. This suggests that most of the material powering the kilonova should have been produced during the viscous evolution of the accretion disk.

This finding is in agreement with the results from long-term GRMHD simulations of postmerger disks that show that up to {\sim}40\% of the accretion disk can be ejected over the viscous timescale (Siegel & Metzger, 2017; Fernández et al., 2018). According to this scenario, M_{\rm disk} should have been larger than at least {\sim}0.1\ M_{\odot}. This, in turn, implies that \tilde{\Lambda} for GW170817 should have been larger than about 400. Another possibility, that, however, we cannot presently test with our data, is that the progenitor binary to GW170817 had a large mass asymmetry. Under these conditions it has been suggested that the merger could still produce a massive disk even for compact configurations (Shibata et al., 2003; Shibata & Taniguchi, 2006; Rezzolla et al., 2010). On the other hand, large mass asymmetries are disfavored on the basis of mass measurements for galactic double pulsars (Ozel & Freire, 2016). Moreover, Shibata & Taniguchi (2006) found that the accretion disk mass only increases by about an order of magnitude for extremely asymmetric binaries.

Figure 16.— Dynamical ejecta M_{\rm ej;dyn} versus secular ejecta masses M_{\rm ej;sec}. With the exception of the prompt BH formation cases that are able to expel at least a few 10^{-4}\,M_{\odot} in dynamical ejecta, the secular ejecta dominate over the dynamical ejecta.

We speculate on the total amount of mass expelled in a BNS merger. We consider three different ejection channels. In addition to the dynamical ejecta, directly extracted from the simulations, we include the possible presence of neutrino- and magnetically-driven winds from the disk, which develop on a time scale of a few tens of ms. We parametrize their contribution to the ejecta as a fraction of the disk mass. The neutrino-driven wind is estimated as M_{\rm ej;wind}=\xi_{\rm wind}M_{\rm disk}, with \xi_{\rm wind}=0.03\pm 0.015. The uncertainty includes variations due to the possible presence of a longer-lived remnant (e.g., Martin et al., 2015; Just et al., 2015). The viscous ejecta is taken to be M_{\rm ej;vis}=\xi_{\rm vis}M_{\rm disk} with \xi_{\rm vis}=0.2\pm 0.1, a range including the results from most postmerger simulations (e.g., Metzger & Fernández, 2014; Just et al., 2015; Siegel & Metzger, 2018; Fujibayashi et al., 2018; Fernández et al., 2018).

In Fig. 16 we compare dynamical and secular ejecta masses, the latter estimated as M_{\rm ej;sec}=M_{\rm ej;wind}+M_{\rm ej;vis}. With the exception of the prompt-BH formation cases that produce at least 10^{-4}\,M_{\odot} of dynamical ejecta, the total ejecta is largely dominated by the disk ejecta. As a consequence of the tight correlation between the disk mass and \tilde{\Lambda}, we expect a correlation between the total ejecta and \tilde{\Lambda}. Indeed, our estimates for M_{\rm ej}+M_{\rm ej;wind}+M_{\rm ej;vis} are reasonably well fitted by the same simple formula used for M_{\rm disk}, Eq. (25), but assuming a floor value of 5\times 10^{-4}M_{\odot} and with fitting coefficients \alpha=0.0202, \beta=0.0341, \gamma=538.8, and \delta=439.4. In most of the cases, the relative error between the total ejecta mass and the fit is below 50%. When prompt BH formation occurs, we do not expect the total ejecta to be larger than \sim 10^{-3}M_{\odot}. On the other hand, for long-lived remnants the mass of the unbound material can span a broad range of masses: from a few times 10^{-3}M_{\odot} to 0.1M_{\odot}, increasing with \tilde{\Lambda}.

In comparison to previous studies (Oechslin et al. 2006; Hotokezaka et al. 2013b, as reported in Wu et al. 2016) we find that disk winds and viscous outflows should contribute a significantly larger fraction of the overall ejecta. The difference can be explained in part by the fact that previous simulations did not include the effects of neutrino cooling, or used approximate treatments for the gravity, and in part by the fact that we assume that a larger fraction of the accretion disk can be unbound secularly compared to Wu et al. (2016). The latter is motivated by recent GRMHD simulations that showed that up to {\sim}40\% of the disk can be unbound secularly (Siegel & Metzger, 2017; Fernández et al., 2018). Our study differs from some of the previous studies also in the numerical setup and analysis methodology.

4. Nucleosynthesis

The discovery of a kilonova counterpart to GW170817 (Abbott et al., 2017e; Arcavi et al., 2017; Chornock et al., 2017; Cowperthwaite et al., 2017; Coulter et al., 2017; Drout et al., 2017; Evans et al., 2017; Kasliwal et al., 2017; Nicholl et al., 2017; Smartt et al., 2017; Soares-Santos et al., 2017; Tanvir et al., 2017; Tanaka et al., 2017) provided compelling evidence that NS mergers are one of the main sources of r-process elements in the Universe (Kasen et al., 2017; Rosswog et al., 2018; Hotokezaka et al., 2018a). However, the question of whether NS mergers are the only source of r-process elements or whether is a contribution from other sources is still not completely settled. Part of the uncertainty is due to the lack of a full theoretical understanding of the nucleosynthetic yields from mergers. Here, we study in detail the dependency of the r-process nucleosynthesis on the properties of the binary, mostly focusing on the dynamical ejecta.

4.1. Dynamical Ejecta

For simplicity, we perform most of our nucleosynthesis calculations using the approach we developed in Radice et al. (2016b) and that we briefly recall. We extract electron fraction Y_{e,R}, specific entropy s_{R}, velocity v_{R}, and rest-mass density \rho_{0,R} of the unbound material crossing a coordinate sphere surface with radius R=300\,G/c^{2}\,M_{\odot}\simeq 443\,{\rm km}. A fluid element is considered to be unbound if its kinetic energy is sufficient to overcome the gravitational potential well, that is, u_{t}\leq-1, where we have assumed a nearly stationary metric. See Kastaun & Galeazzi (2015) and Bovard et al. (2017) for possible alternative criteria. For the nucleosynthesis calculations we further assume that the outflow is undergoing an homologous expansion

\rho_{0}(t)=\rho_{0,R}\left(\frac{v_{R}}{cR}t\right)^{-3}\,. (26)

This density history is matched with the expansion histories used in the parametrized r-process calculations of Lippuner & Roberts (2015):

\rho_{0}(t)=\rho_{0}\big{(}s,Y_{e},T=6\,{\rm GK}\big{)}\left(\frac{3\tau}{et}% \right)^{3}\,, (27)

where e\simeq 2.718 is Euler’s number. From the matching we extract the expansion timescale \tau. To compute the nucleosynthesis yields we bin the ejecta according to their density, entropy, and expansion time scale. The nucleosynthesis yields in each bin can be pre-computed, because, under the assumption of homologous expansion, the r-process outcome depends only on \rho_{0,R}, s_{R} Y_{e,R}, and \tau_{R}. Then we compute the full nucleosynthetic yields by summing up the contributions from all bins. The uncertainty due to this procedure and the comparison with nucleosynthetic calculations performed using Lagrangian tracer particles are discussed in Sec. 4.3.

Figure 17.— Electron fraction of the ejecta before the activation of the r-process (left panel) and final nucleosynthetic abundances (right panel). The histograms show the results from our fiducial subset of simulations. We only include models with total ejecta mass larger than 5\times 10^{-5} M_{\odot}. The green dots in the right panel show the solar abundances from Arlandini et al. (1999). All abundance curves are normalized by fixing the overall fraction of elements with 180\leq A\leq 200.

We discuss first the systematics of the nucleosynthetic yields for our fiducial subset of runs, for which neutrino re-absorption has been neglected. The effect of neutrino re-absorption is discussed below. Fig. 17 reports histograms of the electron fraction Y_{e} of the ejecta and the relative abundances of different isotopes synthesized by the r-process 32 years after the merger. We show the electron fraction since it is the most important variable determining the outcome of the r-process in the conditions relevant for NS mergers (e.g., Lippuner & Roberts, 2015; Radice et al., 2016b). In the figure we compare isotopic abundances of r-process elements estimated from our simulation with the solar abundances from Arlandini et al. (1999). Both abundances are normalized by fixing the overall fraction of elements with 180\leq A\leq 200.

In absence of neutrino heating, the distribution of the ejecta as a function of Y_{e} drops rapidly for Y_{e}>0.15{-}0.2. We find systematic changes in the composition of the outflows as the total mass of the binary increases. However, these variations are inconsistent between the different EOSs. The BHB\Lambda\phi and SFHo binaries produce higher Y_{e} ejecta as the binary mass is increased. Presumably because the larger compactness of the more massive binaries inhibits the production of tidally-driven ejecta. As a consequence, these binaries produce more of the light r-process elements, i.e., the isotopes with 90\lesssim A\lesssim 125, when approaching the prompt BH-formation threshold. Conversely, the relative amount of shocked ejecta from the LS220 binaries progressively decreases as the binary mass increases and the ejecta for the most massive binaries appears to be dominated by the tidal tail. Consequently, the relative abundance of light r-process elements decreases with the total binary mass for the LS220 binaries. Finally, the behavior of the DD2 binaries is non monotonic with mass: for chirp masses up to {\sim}1.25\,M_{\odot} the average electron fraction \langle Y_{e}\rangle increases with the mass. However, the highest mass binary DD2_M160160_LK, has the lowest \langle Y_{e}\rangle and the lowest relative abundance of light r-process elements.

Overall, we find that, irrespective of binary parameters and EOS, the dynamical ejecta robustly produce second- and third-peak elements, i.e., elements with 125\lesssim A\lesssim 145 and 185\lesssim A\lesssim 210, respectively. The reason for this is that the bulk of the outflows are always very neutron rich, with \langle Y_{e}\rangle<0.2. On the other hand, we observe variance in the production of the light r-process elements.

The normalization condition we impose effectively scales the second and third peaks of the solar abundances to the second and third peaks of the calculated abundances (see Figure 17). With this normalization, it is clear that our calculations underproduce the material in the mass {\sim}140 region beyond the second r-process peak. We note that although different models produce a wide range of Y_{e} distributions, this underproduction is robust. Therefore, it is more likely due to the choice of nuclear input data used in our network calculations. In particular, the choice of fission fragment distributions and neutron-induced fission rates can have a significant impact on abundances in this region (Eichler et al., 2015), and the symmetric fission fragment distributions used in our calculations likely underproduce material in this region.

Figure 18.— Dynamical ejecta sensitivity to the neutrino treatment. Left panel: electron fraction. Right panel: nucleosynthetic yields. All abundance curves are normalized by fixing the overall fraction of elements with 180\leq A\leq 200. Neutrino absorption can have a strong impact on the ejecta composition and on the nucleosynthetic yields around the light r-process elements.

Neutrino-matter interaction rates roughly scale with the square of the incoming neutrino energy. Consequently, the composition, nucleosynthetic yield, and EM opacity of the ejecta depend on the details of the neutrino radiation spectra (Foucart et al., 2016b). The determination of the latter is not possible using an energy integrated scheme as is the one adopted for this work. To quantify the relative uncertainty, we simulate the (1.35+1.35)\,M_{\odot} binary using the LS220 EOS and two different schemes for the calculation of the absorption opacities. In addition to the standard treatment in which we approximatively account for the incoming neutrino energy using Eq. (8), we also perform a simulation in which the absorption cross-section is calculated assuming thermal equilibration between matter and neutrinos. The latter simulation is labelled LS220_M135135_M0_LTE.

The results of this comparison are shown in Fig. 18. As anticipated, the inclusion of compositional changes and heating due to the absorption of neutrinos results in a significant shift in the ejecta Y_{e} distribution. The tidal tail is reprocessed to slightly higher values of Y_{e} and a new ejecta component, with Y_{e} up to 0.4 is found. The approximate inclusion of non-LTE effects results in a slight increase in the absorption rates. This is expected because the incoming neutrinos originally decoupled from the central regions of the remnant and the accretion disk at higher temperatures than those found in the ejecta. Overall, we find that the dynamical ejecta, especially in the regions close to the orbital plane, robustly produces heavy r-process elements with relative abundances close to solar. On the other hand, the neutrino re-absorption significantly affects the production of elements in the region close to the first peak of the r-process relative abundance pattern, i.e., around A=80. The treatment of the absorption opacities can result in changes in the relative abundance of elements in this region of up to factors of a few. This points to the need for more sophisticated simulations with spectral neutrino transport.

We do not consider the effect of neutrino oscillations in this work. However, we point out that recent studies have shown that NS merger remnants might host ideal conditions for resonant oscillations and fast flavor conversion (Zhu et al., 2016; Frensel et al., 2017; Deaton et al., 2018). This might have an impact on the nucleosynthesis of light r-process elements (Wu et al., 2017). The investigation of these effects in the context of dynamical simulations is urgent.

Figure 19.— Dynamical ejecta sensitivity to the binary mass ratio. Left panel: electron fraction. Right panel: nucleosynthetic yields. All abundance curves are normalized by fixing the overall fraction of elements with 180\leq A\leq 200. The nucleosynthetic yield from the dynamical ejecta is sensitive to the mass ratio in the region of the light r-process elements A\lesssim 120.

The nucleosynthetic yield of the dynamical ejecta is also sensitive to the mass ratio of the binary, especially in the region of light r-process elements. This is shown in Fig. 19 where we compare composition and nucleosynthesis yeids between the LS220_M135135_M0 and LS220_M140120_M0 binaries which are representative of the general trend with mass ratio. The relative abundances are normalized as in Fig. 17 over the mass range 180\leq A\leq 200. Because of the more abundant tidal tail the relative composition of the ejecta is shifted to lower Y_{e}. This affects the relative abundances of first-peak r-process elements which are reduced by factors of several. A similar trend is also present in the simulations that do not include neutrino re-absorption, see Table 2.

Figure 20.— Dynamical ejecta sensitivity to the EOS. Left panel: electron fraction. Right panel: nucleosynthetic yields. All abundance curves are normalized by fixing the overall fraction of elements with 180\leq A\leq 200. The total ejecta mass shows depends sensitively on the EOS, however electron fraction and nucleosynthetic yields appear to be insensitive to the EOS.

Even though there are differences in the quality and quantity of the dynamical ejecta, the r-process nucleosynthesis appears to be only weakly sensitive to the EOS. In Fig. 20 we show electron fraction distributions and final isotopic abundances for the dynamical ejecta from the (1.35+1.35)\,M_{\odot} binary simulated with different EOSs. Neutrino re-absorption has been included in these simulations using the M0 scheme. Irrespective of the EOS, the dynamical ejecta is neutron rich \langle Y_{e}\rangle\lesssim 0.25, but has a broad distribution in Y_{e}, with a significant amount of ejecta with Y_{e} as large as 0.4. Consequently, both heavy and light r-process elements are produced. For comparable mass binaries, as are those shown in Fig. 20, we find that the dynamical ejecta robustly synthesizes r-process elements with isotopic abundances close to solar. Conversely, the dynamical ejecta from binaries with large mass asymmetry underproduces light r-process elements. These conclusions are not altered with the inclusion of viscosity. Indeed, while viscosity impacts the overall ejecta mass (Table 2 and Radice et al. 2018b), it does not affect the electron fraction and the isotopic abundances of the r-process in the dynamical ejecta. Overall, our simulations suggest that the EOS of dense matter controls the amount of the dynamical ejecta, but that the relative isotopic abundances are most sensitive to weak reaction rates, mass ratio, and total binary mass.

4.2. Secular Ejecta

The nucleosynthetic yields of the secular ejecta have been the subject of several recent studies. In the case of neutrino-driven winds, neutrino irradiation of the expanding ejecta drives the electron fraction towards higher values than those of the original cold weak equilibrium configuration. The electron fraction of the material in the wind depends on the relative timescales for weak re-equilibration and expansion. If the expansion is not too rapid, then the material achieves electron fractions given by the weak equilibrium in optically thin conditions with neutrinos (Qian & Woosley, 1996):

\left(Y_{e}\right)_{\rm eq}=\left[1+\frac{L_{\bar{\nu}_{e}}\left(\epsilon_{% \bar{\nu}_{e}}-2\Delta+1.2\Delta^{2}/\epsilon_{\bar{\nu}_{e}}\right)}{L_{{\nu}% _{e}}\left(\epsilon_{{\nu}_{e}}+2\Delta+1.2\Delta^{2}/\epsilon_{{\nu}_{e}}% \right)}\right]^{-1}\,, (28)

where L_{\nu_{e}} and L_{\bar{\nu}_{e}} are the luminosities for electron neutrinos and antineutrinos, \epsilon_{\nu}\equiv\langle E^{2}_{\nu}\rangle/\langle E_{\nu}\rangle, E_{\nu} is the neutrino energy, \langle x\rangle the average of x over the neutrino distribution function, and \Delta is the mass difference between neutrons and protons in vacuum. During the early post merger phase, L_{\bar{\nu}_{e}}\gtrsim L_{\nu_{e}}, while electron antineutrinos are significantly hotter than neutrinos (\langle E_{\bar{\nu}_{e}}\rangle\approx 15~{}{\rm MeV}>\langle E_{\nu_{e}}% \rangle\approx 10~{}{\rm MeV}, e.g. Dessart et al. 2009; Perego et al. 2014; Foucart et al. 2016b). Thus, \left(Y_{e}\right)_{\rm eq}\lesssim 0.45 and r-process nucleosynthesis can occurs. However, due to the small neutron-to-seed ratio, only nuclei with A\lesssim 130 can be synthetized and this ejecta is expected to contribute to the first r-process peak. If the dynamical ejecta is dominated by the tidal component and is poor of first peak r-process elements, the neutrino-driven wind can complement the nucleosynthesis and lead to the production of all r-process elements (Martin et al., 2015).

If a massive disk forms after the merger, viscously-driven ejection, occurring over the disk lifetime, can become the dominant source of ejecta from a BNS merger (see Sec. 3.2). Viscous hydrodynamics simulations of the long term disk evolution, mainly performed in axisymmetry assuming a BH-torus system, showed that, if the disk becomes transparent to neutrinos, the rapid decrease in temperature makes neutrino cooling inefficient and the disk becomes convectively unstable (e.g. Fernández & Metzger, 2013; Just et al., 2015). The resulting large scale mixing, combined with the long time scale over which neutrino-matter interactions can occur, produces a rather uniform, broad distribution of Y_{e} in the ejecta. In particular, it was found that the resulting distribution has 0.1\lesssim Y_{e}\lesssim 0.45 and all r-process elements from the first to the third peak, as well as Uranium and Thorium, can be synthesized in proportions close to solar (e.g. Wu et al., 2016).

Recent 3D GRMHD simulations of a BH disk torus (Siegel & Metzger, 2017; Fernández et al., 2018) confirmed the presence of a self-regulating mechanism based on electron degeneracy in the disk mid-plane that ensures the presence of a reservoir of neutron rich material (Y_{e}\sim 0.1). This results in the production of neutron rich outflows (\langle Y_{e}\rangle\sim 0.2). The resulting nucleosynthesis yields all r-process elements between the second and the third peak. If the kinetic energy dissipation in the innermost part of the disk results in a significant neutrino luminosity, neutrino absorption increases the electron fraction of the ejecta, producing also elements down to the first peak. Neutrino influence on the properties of the viscous ejecta can be even more relevant in the presence of a long-lived massive neutron star. This is expected to emit a large amount of neutrinos over the diffusion time scale (a few seconds; e.g., Dessart et al., 2009; Perego et al., 2017b). The high neutrino flux is expected to unbind matter in a neutrino-driven wind during the first tens of ms after the merger (Dessart et al., 2009; Perego et al., 2014; Fujibayashi et al., 2017) and can further increase the electron fraction of the viscous ejecta. For a very long lived massive neutron star, the properties of the neutrino- and viscously-driven ejecta could become similar and the resulting r-process nucleosynthesis in the viscous ejecta could be limited to the first and second r-process peaks (Lippuner et al., 2017).

4.3. Effect of the Thermodynamic History of the Ejecta

In principle, nucleosynthesis in the ejecta should be calculated by following the non-equilibrium evolution of the materials composition as it is advected along with the fluid flow and potentially undergoes mixing. Such an approach would require tracking the large number of isotopic abundances needed to follow the r-process flow along with adding stiff, coupled source terms to the composition equations. Such an approach is too computationally expensive, but it would include the possible feedback on the ejecta dynamics of nuclear heating and composition based changes to the EOS (Metzger et al., 2010). The next level of approximation to the nucleosynthesis in the outflows would be following nucleosynthesis in Lagrangian tracers of the flow, while ignoring the backreaction of nucleosynthesis on the flow dynamics. This is likely to be a very reasonable approximation, since the nuclear energy release is only likely to be important to the dynamics of a small amount of marginally bound material. Nevertheless, nuclear burning will produce entropy in these fluid elements which may change the nuclear flow. Therefore, most calculations of nucleosynthesis in binary NS merger ejecta involve taking density histories, \rho(t), of Lagrangian tracers and evolving the composition and entropy of the material in time starting at t_{0}, with entropy s_{0} and electron fraction Y_{e,0} extracted from the simulation output using a self-heating nuclear reaction network (Freiburghaus et al., 1999).

For the nucleosynthesis results discussed above in this section, an even more approximate method has been used in which we assume \rho(t)\approx\rho_{0}(t;\tau_{d},s_{R},Y_{e,R}), where Y_{e,R} and s_{R} are the electron fraction and entropy at of an ejected fluid element when it crosses a radius R and \tau_{d}=eR/(3v_{R})(\rho_{R}/\rho_{0})^{1/3} as described above. This allows us to rapidly calculate nucleosynthesis without including Lagrangian tracer particles in the simulations, but this approximation needs to be tested. Therefore, we have run one simulation with a large number of tracer particles and calculated nucleosynthesis in these tracers using their actual density histories. The results of this calculation compared to our approximate method of calculating nucleosynthesis are shown in Fig. 21. At all mass numbers, the calculations agree to within a factor of two and in most regions to better than 20%. This agreement is good enough to give us confidence in our approximation. Nevertheless, we would like to understand where the variations come from.

Figure 21.— Integrated nucleosnythesis from a simulation using the actual density history from the simulation extracted from tracer particles (dashed line) and using the integrated mass flux approximation described in the text (solid line). The thin dashed line shows the ratio of the two calculations.

The first possible source of error in our approximation is that going from \rho(t)\rightarrow\rho_{0}(t;\tau_{d},s_{R},Y_{e,R}) introduces errors into the calculation. We test this for individual tracer particles by running nucleosynthesis calculations using both the actual density history \rho(t) and using \rho_{0}(t). The nucleosynthesis results for these particles are shown in Fig. 22. Our approximate functional form for the density captures most of the behavior of the actual particles and we recover very similar nucleosynthesis in both cases. All abundances above 10^{-5} agree within a factor of two between the two calculations for all Y_{e}.

Figure 22.— Nucleosynthesis for four selected tracer particles calculated using the actual density history from the simulation (solid lines) and using a parameterized density history extracted as described in the text (dashed lines). The legend labels the Y_{e} of each separate trajectory. We see that there is generally good agreement between the nucleosynthesis predictions with the different density histories across a broad range of Y_{e}. The abundance curves are scaled by arbitrary factors for clarity.

The second possible source of error comes from sampling error in the Lagrangian tracer particles. To test how well they represent the underlying Eulerian flow, we compare the distribution of Y_{e,R} and s_{R} inferred from integrating the flux using the Eulerian outflow and sampling the conditions in the Lagrangian particles at R. These are shown in Fig. 23. The Y_{e} and entropy of the particle may change between R=300\,G/c^{2}\,M_{\odot} and the time nucleosynthesis begins at {\sim}6\,\textrm{GK}. The average electron for the Eulerian outflow is \langle Y_{e}\rangle=0.22, while it is \langle Y_{e}\rangle=0.20 at both measurement points in the tracer particles. By comparing the tracer Y_{e} distributions at these two points, we can see that this approximation introduces at most a 25\% error in the Y_{e} distribution (at Y_{e}=0.04) and substantially less error at most Y_{e}. The difference between either tracer Y_{e} distribution and the integrated Eulerian outflow Y_{e} distribution is substantially larger, almost 40\% at Y_{e}=0.25. We have checked that this error is not due to undersampling by the Lagrangian tracer particles.

Figure 23.— Distribution of the electron fraction in the ejecta as inferred from the tracer particles when they reach a radius of 300\,G/c^{2}\,M_{\odot}, when they reach a temperature of 6\,\textrm{GK}, and from the integrated mass outflow at 300\,G/c^{2}\,M_{\odot} calculated on the Eulerian grid as described above.

5. Electromagnetic Signatures

The radioactive decay of freshly synthesized r-process nuclei in the expanding ejecta powers a quasi-thermal emission known as kilonova (Li & Paczynski, 1998; Kulkarni, 2005; Metzger et al., 2010). The properties of the emission crucially depend on the amount of mass, on the expansion velocity, and on the detailed composition of the ejecta. The latter, in particular, determines the photon opacity, which can substantially differ from the opacity of iron group elements in the presence of a significant fraction of lanthanides (e.g., Kasen et al., 2013; Tanaka & Hotokezaka, 2013).

Eleven hours after the detection of GW170817, a kilonova transient consistent with a binary NS merger was detected and intensively followed up for a few weeks (AT2017gfo; Abbott et al., 2017e; Arcavi et al., 2017; Coulter et al., 2017; Drout et al., 2017; Evans et al., 2017; Kasliwal et al., 2017; Nicholl et al., 2017; Smartt et al., 2017; Soares-Santos et al., 2017; Tanvir et al., 2017). The emission peaked in the UV and visible bands within the first day. Subsequently, the kilonova reddened and peaked in the near infrared (NIR) on a time scale of a few days.

A long-lasting synchrotron remnant is suggested as another promising electromagnetic counterparts to neutron star mergers (Nakar & Piran, 2011; Hotokezaka et al., 2016). Previous studies showed that such signals can be produced by several different ejecta components, e.g., short GRB jet, cocoon, and dynamical ejecta (van Eerten & MacFadyen, 2011; Piran et al., 2013; Hotokezaka & Piran, 2015). A synchrotron remnant has been detected in X-ray, optical, and radio bands for GW170817 (Haggard et al., 2017; Hallinan et al., 2017; Margutti et al., 2017; Mooley et al., 2017; Ruan et al., 2018; Troja et al., 2017; Alexander et al., 2018; Dobie et al., 2018; Lyman et al., 2018; Margutti et al., 2018). Recently, Mooley et al. (2018) and Ghirlanda et al. (2018) observed the superluminal motion of a compact radio emission region in GW170817 and provided us with direct evidence of the existence of a narrowly collimated off-axis jet.

5.1. Kilonovae

We compute synthetic kilonova light curves for the non-viscous, standard resolution models presented in Table 2 for which disk masses are available. We use the semi-analytical kilonova model presented in Perego et al. (2017a). The model accounts for different ejecta components and their possible anisotropies. It assumes azimuthal symmetry around the rotational axis of the binary and reflection symmetry with respect to the orbital plane. The polar angle \theta is discretized in 30 angular bins equally spaced in \cos{\theta}. Within each polar ray, the luminosity, radius, and temperature at the outer photosphere are related through Stefan-Boltzmann’s law. Finally, we compute AB magnitudes in different UV/visible/IR bands assuming the distance of the source to be 40 Mpc.

Following Villar et al. (2017), we include a floor temperature T_{\rm c} for the photosphere. We assume all material with temperature less than T_{\rm c} to be transparent so that the effective photosphere temperature is always larger or equal to T_{\rm c} (see, e.g., Barnes & Kasen, 2013). Guided by the values of T_{\rm c} obtained by Villar et al. (2017) for AT2017gfo for different ejecta components, we assume a composition dependent T_{\rm c}. For material having Y_{e}>0.25 we set T_{\rm c}=T_{\rm c,blue}=3000\,{\rm K}, otherwise we set T_{\rm c}=T_{\rm c,red}=1000\,{\rm K}.

Angular profile, expansion velocity, and electron fraction of the dynamical ejecta are directly extracted from the simulations. Where Y_{e} is larger than 0.25, the kilonova model assumes a Lanthanide-free gray photon opacity \kappa=\kappa_{\rm blue}=1.0~{}{\rm cm^{2}~{}g^{-1}}. Otherwise we set \kappa=\kappa_{\rm red}=30~{}{\rm cm^{2}~{}g^{-1}}.

For the secular ejecta we include both neutrino-driven and viscously-driven winds. The mass of these two components are assumed to be fixed fractions of the remnant accretion disk, as outlined in Sec. 3.2.

The neutrino-driven wind part of the secular ejecta is assumed to be uniformly distributed from the pole (\theta=0) to \theta=60^{\circ}. This ejecta component is assumed to expand radially with v_{\rm rms}=0.08~{}c. We assume low opacity \kappa=\kappa_{\rm blue} for the part of the neutrino driven wind at \theta\leq 45^{\circ}, which is found to be less neutron rich in postmerger simulations (e.g., Perego et al., 2014; Martin et al., 2015), while we set the opacity of the rest of the neutrino-driven wind to be \kappa=\kappa_{\rm purple}=5.0\,{\rm cm^{2}\,g^{-1}}, since this part of the wind is expected to have a broader distribution of Y_{e}’s.

The viscously-driven wind is assumed to have a \sin^{2}{\theta} distribution in mass as a function of the polar angle and to expand homologously with v_{\rm rms}=0.06~{}c. Due to the expected broad and homogeneous Y_{e} distribution of this outflow component (e.g., Metzger & Fernández, 2014; Fujibayashi et al., 2018), we assume \kappa=\kappa_{\rm purple} at all angles.

The model parameters discussed above, as well as other parameters not explicitly mentioned, have been calibrated using AT2017gfo in Perego et al. (2017a). We refer to that study for a more detailed account of the procedure used to generate synthetic kilonova lightcurves.

Figure 24.— Kilonova synthetic light curves in three different bands (g, z, and K_{s}) for three representative models. The left panel shows a binary with prompt BH formation. The middle panel shows a binary forming an hypermassive massive NS (t_{\rm BH}=7.7~{}{\rm ms}). Finally, the right panel shows a binary forming a long-lived supramassive NS (t_{\rm BH}>34.2~{}{\rm ms}). Solid and dashed lines correspond to polar and equatorial viewing angles, respectively. Prompt BH formation binaries do not form massive accretion disk and result in faint and fast transients. Longer-lived NS remnants are associated with the formation of more massive disks that are source of more abundant secular outflows. They result in brighter kilonovae that evolve on longer timescales.

In Fig. 24 we present synthetic light curves for three representative binaries in three different bands. The left panel shows a prompt BH formation case, SFHo_M144139_LK, the middle panel shows a short-lived HMNS case, SFHo_M135135_M0, and the right panel shows a case where a black hole has not formed by the end of the simulation, DD2_M140120_M0. In the latter case we assume that the massive NS will not collapse within the accretion disk lifetime. In a previous work, we have already speculated on the possible implications of a long-lived or even stable massive NS on the kilonova light curves (Radice et al., 2018a). We compute light curves in a band in the visible part of the EM spectrum (g) and in two NIR bands (z and K_{s}). These are chosen because their effective wavelength midpoints satisfy the proportion 1:2:4, and because they span a significant fraction of expected detection range.

The SFHo_M144139_LK binary undergoes prompt collapse and forms a BH surrounded by a light accretion disk with M_{\rm disk}=9\times 10^{-4}~{}M_{\odot}. In this case the ejecta is predominantly of dynamical origin (M_{\rm ej}=4\times 10^{-4}M_{\odot}). The outflow is neutron rich (\langle Y_{e}\rangle=0.18) and expands rapidly (v_{\rm ej}=0.33\,c). Despite the large effective photon opacity of the ejecta, because of the small mass and fast expansion of the outflows, this binary produces a rapid, but faint transient peaking within the first day of the merger in all bands. Polar observers preferentially receive radiation from lower opacity, but also lower density material. For these observers the kilonova fades even more rapidly and is bluer.

The SFHo_M135135_M0 binary produces a short-lived HMNS. After its collapse a BH-torus system with a disk of 0.0123~{}M_{\odot} is formed. We estimate that this binary will eject \lesssim 3\times 10^{-3}~{}M_{\odot} of material in the form of secular ejecta, an amount comparable to that of the dynamical ejecta, {\sim}4.2\times 10^{-3}M_{\odot}. The resulting kilonova transient peaks within a day in both g and z bands and subsequently reddens. Because of the higher Y_{e} of the dynamical ejecta in the polar region, and because of the presence of the neutrino-driven wind from the disk, polar observers will see a marginally brighter kilonova in the UV/visible bands compared to equatorial observers. Conversely, due to the larger amount of dynamical and viscous ejecta close the orbital plane, equatorial observers will receive a larger NIR flux.

The results are qualitatively different for the third binary, DD2_M140120_M0, which forms a long-lived remnant. While the mass and the characteristics of the dynamical ejecta from this binary are in line with that of the SFHo_M135135_M0 binary, the former produces a much more massive disk (M_{\rm disk}\approx 0.19~{}M_{\odot}). Consequently, the secular ejecta dominates over the dynamical ejecta in this case and the outflow is overall more massive and has a lower expansion velocity. The ejecta are optically thick for days and this results in more slowly evolving light curves in all bands. The light curves in both the z and the K_{s} bands present a double peak structure with a first peak around 1 day and a second one at around 5-7 days depending on the inclination. The z band shows a kink instead of a second peak, and both peak and kink are shifted to earlier times compared to the two peaks redder filters. These features are the result of the presence of different ejecta components: the dynamical and neutrino-driven wind ejecta, which mostly determine the early light curve, and the viscously-driven wind, which dominates at late times.

Figure 25.— Color evolution of the kilonova light curves computed as the difference in the AB magnitude between g and K_{s} bands for three representative models. Binaries forming HMNSs or dynamically stable remnants power significantly bluer transients.

The color evolutions for these models, exemplified by the difference in magnitude between the g and K_{s} bands, are shown in Fig. 25. All of the synthetic kilonova signals show a fast evolution towards the infrared after about a day from the merger. However, there are are significant difference, between them. This suggests that the color of the kilonova signal could be used to infer the outcome of the merger. We remark that a similar point was also made by Metzger & Fernández (2014) and Lippuner et al. (2017). However, in our case the difference in the color is purely due to the differences in dynamical ejecta and disk masses, since we do not assume that long-lived remnant produce disk winds with different compositions, as instead argued by those authors.

Figure 26.— Synthetic kilonova light curves for an observer located at 45^{\circ} in the g (left), z (middle) and K_{s} (right) bands for the DD2_M140120_M0 binary computed using different combinations of ejecta components: dynamical, “d”, neutrino-driven wind, “w”, and viscously-driven wind, “s”. The comparison between the light curve obtained including all ejecta componetns (d+w+s) and light curves obtained by neglecting in turn one of the components reveals the different roles of the components in shaping the light curves, at different times and in different bands.

We study the relative impact that each of the ejecta components has on the kilonova light curve in the case of the DD2_M140120_M0 binary. In addition to the “full” light curve, which includes all ejecta components, we generate light curves in which, in turn, one of the ejecta components has been removed. In doing so we assume the viewing angle to be 45^{\circ}. The results are shown in Fig. 26. We remind that our model does not simply add contributions from the different components, but also accounts for irradiation and reprocessing effects. Starting from a few days from the merger the light curves are dominated by the viscous ejecta, especially in the NIR bands. This is not surprising, since this component dominates the overall mass ejection. However, we find that both the neutrino-driven wind and the dynamical ejecta play an important role in determining the visible and NIR light curves in the first days. The dynamical ejecta is most important in the very first day, when the viscous ejecta is still too opaque to contribute significantly. All ejecta components are necessary to reproduce the light curve properties in the UV and visible bands.

Figure 27.— Synthetic kilonova light curves in the g, z and K_{s} bands for the BHBlp_M135135_M0 and BHBlp_M135135_LK binaries, both in the polar (top panel) and in the equatorial (bottom panel) directions. Because neutrino-driven winds from the disk dominate the early light curve, the presence of a neutrino-irradiate component of the dynamical ejecta only yields a modest \lesssim 0.5 mag change of the color light curve in the first few hours after merger.

The presence of a polar component of the dynamical ejecta irradiated by neutrinos could, in principle, impact the properties of the kilonova emission. We study this in Fig. 27 where we compare light curves obtained from two binaries having the same NS masses and EOS, but different neutrino treatments: BHBlp_M135135_LK and BHBlp_M135135_M0. The angular distribution and composition of dynamical ejecta of these binaries are shown in Fig. 4. The inclusion of neutrino re-absorption results in an increase of the magnitude in the UV and visible bands within the first few hours of the merger. However, these differences amount only to less than 0.5 mag, even for an observer located along the polar direction. The reason is that the overall UV/optical signal is actually dominated by the secular neutrino-driven wind, while the dynamical ejecta is less important.

In general, we find that the high-latitude neutrino-irradiated part of the dynamical ejecta has a very modest impact on the kilonova light curve. Instead the UV/optical signal is typically dominated by the secular neutrino-driven winds. There are two reasons for this. First, the secular neutrino-driven wind dominates the overall outflow of high-Y_{e} material at high latitudes. Indeed, only a small fraction of the dynamical ejecta is directed far from the orbital plane and experiences significant neutrino irradiation. Second, because of its fast expansion, the lanthanides-free part of the dynamical ejecta becomes transparent already a few hours after merger and does not contribute significantly to the kilonova light curve afterwards. The situation is slightly different for prompt collapse binaries for which the accretion disk masses are modest. However, in these cases the effect of neutrino irradiation of the dynamical ejecta is expected to be negligible.

Figure 28.— Peak times (top), AB magnitudes (middle), and widths (bottom) in three phomometric bands (g, z, and K_{s}, moving from left to right) of the kilonova synthetic light curves for our fiducial subset of simulations assuming a distance of 40 Mpc for the sources. The peak width is defined as the time interval over which the AB magnitude increases by one across the peak. The dashed black lines correspond to the values (or limits) obtained for AT2017gfo, the EM counterpart to GW170817. The horizontal extension of the lines correspond to the 90% highest posterior density interval 300^{+420}_{-230} for \tilde{\Lambda} obtained assuming low-spin priors by Abbott et al. (2018b). The peak times, magnitudes, and widths of AT2017gfo are obtained from Villar et al. (2017), Coulter et al. (2017), Smartt et al. (2017), and Tanvir et al. (2017). BH formation (small \tilde{\Lambda}) results in faint and rapid transients. Longer-lived remnants (large \tilde{\Lambda}) are associated with brighter luminosity in the UV/visible bands and longer lasting, more slowly evolving NIR signals.

A summary of the most important features of the kilonova light curves in different bands is presented in Fig. 28. For each binary in our fiducial subset of simulations we present the time, the magnitude, and the width of the kilonova peak as a function of \tilde{\Lambda} in three different bands: g, z, and K_{s}. The peak width is defined as the time interval around the peak over which \Delta M=+1. We further distinguish between polar and equatorial observers. Different symbols denote different merger outcomes. There is a clear dependence of the kilonova properties on the nature of the merger remnant and on \tilde{\Lambda}. These trends, already anticipated by the three representative cases of Fig. 24, are a consequence of the dependence of the disk mass on \tilde{\Lambda} (see Fig. 15) and of the dominant role of the secular ejecta in determining the properties of the kilonova. Faint and rapidly decreasing kilonovae indicate the formation of BH via prompt collapse, while bright and long lasting light curves in the NIR bands are signatures of longer-lived remnants.

Figure 28 also reports the values or upper limits inferred for AT2017gfo. The properties of the kilonova are compatible with the formation of a massive NS that survived for at least several milliseconds after merger, as also argued by Shibata et al. (2017a). However, it is important to emphasize that our theoretical models are based on a number of assumptions that still need to be verified with first-principle simulations, most notably the assumption that a fixed fraction of the accretion disk is unbound by winds.

The differences in the kilonova brightnesses have important consequences for their detectability. We consider the limiting magnitudes reported in (Rosswog et al., 2017) for the Large Synoptic Survey Telescope (LSST) in the g and z bands, and for the Visual and Infrared Survey Telescope for Astronomy (VISTA) in the K_{\rm s} band and assume exposure times of 60 seconds. Considering the detection horizon for Advanced LIGO to be of 120-170 Mpc (Abbott et al., 2018c), we estimate that kilonovae counterparts would be detectable for a large fraction of the NS merger GW events. If a long-lived massive NS is formed, then the kilonova would be dectable for all GW events in both the g and z bands up to the horizon distance of Advanced LIGO, and up to a luminosity distance of {\sim}90 Mpc in the K_{\rm s} band. If, on the other hand, a BH is rapidly formed after merger, then the kilonova would only be detectable in the K_{\rm s} band to distances of 30-50 Mpc, depending on the orientation of the binary. Kilonovae associated with prompt BH formation are detactable in the z and g bands in the entire Advanced LIGO and Virgo volume, however their observation will be challenging given their fast decline.

Figure 29.— Color of the kilonova light curves for all the fiducial subset of simulations at three different times, as a function of the binary tidal parameter \tilde{\Lambda}. Light curves are computed for an observer located at a polar angle of 45^{\circ}. The color of the transient at different times shows a significant correlation with \tilde{\Lambda}.

In Fig. 29 we show the difference in magnitude between the g and K_{s} bands for kilonovae at three different epocs and as a function of \tilde{\Lambda}. Our results show that the color evolution and properties previously described for the three representative models are generic. Moreover, there is a clear correlation between g-K_{s} and \tilde{\Lambda} which suggests that it might be possible to constrain the EOS of NS matter with EM observations in different bands.

We want to stress that, while the correlations between kilonova light curve properties and \tilde{\Lambda} document in Figs. 28 and 29 appear robust, their quantitative determination must wait until simulations self-consistently including dynamical and secular ejecta from NS mergers become available. The impact of viscosity on the kilonova lightcurves is discussed in a companion paper (Radice et al., 2018b).

5.2. Synchrotron remnants

We calculate the light curve of the synchrotron radiation arising from the dynamical ejecta with the semi-analytic method of Hotokezaka & Piran (2015). According to this model electrons accelerated in the shock between the ejecta and the ISM emit synchrotron radiation in the amplified magnetic filed. The total flux density is calculated by integrating the radiation flux from each solid angle over the equal-arrival time surface. The ISM number density n is a parameter of the model. The conversion efficiencies of the internal energy of the shock to the energy of the accelerated electrons and amplified magnetic field are assumed to be \epsilon_{e} and \epsilon_{B}, respectively. The initial velocity profile of the ejecta is given by mapping the three dimensional structure of the ejecta to a one dimensional structure, then the ejecta are evolved adiabatically.

Before discussing the numerical results, here we give briefly the scalings of the peak time and peak flux with the relevant physical quantities. The peak time of the light curve can be estimated from the deceleration time of the ejecta (Nakar & Piran, 2011):

\displaystyle t_{\rm dec} \displaystyle\sim \displaystyle 3\cdot 10^{3}\,{\rm day}\,\left(\frac{T_{\rm ej}}{10^{50}{\rm erg% }}\right)^{1/3}
\displaystyle~{}~{}~{}\times\left(\frac{n}{10^{-3}}{\rm cm^{-3}}\right)^{-1/3}% \left(\frac{v_{\rm ej}}{0.3c}\right)^{-5/3},

where E and v are the ejecta kinetic energy and initial velocity222This deceleration time is measured in the merger rest frame. The difference between this time and observer time is about a factor of a few.. The peak flux can be estimated as

\displaystyle F_{\nu} \displaystyle\sim 10\,{\mu}{\rm Jy}\,\left(\frac{T_{\rm ej}}{10^{50}{\rm erg}}% \right)\left(\frac{n}{10^{-3}{\rm cm^{-3}}}\right)^{\frac{p+2}{4}}\left(\frac{% \epsilon_{B}}{10^{-2}}\right)^{\frac{p+2}{4}}
\displaystyle\times\left(\frac{\epsilon_{e}}{10^{-1}}\right)^{p-1}\left(\frac{% v_{\rm ej}}{0.3c}\right)^{\frac{5p-7}{2}}\left(\frac{D}{40\,{\rm Mpc}}\right)^% {-2}\left(\frac{\nu}{3\,{\rm GHz}}\right)^{-\frac{p-1}{2}},

where p is the spectral index of the non-thermal electrons.

Figure 30.— Radio light curves of the dynamical ejecta of SFHo_M135135_LK at 3 GHz. Here we assume the microphysics parameters to be \epsilon_{e}=0.1, \epsilon_{B}=0.01, and p=2.16. Also shown as open circles are the observed flux densities at 3 GHz of the afterglow in GW170817 (Hallinan et al., 2017; Mooley et al., 2017, 2018).

As in the previous work by Hotokezaka et al. (2018b), who used the results of high resolution merger simulations performed by Kiuchi et al. (2017), we also find that it is the fast component of the dynamical ejecta with velocities of {\sim}0.3{-}0.8\,c to predominantly produce a synchrotron signal. Figure 30 shows the expected radio signals of the dynamical ejecta of SFHo_M135135_LK compared with GW170817. Here we calculate the light curve with a similar method to Hotokezaka et al. (2018b) and choose the number density of the ISM to be 10^{-4}5\cdot 10^{-3}\,{\rm cm^{-3}} as suggested by Mooley et al. (2018). While the dynamical ejecta component is fainter than the observed flux densities until {\sim}300 days, this component can be detectable in radio and X-rays on time scales of a year to ten years in the optimistic case. The peak flux density depends also on the EOS and it is fainter for the cases of DD2, BHB\Lambda\phi, and LS200 because the kinetic energy in the high velocity component is lower than that of SFHo (see Eq. 5.2 for the scaling).

Figure 31.— Radio light curves of DD2_M135135_LK and BHBlp_M135135_LK at 3 GHz. The distance to the source and ISM density are assumed to be 100 Mpc and 0.1\,{\rm cm^{-3}}, respectively. Here the microphysics parameters are set to be \epsilon_{e}=\epsilon_{B}=0.1 and p=2.5.

The synchrotron remnant arising from the dynamical ejecta may be detectable in future GW events if the ISM density is large enough. For instance, radio counterparts with a flux of \gtrsim 100\,{\rm\mu Jy} can be detectable by blind survey with VLA, ASKAP, and MeerKAT when the GW localization area is better than {\sim}30\,{\rm deg^{2}} (Hotokezaka et al., 2016). Of course, if the host galaxy is identified by finding other EM counterparts, the detection limit is reduced to a few tens of {\rm\mu Jy}. Figure 31 shows the expected light curves for DD2_M135135_LK and BHBlp_M135135_LK at 3 GHz, assuming a distance of 100 Mpc and density of 0.1\,{\rm cm^{-3}}. Note that the radio flux of SFHo_M135135_LK is brighter than both models under the same condition (density, distance and microphysics parameters). Therefore, the radio afterglow arising from the dynamical ejecta can be detectable for SFHo_M135135_LK and BHBlp_M135135_LK if a merger occurs within 100 Mpc and the surrounding ISM density of \gtrsim 0.1\,{\rm cm^{-3}}.

It is worth noting that differences in the high density part of the EOSs lead to large differences in the expected light curves. For example, as shown in Fig. 31, BHBlp_M135135_LK results in the radio flux that is much brighter than that of DD2_M135135_LK, even though these EOSs share the same form up to around the nuclear saturation density. The reason is that BHBlp_M135135_LK produces a larger amount of fast-moving ejecta because of the more violent merger dynamics due to the appearance of \Lambda hyperons at high densities (see Fig. 9). Therefore, we may be able to constrain the neutron star EOS using the observed light curves of the synchrotron remnants of future GW events.

The impact of viscosity on the synchrotron lightcurves is discussed in a companion paper (Radice et al., 2018b).

6. Conclusions

We have performed 59 full-GR NS merger simulations employing a microphysical treatment of the NS matter and including compositional and energy changes due to the emission of neutrinos. A subset of our simulations also included a treatment of neutrino re-absorption and/or of the effective viscosity due to MHD turbulence in the stars. This is the largest set of NS simulations with realistic microphysics to date. Our studies focused on the ejection of material during and after the mergers, and on the associated nucleosynthetic and EM signatures.

6.1. Mass Ejection

We find that material is ejected on a dynamical timescale during the mergers due to the combined effects of tidal torques, shocks, and neutrino heating. Tidally-driven ejecta flow close to the orbital plane of the binary, have low entropies ({\sim}10\,k_{\rm B} per baryon), and are very neutron rich, with electron fraction {\sim}0.1. Shock-driven ejecta have broad distributions in entropy and electron fractions, and are distributed over a broad {\sim}60^{\circ} angle from the orbital plane. Neutrino driven winds are emitted preferentially close to the polar axis and have electron fractions in excess of 0.25. In accordance with previous studies using approximate GR and/or approximate microphysics (e.g., Hotokezaka et al., 2013b; Bauswein et al., 2013), we find that the most conspicuous episode of mass ejection is triggered by the centrifugal bounce of the merger remnant shortly after merger. Overall, we find dynamical ejecta masses of few times 10^{-3}\,M_{\odot} with average electron fractions \lesssim 0.25, and average entropies in the range 10{-}30\,k_{\rm B} per baryon.

The mass of the dynamical ejecta has a large numerical uncertainty. From the comparison of low- and high-resolution data, we estimate the relative error in its determination to be as large as 50%. Even larger discrepancies of a factor of a few are found when comparing with and among results published by different groups (Bauswein et al., 2013; Sekiguchi et al., 2015; Lehner et al., 2016; Bovard et al., 2017). On the other hand, the intensive properties of the ejecta, e.g., asymptotic velocity, composition, etc., appear to be robust with resolution.

We have studied the dependency of the outflow properties on the NS masses and EOS. Following Dietrich & Ujevic (2017) we have constructed fits of the ejecta mass and velocity in terms of the NS masses and compactnesses. These capture some of the qualitative trend of our data reasonably well, especially for the ejecta velocity. However, the fitting coefficients we infer from our data are significantly different from those of Dietrich & Ujevic (2017). This is not too surprising given that most of the simulations used by Dietrich & Ujevic (2017) employed idealized microphysics and neglected weak reactions, which likely resulted in the overestimation of the ejecta mass (Radice et al., 2016b). We also find that there are systematic effects not captured by these fits. One of the main reasons is that the strength of the bounce of the merger remnant, an important parameter determining the ejecta mass, depends on details of the EOS at higher densities and temperatures than those determining the fitting parameters. Accurate models of the dynamical ejecta mass could have important applications to multimessenger astronomy and, indeed, the fits of Dietrich & Ujevic (2017) have already been applied to GW170817 (Abbott et al., 2017b; Coughlin et al., 2018). However, our results indicate that more work and better simulations will be required to build truly quantitative ejecta models.

Our results indicate that soft NS EOSs predict larger ejecta masses, at least in the range of mass ratios we have probed, in accordance with previous results (Hotokezaka et al., 2013b; Bauswein et al., 2013; Lehner et al., 2016; Sekiguchi et al., 2016; Dietrich et al., 2017b; Bovard et al., 2017). EOS effects are also imprinted in the ratio between the masses of the tidal and shock-heated components of the ejecta. However, we do not find any clear correlation between ejecta mass and velocity and the tidal deformability of the binary \tilde{\Lambda}. Consequently, it appears to be impossible to directly constrain the NS EOS from measurements of the dynamical ejecta. On the other hand, dynamical ejecta measurements might help to break degeneracies in the other binary parameters and, in this way, they might improve the constraints on the tidal deformability of NSs indirectly.

We have analyzed the geometry of the outflows and found that it depends on the relative amount of tidal- and shock-driven ejecta, as does the average electron fraction. Consequently, we find that there is a correlation between the rms opening angle of the ejecta and their average electron fraction. This suggests that it might be possible to indirectly constrain the composition, and hence the nucleosynthetic yields, of the dynamical ejecta by combining observations of similar binary systems with different orientations.

We have computed the velocity distribution of the dynamical ejecta and found that shocks during and shortly after merger can accelerate a small fraction of the ejecta {\sim}10^{-6}\,M_{\odot} to asymptotic velocities in excess to 0.6\,c. This fast moving material is preferentially located close to the orbital plane, presumably because of the oblate shape of the merger debris cloud into which the accelerating shocks propagates. It is also distributed in a very anisotropic fashion. The amount of material reaching these high velocities depends on the binary parameters and on the EOS. Binaries with compact NSs typically produce significantly more of this fast-moving component of the outflow. When comparing binaries simulated with the BHB\Lambda\phi and DD2 EOSs, we find that the former systematically predicts a larger mass of the fast moving ejecta by a factor of a few. This is because, even though BHB\Lambda\phi and DD2 predict identical NS structures for most of the binaries we have considered, the former softens due to the appearance of \Lambda-hyperons after merger and predicts more violent bounces (Radice et al., 2017).

We find that binaries with remnants that are stable for at least several rotation periods after merger result in the formation of massive {\sim}0.1{-}0.2\,M_{\odot} accretion disks. Conversely, BHs formed promptly after merger, i.e., within {\sim}1\,{\rm ms}, are endowed with rather light accretion disks {\sim}10^{-3}\,M_{\odot}. More in general, the accretion disk masses are found to depend on the lifetime of the remnant and to correlate with the binary’s tidal deformability (Radice et al., 2018c). It is expected that neutrino- and viscously-driven outflows will carry away a significant fraction {\sim}10{-}50\% of these accretion disks on a timescale of a few seconds after the end of our calculations (Metzger & Fernández, 2014; Just et al., 2015; Siegel & Metzger, 2018; Fujibayashi et al., 2018; Fernández et al., 2018). Accordingly, these secular outflows are expected to be the dominant component of the ejecta. For this reason, kilonova observations can be used to constrain the accretion disk masses (Perego et al., 2017a) and binary tidal deformabilities (Radice et al., 2018c).

6.2. Nucleosynthesis

We have computed the nucleosynthetic yields of the dynamical ejecta. We find that second and third r-process peak isotopes, i.e., with A\gtrsim 125, are robustly produced with relative abundance close to solar. However, the relative abundance of light r-process elements, i.e., with 90\lesssim A\lesssim 125, and second and third peak isotopes are sensitive to the binary masses, weak reactions, and, to a lesser extend, to the NS EOS.

We find that, for binaries close to the prompt BH formation threshold, different EOS predict different ratios of tidal- and shock-driven dynamical ejecta and, consequently, different relative abundances of light r-process elements and second and third peak r-process elements. In particular, the LS220 and DD2 EOSs predicts that high mass binaries should produce a smaller relative fraction of light elements compared to low mass binaries, while the BHB\Lambda\phi and SFHo EOSs have the opposite trend.

The absorption of neutrinos has a large impact on the composition and yields of the outflow. When neutrino absorption is included, the ejecta distribution in electron fraction become broad and Y_{e}’s of up to 0.4 are reached. The dynamical ejecta from comparable mass binaries in simulations that included neutrino re-absorption produce r-process isotopic abundances close to solar. If neutrino re-absorption is not included, instead, light r-process elements are underproduced in our simulations. The isotopic abundances in the region of the first r-process peak are also sensitive to details of the neutrino treatment, such as the assumptions made for the incoming neutrino energy. This points to the need for simulations including energy dependent neutrino-radiation treatments. This will be the object of future work.

The binary mass ratio also impacts the production of first-peak elements, with unequal mass systems underproducing elements with 90\lesssim A\lesssim 125. We find the relative yields to depend on the mass ratio almost as sensitively as on the neutrino treatment.

The secular ejecta is expected to be an important, if not dominant, component of the overall outflow, especially in the case of long-lived remnants or when massive disks are formed. Consequently, the secular ejecta should be taken into account when estimating the nucleosynthetic yields from NS mergers. Presently, however, we can only speculate on the nucleosynthesis from this component. Long term merger and postmerger simulations will be required to construct self-consistent models. We leave this to future work.

We have estimated the r-process nucleosynthesis yields using the simulation data extracted at a fixed radius of 300\,G/c^{2}\,M_{\odot}\simeq 450\,{\rm km} and used precomputed yields from parametrized trajectories. This is the same strategy we employed in Radice et al. (2016b). As part of our analysis, we have now validated this procedure by comparing with a more computationally expensive analysis done with Lagrangian tracer particles. We find the nucleosynthesis to be rather insensitive to the detailed thermodynamical history of the tracer particles. The electron fraction of the material does not significantly evolve between the time we record it at 300\,G/c^{2}\,M_{\odot} and when the temperature drops below 6\,{\rm GK} and the r-process begins. On the other hand, we find that, due to advection errors in the position of the tracer particles, the Y_{e} distribution inferred from the tracers can differ quantitatively from that inferred from the grid data at the same radius. Overall, however, the differences between the results obtained with the two methods are well below those associated with other sources of uncertainty, most notably the treatment of neutrino radiation.

6.3. Electromagnetic Counterparts

We have computed synthetic kilonova light curves for all our models using the model of Perego et al. (2017a). Our light curve model is informed using the detailed angular structure of the outflows, density, velocity, and composition, as extracted from the simulations. These have been augmented with the inclusion of secular ejecta composed of neutrino-driven and viscously-driven winds which we assume to entrain 3% and 20% of the disk, respectively. These values are motivated by recent long-term postmerger simulations (Perego et al., 2014; Martin et al., 2015; Just et al., 2015; Metzger & Fernández, 2014; Just et al., 2015; Siegel & Metzger, 2018; Fujibayashi et al., 2017, 2018; Fernández et al., 2018).

Binaries resulting in prompt BH formation have kilonovae dominated by the dynamical ejecta. These peak at about a day after merger and then rapidly decline. They are also very red, with g-K_{s}\simeq 10 mag three days after the merger. As the merger remnant lifetime increases, so does the disk mass, and hence the associated kilonovae become increasingly dominated by radiation coming from the secular ejecta. These kilonovae peak on longer timescales of few days to a week and are bluer, with g-K_{s}\simeq 3{-}7 mag three days after the merger. This suggests that kilonovae observations could be used to probe the lifetime of the merger remnant. A similar suggestion was also put forward by Metzger & Fernández (2014) and, more recently, by Lippuner et al. (2017). They argued that longer lived remnants should result in bluer kilonovae because of the neutrino irradiation of the ejecta from the central remnant. This is expected to lower the electron fraction of the outflow below the threshold needed for the production of lanthanides resulting in a drop of the photon opacity of the material. The phenomenon we find goes in the same direction, but is physically distinct since it arises from the correlation between disk masses and remnant lifetimes and not from changes in the ejecta photon opacity, which we have kept constant instead.

Since for most binaries the kilonova signal is dominated by the secular ejecta, the effect of the inclusion of neutrino re-absorption in the modeling of the dynamical ejecta has only a modest impact on the kilonova signal. Moreover, since remnants disk masses strongly correlate with binary tidal deformabilities, the peak properties of the kilonovae, peak time, magnitude, width, are found to depend sensitively on \tilde{\Lambda}. This suggests that kilonovae observations could be used directly to probe the EOS of NSs. However, it is important to stress than this correlation is in part due to our assumption that a fixed fraction of the accretion disk is ejected after merger. Long-term simulations of the evolution of postmerger remnants are necessary to verify if this assumption is valid, or whether light and massive postmerger disks evolve in quantitatively different ways.

We have computed the synchrotron radio signal expected from the interaction of the ejecta with the ISM using the model of Hotokezaka & Piran (2015). The radio fluence and the timescales over which the radio remnant evolves depend sensitively on the ISM density, and on the kinetic energy and the detailed velocity distribution of the ejecta.

We find that some of our models predict a rebrightening of the synchrotron signal from GW170817 on a timescale of months to years from the merger, when the emission from the ejecta will start to dominate over the emission coming from the short \gamma-ray burst (SGRB) jet.

Depending on the ISM density and the orientation of the binary, it could be possible to detect the radio signal due to the interaction of the ejecta with the ISM on timescales of weeks to months after merger at distances of up to {\sim}100 Mpc. These observations would allow us to probe the velocity distribution of the ejecta and, hence, the violence with which the merged objects bounces back after the NS collision. The latter, in turn, depends on the EOS of matter at several times nuclear density. For instance, due to the appearance of \Lambda hyperons after the merger, the BHB\Lambda\phi EOS predicts more violent mergers and brighter radio flares by up to two orders of magnitude compared with the DD2 EOS from which it only differs at high densities due to the inclusion of hyperons. This suggests that radio flares could be used to probe the EOS of NSs in a regime that is not accessible with GW observations of the inspiral alone333The high density EOS could also be constrained with GW observations of the postmerger signal (Bernuzzi et al., 2015; Radice et al., 2017; Most et al., 2018a; Bauswein et al., 2018)..

It is a pleasure to acknowledge A. Burrows for the many stimulating discussions. DR gratefully acknowledges support from the Schmidt Fellowship and the Max-Planck/Princeton Center (MPPC) for Plasma Physics (NSF PHY-1523261). AP acknowledges support from the INFN initiative ”High Performance data Network” funded by CIPE. DR and AP acknowledge support from the Institute for Nuclear Theory (17-2b program). AP thanks the Institute for Advanced Study for its hospitality and support. SAF acknowledges support from the United States Department of Energy through the Computational Science Graduate Fellowship. SB acknowledge support by the EU H2020 under ERC Starting Grant, no. BinGraSp-714626. LFR acknowledges support from U.S. Department of Energy through the award number DE-SC0017955. The simulations were performed on BlueWaters, Bridges, Comet, and Stampede, and were enabled by the NSF PRAC program (ACI-1440083 and AWD-1811236) and the NSF XSEDE program (TG-PHY160025). The analysis employed computational resources provided by both the TIGRESS high performance computer center at Princeton University, which is jointly supported by the Princeton Institute for Computational Science and Engineering (PICSciE) and the Princeton University Office of Information Technology, and the Institute for Cyber-Enabled Research, which is supported by Michigan State University.


  • Abbott et al. (2017a) Abbott, B. P., et al. 2017a, Nature, 1710.05835
  • Abbott et al. (2017b) ——. 2017b, Astrophys. J., 850, L39, 1710.05836
  • Abbott et al. (2017c) ——. 2017c, Astrophys. J., 848, L13, 1710.05834
  • Abbott et al. (2017d) ——. 2017d, Phys. Rev. Lett., 119, 161101, 1710.05832
  • Abbott et al. (2017e) ——. 2017e, Astrophys. J., 848, L12, 1710.05833
  • Abbott et al. (2018a) ——. 2018a, 1805.11581
  • Abbott et al. (2018b) ——. 2018b, 1805.11579
  • Abbott et al. (2018c) ——. 2018c, Living Rev. Rel., 21, 3, 1304.0670, [Living Rev. Rel.19,1(2016)]
  • Alexander et al. (2018) Alexander, K. D., et al. 2018, Astrophys. J., 863, L18, 1805.02870
  • Annala et al. (2018) Annala, E., Gorda, T., Kurkela, A., & Vuorinen, A. 2018, Phys. Rev. Lett., 120, 172703, 1711.02644
  • Arcavi et al. (2017) Arcavi, I. et al. 2017, The Astrophysical Journal, 848, L33
  • Ardevol-Pulpillo et al. (2018) Ardevol-Pulpillo, R., Janka, H. T., Just, O., & Bauswein, A. 2018, 1808.00006
  • Arlandini et al. (1999) Arlandini, C., Kaeppeler, F., Wisshak, K., Gallino, R., Lugaro, M., Busso, M., & Straniero, O. 1999, Astrophys. J., 525, 886, astro-ph/9906266
  • Baiotti et al. (2008) Baiotti, L., Giacomazzo, B., & Rezzolla, L. 2008, Phys. Rev., D78, 084033, 0804.0594
  • Baker et al. (2017) Baker, T., Bellini, E., Ferreira, P. G., Lagos, M., Noller, J., & Sawicki, I. 2017, Phys. Rev. Lett., 119, 251301, 1710.06394
  • Banik et al. (2014) Banik, S., Hempel, M., & Bandyopadhyay, D. 2014, Astrophys. J. Suppl., 214, 22, 1404.6173
  • Barkov et al. (2018) Barkov, M. V., Kathirgamaraju, A., Luo, Y., Lyutikov, M., & Giannios, D. 2018, 1805.08338
  • Barnes & Kasen (2013) Barnes, J., & Kasen, D. 2013, Astrophys. J., 775, 18, 1303.5787
  • Baumgarte et al. (2000) Baumgarte, T. W., Shapiro, S. L., & Shibata, M. 2000, Astrophys. J., 528, L29, astro-ph/9910565
  • Bauswein et al. (2018) Bauswein, A., Bastian, N.-U. F., Blaschke, D. B., Chatziioannou, K., Clark, J. A., Fischer, T., & Oertel, M. 2018, 1809.01116
  • Bauswein et al. (2013) Bauswein, A., Goriely, S., & Janka, H. T. 2013, Astrophys. J., 773, 78, 1302.6530
  • Bauswein et al. (2017) Bauswein, A., Just, O., Janka, H.-T., & Stergioulas, N. 2017, Astrophys. J., 850, L34, 1710.06843
  • Beniamini et al. (2018) Beniamini, P., Petropoulou, M., Duran, R. B., & Giannios, D. 2018, 1808.04831
  • Berger & Colella (1989) Berger, M., & Colella, P. 1989, J. Comput. Phys., 82, 64
  • Berger & Oliger (1984) Berger, M. J., & Oliger, J. 1984, J. Comput. Phys., 53, 484
  • Bernuzzi et al. (2015) Bernuzzi, S., Dietrich, T., & Nagar, A. 2015, Phys. Rev. Lett., 115, 091101, 1504.01764
  • Bernuzzi & Hilditch (2010) Bernuzzi, S., & Hilditch, D. 2010, Phys. Rev., D81, 084003, 0912.2920
  • Bernuzzi et al. (2016) Bernuzzi, S., Radice, D., Ott, C. D., Roberts, L. F., Moesta, P., & Galeazzi, F. 2016, Phys. Rev., D94, 024023, 1512.06397
  • Bovard et al. (2017) Bovard, L., Martin, D., Guercilena, F., Arcones, A., Rezzolla, L., & Korobkin, O. 2017, Phys. Rev., D96, 124005, 1709.09630
  • Bruenn (1985) Bruenn, S. W. 1985, Astrophys. J. Suppl., 58, 771
  • Burrows et al. (2006) Burrows, A., Reddy, S., & Thompson, T. A. 2006, Nucl. Phys., A777, 356, astro-ph/0404432
  • Chornock et al. (2017) Chornock, R., et al. 2017, Astrophys. J., 848, L19, 1710.05454
  • Ciolfi et al. (2017) Ciolfi, R., Kastaun, W., Giacomazzo, B., Endrizzi, A., Siegel, D. M., & Perna, R. 2017, Phys. Rev., D95, 063016, 1701.08738
  • Coughlin et al. (2018) Coughlin, M. W., et al. 2018, Mon. Not. Roy. Astron. Soc., 1805.09371
  • Coulter et al. (2017) Coulter, D. A., et al. 2017, Science, 1710.05452
  • Cowperthwaite et al. (2017) Cowperthwaite, P. S., et al. 2017, Astrophys. J., 848, L17, 1710.05840
  • De et al. (2018) De, S., Finstad, D., Lattimer, J. M., Brown, D. A., Berger, E., & Biwer, C. M. 2018, Phys. Rev. Lett., 121, 091102, 1804.08583
  • Deaton et al. (2018) Deaton, M. B., O’Connor, E., Zhu, Y. L., Bohn, A., Jesse, J., Foucart, F., Duez, M. D., & McLaughlin, G. C. 2018, 1806.10255
  • Dessart et al. (2009) Dessart, L., Ott, C., Burrows, A., Rosswog, S., & Livne, E. 2009, Astrophys. J., 690, 1681, 0806.4380
  • Dietrich et al. (2015) Dietrich, T., Bernuzzi, S., Ujevic, M., & Brügmann, B. 2015, Phys. Rev., D91, 124041, 1504.01266
  • Dietrich et al. (2017a) Dietrich, T., Bernuzzi, S., Ujevic, M., & Tichy, W. 2017a, Phys. Rev., D95, 044045, 1611.07367
  • Dietrich et al. (2018) Dietrich, T. et al. 2018, 1806.01625
  • Dietrich & Ujevic (2017) Dietrich, T., & Ujevic, M. 2017, Class. Quant. Grav., 34, 105014, 1612.03665
  • Dietrich et al. (2017b) Dietrich, T., Ujevic, M., Tichy, W., Bernuzzi, S., & Bruegmann, B. 2017b, Phys. Rev., D95, 024029, 1607.06636
  • Dobie et al. (2018) Dobie, D. et al. 2018, ApJ, 858, L15, 1803.06853
  • Drago et al. (2018) Drago, A., Pagliara, G., Popov, S. B., Traversi, S., & Wiktorowicz, G. 2018, Universe, 4, 50, 1802.02495
  • Drout et al. (2017) Drout, M. R., et al. 2017, Science, 1710.05443
  • East et al. (2016) East, W. E., Paschalidis, V., Pretorius, F., & Shapiro, S. L. 2016, Phys. Rev., D93, 024011, 1511.01093
  • Eichler et al. (1989) Eichler, D., Livio, M., Piran, T., & Schramm, D. N. 1989, Nature, 340, 126
  • Eichler et al. (2015) Eichler, M., et al. 2015, Astrophys. J., 808, 30, 1411.0974
  • Einfeldt (1988) Einfeldt, B. 1988, in Shock tubes and waves; Proceedings of the Sixteenth International Symposium, Aachen, Germany, July 26–31, 1987. VCH Verlag, Weinheim, Germany, 671
  • Evans et al. (2017) Evans, P. A., et al. 2017, Science, 1710.05437
  • Fattoyev et al. (2018) Fattoyev, F. J., Piekarewicz, J., & Horowitz, C. J. 2018, Phys. Rev. Lett., 120, 172702, 1711.06615
  • Favata (2014) Favata, M. 2014, Phys. Rev. Lett., 112, 101101, 1310.8288
  • Fernández & Metzger (2013) Fernández, R., & Metzger, B. D. 2013, Mon. Not. Roy. Astron. Soc., 435, 502, 1304.6720
  • Fernández et al. (2018) Fernández, R., Tchekhovskoy, A., Quataert, E., Foucart, F., & Kasen, D. 2018, 1808.00461
  • Finstad et al. (2018) Finstad, D., De, S., Brown, D. A., Berger, E., & Biwer, C. M. 2018, Astrophys. J., 860, L2, 1804.04179
  • Flanagan & Hinderer (2008) Flanagan, E. E., & Hinderer, T. 2008, Phys. Rev., D77, 021502, 0709.1915
  • Foucart (2012) Foucart, F. 2012, Phys. Rev., D86, 124007, 1207.6304
  • Foucart et al. (2018a) Foucart, F., Duez, M. D., Kidder, L. E., Nguyen, R., Pfeiffer, H. P., & Scheel, M. A. 2018a, 1806.02349
  • Foucart et al. (2016a) Foucart, F. et al. 2016a, Phys. Rev., D93, 044019, 1510.06398
  • Foucart et al. (2018b) Foucart, F., Hinderer, T., & Nissanke, S. 2018b, 1807.00011
  • Foucart et al. (2015) Foucart, F. et al. 2015, Phys. Rev., D91, 124021, 1502.04146
  • Foucart et al. (2016b) Foucart, F., O’Connor, E., Roberts, L., Kidder, L. E., Pfeiffer, H. P., & Scheel, M. A. 2016b, Phys. Rev., D94, 123016, 1607.07450
  • Freiburghaus et al. (1999) Freiburghaus, C., Rosswog, S., & Thielemann, F.-K. 1999, Astrophys. J., 525, L121
  • Frensel et al. (2017) Frensel, M., Wu, M.-R., Volpe, C., & Perego, A. 2017, Phys. Rev., D95, 023011, 1607.05938
  • Fujibayashi et al. (2018) Fujibayashi, S., Kiuchi, K., Nishimura, N., Sekiguchi, Y., & Shibata, M. 2018, Astrophys. J., 860, 64, 1711.02093
  • Fujibayashi et al. (2017) Fujibayashi, S., Sekiguchi, Y., Kiuchi, K., & Shibata, M. 2017, Astrophys. J., 846, 114, 1703.10191
  • Galeazzi et al. (2013) Galeazzi, F., Kastaun, W., Rezzolla, L., & Font, J. A. 2013, Phys. Rev., D88, 064009, 1306.4953
  • Ghirlanda et al. (2018) Ghirlanda, G., et al. 2018, 1808.00469
  • Giacomazzo & Perna (2013) Giacomazzo, B., & Perna, R. 2013, Astrophys. J., 771, L26, 1306.1608
  • Goriely et al. (2011) Goriely, S., Bauswein, A., & Janka, H. T. 2011, Astrophys. J., 738, L32, 1107.0899
  • Gottlieb et al. (2008) Gottlieb, S., Ketcheson, D. I., & Shu, C.-W. 2008, Journal of Scientific Computing, 38, 251
  • Gourgoulhon et al. (1999) Gourgoulhon, E., Grandclément, P., Marck, J.-A., Novak, J., & Taniguchi, K. 1999, LORENE, paris Observatory, Meudon section - LUTH laboratory
  • Haggard et al. (2017) Haggard, D., Nynka, M., Ruan, J. J., Kalogera, V., Bradley Cenko, S., Evans, P., & Kennea, J. A. 2017, Astrophys. J., 848, L25, 1710.05852
  • Hallinan et al. (2017) Hallinan, G., et al. 2017, Science, 1710.05435
  • Hanauske et al. (2017) Hanauske, M., Takami, K., Bovard, L., Rezzolla, L., Font, J. A., Galeazzi, F., & Stöcker, H. 2017, Phys. Rev., D96, 043004, 1611.07152
  • Hempel & Schaffner-Bielich (2010) Hempel, M., & Schaffner-Bielich, J. 2010, Nucl. Phys., A837, 210, 0911.4073
  • Hilditch et al. (2013) Hilditch, D., Bernuzzi, S., Thierfelder, M., Cao, Z., Tichy, W., & Bruegmann, B. 2013, Phys. Rev., D88, 084057, 1212.2901
  • Hinderer et al. (2018) Hinderer, T., et al. 2018, 1808.03836
  • Hotokezaka et al. (2018a) Hotokezaka, K., Beniamini, P., & Piran, T. 2018a, 1801.01141
  • Hotokezaka et al. (2013a) Hotokezaka, K., Kiuchi, K., Kyutoku, K., Muranushi, T., Sekiguchi, Y., Shibata, M., & Taniguchi, K. 2013a, Phys. Rev., D88, 044026, 1307.5888
  • Hotokezaka et al. (2013b) Hotokezaka, K., Kiuchi, K., Kyutoku, K., Okawa, H., Sekiguchi, Y.-i., Shibata, M., & Taniguchi, K. 2013b, Phys. Rev., D87, 024001, 1212.0905
  • Hotokezaka et al. (2018b) Hotokezaka, K., Kiuchi, K., Shibata, M., Nakar, E., & Piran, T. 2018b, 1803.00599
  • Hotokezaka et al. (2011) Hotokezaka, K., Kyutoku, K., Okawa, H., Shibata, M., & Kiuchi, K. 2011, Phys. Rev., D83, 124008, 1105.4370
  • Hotokezaka et al. (2018c) Hotokezaka, K., Nakar, E., Gottlieb, O., Nissanke, S., Masuda, K., Hallinan, G., Mooley, K. P., & Deller, A. 2018c, 1806.10596
  • Hotokezaka et al. (2016) Hotokezaka, K., Nissanke, S., Hallinan, G., Lazio, T. J. W., Nakar, E., & Piran, T. 2016, Astrophys. J., 831, 190, 1605.09395
  • Hotokezaka & Piran (2015) Hotokezaka, K., & Piran, T. 2015, Mon. Not. Roy. Astron. Soc., 450, 1430, 1501.01986
  • Ishii et al. (2018) Ishii, A., Shigeyama, T., & Tanaka, M. 2018, Astrophys. J., 861, 25, 1805.04909
  • Just et al. (2015) Just, O., Bauswein, A., Pulpillo, R. A., Goriely, S., & Janka, H. T. 2015, Mon. Not. Roy. Astron. Soc., 448, 541, 1406.2687
  • Kasen et al. (2013) Kasen, D., Badnell, N. R., & Barnes, J. 2013, Astrophys. J., 774, 25, 1303.5788
  • Kasen et al. (2017) Kasen, D., Metzger, B., Barnes, J., Quataert, E., & Ramirez-Ruiz, E. 2017, Nature, 1710.05463, [Nature551,80(2017)]
  • Kasliwal et al. (2017) Kasliwal, M. M., et al. 2017, Science, 1710.05436
  • Kastaun et al. (2016) Kastaun, W., Ciolfi, R., & Giacomazzo, B. 2016, Phys. Rev., D94, 044060, 1607.02186
  • Kastaun & Galeazzi (2015) Kastaun, W., & Galeazzi, F. 2015, Phys. Rev., D91, 064027, 1411.7975
  • Kawaguchi et al. (2016) Kawaguchi, K., Kyutoku, K., Shibata, M., & Tanaka, M. 2016, Astrophys. J., 825, 52, 1601.07711
  • Kawaguchi et al. (2018) Kawaguchi, K., Shibata, M., & Tanaka, M. 2018, 1806.04088
  • Kiuchi et al. (2017) Kiuchi, K., Kawaguchi, K., Kyutoku, K., Sekiguchi, Y., Shibata, M., & Taniguchi, K. 2017, Phys. Rev., D96, 084060, 1708.08926
  • Kiuchi et al. (2018) Kiuchi, K., Kyutoku, K., Sekiguchi, Y., & Shibata, M. 2018, Phys. Rev., D97, 124039, 1710.01311
  • Kiuchi et al. (2010) Kiuchi, K., Sekiguchi, Y., Shibata, M., & Taniguchi, K. 2010, Phys. Rev. Lett., 104, 141101, 1002.2689
  • Korobkin et al. (2012) Korobkin, O., Rosswog, S., Arcones, A., & Winteler, C. 2012, Mon. Not. Roy. Astron. Soc., 426, 1940, 1206.2379
  • Kreiss & Oliger (1973) Kreiss, H. O., & Oliger, J. 1973, Methods for the approximate solution of time dependent problems (Geneva: GARP publication series No. 10)
  • Kulkarni (2005) Kulkarni, S. R. 2005, astro-ph/0510256
  • Kurganov & Tadmor (2000) Kurganov, A., & Tadmor, E. 2000, Journal of Computational Physics, 160, 241
  • Kyutoku et al. (2014) Kyutoku, K., Ioka, K., & Shibata, M. 2014, Mon. Not. Roy. Astron. Soc., 437, L6, 1209.5747
  • Lattimer (2012) Lattimer, J. M. 2012, Ann. Rev. Nucl. Part. Sci., 62, 485, 1305.3510
  • Lattimer & Schramm (1974) Lattimer, J. M., & Schramm, D. N. 1974, Astrophys. J., 192, L145
  • Lattimer & Swesty (1991) Lattimer, J. M., & Swesty, F. D. 1991, Nucl. Phys., A535, 331
  • Lazzati et al. (2018) Lazzati, D., Perna, R., Morsony, B. J., López-Cámara, D., Cantiello, M., Ciolfi, R., Giacomazzo, B., & Workman, J. C. 2018, Phys. Rev. Lett., 120, 241103, 1712.03237
  • Lee et al. (2009) Lee, W. H., Ramirez-Ruiz, E., & López-Cámara, D. 2009, ApJ, 699, L93, 0904.3752
  • Lehner et al. (2016) Lehner, L., Liebling, S. L., Palenzuela, C., Caballero, O. L., O’Connor, E., Anderson, M., & Neilsen, D. 2016, Class. Quant. Grav., 33, 184002, 1603.00501
  • Li & Paczynski (1998) Li, L.-X., & Paczynski, B. 1998, Astrophys. J., 507, L59, astro-ph/9807272
  • Lippuner et al. (2017) Lippuner, J., Fernández, R., Roberts, L. F., Foucart, F., Kasen, D., Metzger, B. D., & Ott, C. D. 2017, Mon. Not. Roy. Astron. Soc., 472, 904, 1703.06216
  • Lippuner & Roberts (2015) Lippuner, J., & Roberts, L. F. 2015, Astrophys. J., 815, 82, 1508.03133
  • Löffler et al. (2012) Löffler, F., et al. 2012, Class. Quant. Grav., 29, 115001, 1111.3344
  • Lowrie & Morel (2001) Lowrie, R., & Morel, J. 2001, J. Quant. Spectrosc. Radiat. Transfer, 69, 475
  • Lyman et al. (2018) Lyman, J. D., et al. 2018, Nat. Astron., 2, 751, 1801.02669
  • Malik et al. (2018) Malik, T., Alam, N., Fortin, M., Providência, C., Agrawal, B. K., Jha, T. K., Kumar, B., & Patra, S. K. 2018, 1805.11963
  • Margalit & Metzger (2017) Margalit, B., & Metzger, B. D. 2017, Astrophys. J., 850, L19, 1710.05938
  • Margutti et al. (2017) Margutti, R., et al. 2017, Astrophys. J., 848, L20, 1710.05431
  • Margutti et al. (2018) ——. 2018, Astrophys. J., 856, L18, 1801.03531
  • Martin et al. (2015) Martin, D., Perego, A., Arcones, A., Thielemann, F.-K., Korobkin, O., & Rosswog, S. 2015, Astrophys. J., 813, 2, 1506.05048
  • Metzger et al. (2010) Metzger, B. D., Arcones, A., Quataert, E., & Martínez-Pinedo, G. 2010, Mon. Not. Roy. Astron. Soc., 402, 2771, 0908.0530
  • Metzger et al. (2015) Metzger, B. D., Bauswein, A., Goriely, S., & Kasen, D. 2015, Mon. Not. Roy. Astron. Soc., 446, 1115, 1409.0544
  • Metzger & Fernández (2014) Metzger, B. D., & Fernández, R. 2014, Mon. Not. Roy. Astron. Soc., 441, 3444, 1402.4803
  • Metzger et al. (2010) Metzger, B. D. et al. 2010, Mon. Not. Roy. Astron. Soc., 406, 2650, 1001.5029
  • Metzger et al. (2008) Metzger, B. D., Piro, A. L., & Quataert, E. 2008, Mon. Not. Roy. Astron. Soc., 390, 781, 0805.4415
  • Metzger et al. (2009) ——. 2009, Mon. Not. Roy. Astron. Soc., 396, 304, 0810.2535
  • Metzger et al. (2018) Metzger, B. D., Thompson, T. A., & Quataert, E. 2018, Astrophys. J., 856, 101, 1801.04286
  • Meyer (1989) Meyer, B. S. 1989, Astrophys. J., 343, 254
  • Mooley et al. (2018) Mooley, K. P. et al. 2018, 1806.09693
  • Mooley et al. (2017) Mooley, K. P., et al. 2017, Nature
  • Most et al. (2018a) Most, E. R., Papenfort, L. J., Dexheimer, V., Hanauske, M., Schramm, S., Stöcker, H., & Rezzolla, L. 2018a, 1807.03684
  • Most et al. (2018b) Most, E. R., Weih, L. R., Rezzolla, L., & Schaffner-Bielich, J. 2018b, 1803.00549
  • Nakar & Piran (2011) Nakar, E., & Piran, T. 2011, Nature, 478, 82, 1102.1020
  • Nakar & Piran (2018) ——. 2018, Mon. Not. Roy. Astron. Soc., 478, 407, 1801.09712
  • Nandi et al. (2018) Nandi, R., Char, P., & Pal, S. 2018, 1809.07108
  • Neilsen et al. (2014) Neilsen, D., Liebling, S. L., Anderson, M., Lehner, L., O’Connor, E., & Palenzuela, C. 2014, Phys. Rev., D89, 104029, 1403.3680
  • Nicholl et al. (2017) Nicholl, M., et al. 2017, Astrophys. J., 848, L18, 1710.05456
  • O’Connor & Ott (2010) O’Connor, E., & Ott, C. D. 2010, Class. Quant. Grav., 27, 114103, 0912.2393
  • Oechslin et al. (2006) Oechslin, R., Janka, H. T., & Marek, A. 2006, Astron. Astrophys., astro-ph/0611047, [Astron. Astrophys.467,395(2007)]
  • Ozel & Freire (2016) Ozel, F., & Freire, P. 2016, Ann. Rev. Astron. Astrophys., 54, 401, 1603.02698
  • Palenzuela et al. (2013) Palenzuela, C., Lehner, L., Ponce, M., Liebling, S. L., Anderson, M., Neilsen, D., & Motl, P. 2013, Phys. Rev. Lett., 111, 061105, 1301.7074
  • Palenzuela et al. (2015) Palenzuela, C., Liebling, S. L., Neilsen, D., Lehner, L., Caballero, O. L., O’Connor, E., & Anderson, M. 2015, Phys. Rev., D92, 044045, 1505.01607
  • Pardo et al. (2018) Pardo, K., Fishbach, M., Holz, D. E., & Spergel, D. N. 2018, JCAP, 1807, 048, 1801.08160
  • Paschalidis et al. (2015) Paschalidis, V., East, W. E., Pretorius, F., & Shapiro, S. L. 2015, Phys. Rev., D92, 121502, 1510.03432
  • Paschalidis et al. (2018) Paschalidis, V., Yagi, K., Alvarez-Castillo, D., Blaschke, D. B., & Sedrakian, A. 2018, Phys. Rev., D97, 084038, 1712.00451
  • Perego et al. (2016) Perego, A., Cabezón, R., & Käppeli, R. 2016, Astrophys. J. Suppl., 223, 22, 1511.08519
  • Perego et al. (2017a) Perego, A., Radice, D., & Bernuzzi, S. 2017a, Astrophys. J., 850, L37, 1711.03982
  • Perego et al. (2014) Perego, A., Rosswog, S., Cabezón, R. M., Korobkin, O., Käppeli, R., Arcones, A., & Liebendörfer, M. 2014, Mon. Not. Roy. Astron. Soc., 443, 3134, 1405.6730
  • Perego et al. (2017b) Perego, A., Yasin, H., & Arcones, A. 2017b, J. Phys., G44, 084007, 1701.02017
  • Piran et al. (2013) Piran, T., Nakar, E., & Rosswog, S. 2013, MNRAS, 430, 2121, 1204.6242
  • Plewa & Müller (1999) Plewa, T., & Müller, E. 1999, Astron. Astrophys., 342, 179, astro-ph/9807241
  • Pollney et al. (2011) Pollney, D., Reisswig, C., Schnetter, E., Dorband, N., & Diener, P. 2011, Phys. Rev., D83, 044045, 0910.3803
  • Qian & Woosley (1996) Qian, Y. Z., & Woosley, S. E. 1996, Astrophys. J., 471, 331, astro-ph/9611094
  • Radice (2017) Radice, D. 2017, Astrophys. J., 838, L2, 1703.02046
  • Radice et al. (2017) Radice, D., Bernuzzi, S., Del Pozzo, W., Roberts, L. F., & Ott, C. D. 2017, Astrophys. J., 842, L10, 1612.06429
  • Radice et al. (2016a) Radice, D., Bernuzzi, S., & Ott, C. D. 2016a, Phys. Rev., D94, 064011, 1603.05726
  • Radice et al. (2016b) Radice, D., Galeazzi, F., Lippuner, J., Roberts, L. F., Ott, C. D., & Rezzolla, L. 2016b, Mon. Not. Roy. Astron. Soc., 460, 3255, 1601.02426
  • Radice et al. (2018a) Radice, D., Perego, A., Bernuzzi, S., & Zhang, B. 2018a, 1803.10865
  • Radice et al. (2018b) Radice, D., Perego, A., Hotokezaka, K., Bernuzzi, S., Fromm, S., & Roberts, L. F. 2018b, In prep.
  • Radice et al. (2018c) Radice, D., Perego, A., Zappa, F., & Bernuzzi, S. 2018c, Astrophys. J., 852, L29, 1711.03647
  • Radice & Rezzolla (2012) Radice, D., & Rezzolla, L. 2012, Astron. Astrophys., 547, A26, 1206.6502
  • Radice et al. (2014a) Radice, D., Rezzolla, L., & Galeazzi, F. 2014a, Mon. Not. Roy. Astron. Soc., 437, L46, 1306.6052
  • Radice et al. (2014b) ——. 2014b, Class. Quant. Grav., 31, 075012, 1312.5004
  • Radice et al. (2015) ——. 2015, ASP Conf. Ser., 498, 121, 1502.00551
  • Reisswig et al. (2013a) Reisswig, C., Haas, R., Ott, C. D., Abdikamalov, E., Mösta, P., Pollney, D., & Schnetter, E. 2013a, Phys. Rev., D87, 064023, 1212.1191
  • Reisswig et al. (2013b) Reisswig, C., Ott, C. D., Abdikamalov, E., Haas, R., Moesta, P., & Schnetter, E. 2013b, Phys. Rev. Lett., 111, 151101, 1304.7787
  • Rezzolla et al. (2010) Rezzolla, L., Baiotti, L., Giacomazzo, B., Link, D., & Font, J. A. 2010, Class. Quant. Grav., 27, 114105, 1001.3074
  • Rezzolla et al. (2011) Rezzolla, L., Giacomazzo, B., Baiotti, L., Granot, J., Kouveliotou, C., & Aloy, M. A. 2011, Astrophys. J., 732, L6, 1101.4298
  • Rezzolla et al. (2018) Rezzolla, L., Most, E. R., & Weih, L. R. 2018, Astrophys. J., 852, L25, 1711.00314
  • Rosswog & Davies (2003) Rosswog, S., & Davies, M. B. 2003, Mon. Not. Roy. Astron. Soc., 345, 1077, astro-ph/0110180
  • Rosswog et al. (2017) Rosswog, S., Feindt, U., Korobkin, O., Wu, M. R., Sollerman, J., Goobar, A., & Martinez-Pinedo, G. 2017, Class. Quant. Grav., 34, 104001, 1611.09822
  • Rosswog & Liebendoerfer (2003) Rosswog, S., & Liebendoerfer, M. 2003, Mon. Not. Roy. Astron. Soc., 342, 673, astro-ph/0302301
  • Rosswog et al. (1999) Rosswog, S., Liebendoerfer, M., Thielemann, F. K., Davies, M. B., Benz, W., & Piran, T. 1999, Astron. Astrophys., 341, 499, astro-ph/9811367
  • Rosswog et al. (2013) Rosswog, S., Piran, T., & Nakar, E. 2013, Mon. Not. Roy. Astron. Soc., 430, 2585, 1204.6240
  • Rosswog et al. (2003) Rosswog, S., Ramirez-Ruiz, E., & Davies, M. B. 2003, Mon. Not. Roy. Astron. Soc., 345, 1077, astro-ph/0306418
  • Rosswog et al. (2018) Rosswog, S., Sollerman, J., Feindt, U., Goobar, A., Korobkin, O., Fremling, C., & Kasliwal, M. 2018, Astron. Astrophys., 615, A132, 1710.05445
  • Ruan et al. (2018) Ruan, J. J., Nynka, M., Haggard, D., Kalogera, V., & Evans, P. 2018, Astrophys. J., 853, L4, 1712.02809
  • Ruffert et al. (1996) Ruffert, M. H., Janka, H. T., & Schaefer, G. 1996, Astron. Astrophys., 311, 532, astro-ph/9509006
  • Ruiz et al. (2016) Ruiz, M., Lang, R. N., Paschalidis, V., & Shapiro, S. L. 2016, Astrophys. J., 824, L6, 1604.02455
  • Ruiz et al. (2018) Ruiz, M., Shapiro, S. L., & Tsokaros, A. 2018, Phys. Rev., D97, 021501, 1711.00473
  • Schnetter et al. (2004) Schnetter, E., Hawley, S. H., & Hawke, I. 2004, Class. Quant. Grav., 21, 1465, gr-qc/0310042
  • Sekiguchi (2010) Sekiguchi, Y. 2010, Prog. Theor. Phys., 124, 331, 1009.3320
  • Sekiguchi et al. (2011) Sekiguchi, Y., Kiuchi, K., Kyutoku, K., & Shibata, M. 2011, Phys. Rev. Lett., 107, 051102, 1105.2125
  • Sekiguchi et al. (2015) ——. 2015, Phys. Rev., D91, 064059, 1502.06660
  • Sekiguchi et al. (2016) Sekiguchi, Y., Kiuchi, K., Kyutoku, K., Shibata, M., & Taniguchi, K. 2016, Phys. Rev., D93, 124046, 1603.01918
  • Shakura & Sunyaev (1973) Shakura, N. I., & Sunyaev, R. A. 1973, Astron. and Astrophys., 24, 337
  • Shapiro & Teukolsky (1983) Shapiro, S. L., & Teukolsky, S. A. 1983, Black holes, white dwarfs, and neutron stars: The physics of compact objects
  • Shibata (2016) Shibata, M. 2016, Numerical Relativity (World Scientific Publishing Co)
  • Shibata et al. (2017a) Shibata, M., Fujibayashi, S., Hotokezaka, K., Kiuchi, K., Kyutoku, K., Sekiguchi, Y., & Tanaka, M. 2017a, Phys. Rev., D96, 123012, 1710.07579
  • Shibata et al. (2017b) Shibata, M., Kiuchi, K., & Sekiguchi, Y.-i. 2017b, Phys. Rev., D95, 083005, 1703.10303
  • Shibata & Taniguchi (2006) Shibata, M., & Taniguchi, K. 2006, Phys. Rev., D73, 064027, astro-ph/0603145
  • Shibata et al. (2003) Shibata, M., Taniguchi, K., & Uryu, K. 2003, Phys. Rev., D68, 084020, gr-qc/0310030
  • Shibata et al. (2005) ——. 2005, Phys. Rev., D71, 084021, gr-qc/0503119
  • Shibata & Uryu (2000) Shibata, M., & Uryu, K. 2000, Phys. Rev., D61, 064001, gr-qc/9911058
  • Siegel et al. (2014) Siegel, D. M., Ciolfi, R., & Rezzolla, L. 2014, Astrophys. J., 785, L6, 1401.4544
  • Siegel & Metzger (2017) Siegel, D. M., & Metzger, B. D. 2017, Phys. Rev. Lett., 119, 231102, 1705.05473
  • Siegel & Metzger (2018) ——. 2018, Astrophys. J., 858, 52, 1711.00868
  • Smartt et al. (2017) Smartt, S. J., et al. 2017, Nature, 1710.05841
  • Soares-Santos et al. (2017) Soares-Santos, M., et al. 2017, Astrophys. J., 848, L16, 1710.05459
  • Steiner et al. (2013) Steiner, A. W., Hempel, M., & Fischer, T. 2013, Astrophys. J., 774, 17, 1207.2184
  • Suresh (1997) Suresh, A. 1997, J. Comp. Phys., 136, 83
  • Symbalisty & Schramm (1982) Symbalisty, E., & Schramm, D. N. 1982, Astrophys. J. Lett., 22, 143
  • Tanaka & Hotokezaka (2013) Tanaka, M., & Hotokezaka, K. 2013, Astrophys. J., 775, 113, 1306.3742
  • Tanaka et al. (2017) Tanaka, M., et al. 2017, Publ. Astron. Soc. Jap., 1710.05850
  • Tanvir et al. (2017) Tanvir, N. R., et al. 2017, Astrophys. J., 848, L27, 1710.05455
  • Tews et al. (2018) Tews, I., Margueron, J., & Reddy, S. 2018, 1804.02783
  • Thielemann et al. (2017) Thielemann, F. K., Eichler, M., Panov, I. V., & Wehmeyer, B. 2017, Ann. Rev. Nucl. Part. Sci., 67, 253, 1710.02142
  • Troja et al. (2017) Troja, E., et al. 2017, Nature, 1710.05433
  • Tsang et al. (2018) Tsang, C. Y., Tsang, M. B., Danielewicz, P., Lynch, W. G., & Fattoyev, F. J. 2018, 1807.06571
  • Typel et al. (2010) Typel, S., Ropke, G., Klahn, T., Blaschke, D., & Wolter, H. H. 2010, Phys. Rev., C81, 015803, 0908.2344
  • van Eerten & MacFadyen (2011) van Eerten, H. J., & MacFadyen, A. I. 2011, ApJ, 733, L37, 1102.4571
  • van Riper & Lattimer (1981) van Riper, K. A., & Lattimer, J. M. 1981, Astrophys. J., 249, 270
  • Villar et al. (2017) Villar, V. A., et al. 2017, Astrophys. J., 851, L21, 1710.11576
  • Wanajo et al. (2014) Wanajo, S., Sekiguchi, Y., Nishimura, N., Kiuchi, K., Kyutoku, K., & Shibata, M. 2014, Astrophys. J., 789, L39, 1402.7317
  • Waxman et al. (2017) Waxman, E., Ofek, E., Kushnir, D., & Gal-Yam, A. 2017, 1711.09638
  • Wei et al. (2018) Wei, J. B., Figura, A., Burgio, G. F., Chen, H., & Schulze, H. J. 2018, 1809.04315
  • Wu et al. (2016) Wu, M.-R., Fernández, R., Martínez-Pinedo, G., & Metzger, B. D. 2016, Mon. Not. Roy. Astron. Soc., 463, 2323, 1607.05290
  • Wu et al. (2017) Wu, M.-R., Tamborra, I., Just, O., & Janka, H.-T. 2017, 1711.00477
  • Zhu et al. (2016) Zhu, Y.-L., Perego, A., & McLaughlin, G. C. 2016, Phys. Rev., D94, 105006, 1607.04671

Appendix A Finite-Resolution Effects

Here we discuss uncertainties in the ejecta mass, velocity, and composition due to finite-resolution effects. The mass of the dynamical ejecta has not yet reached convergence at the resolution adopted in our simulations. As can be seen in Table 2, the dynamical ejecta mass shows moderate variations and changes in a non monotonic way as we vary the grid resolution. Consequently, our estimate for the ejecta mass should only be considered as semi-quantitative. We estimate the relative uncertainty in the ejecta mass due to finite-resolution effects alone to be as large as 50%. We remark that, as discussed in Sec. 3.1, the differences between ours and other published results, as well as between results published by different groups, are even larger (see also Table 3).

Figure 32.— Dynamical ejecta sensitivity to resolution. Left panel: electron fraction. Right panel: asymptotic velocity. The total ejecta mass shows large fractional variation with resolution, however electron fraction, velocity distribution, and nucleosynthetic yields appear to be robust.

Despite the significant uncertainty in the total ejecta mass, we find the other properties of the outflow, velocity, composition, entropy, to be robust with resolution. As an example, we show in Fig. 32 histograms of the ejecta as a function of Y_{e} and of the asymptotic velocity for the SFHo_M135135_LK binary. The robustness of the ejecta properties suggests that, while the overall morphology of the outflow streams are well captured by the simulations, there are stochastic variations, for example resulting in the fallback of some of the material in these streams due to interaction with the merger debris, that affect the overall amount of matter entrained by these outflows (see also the discussion in Bauswein et al. 2013).

Appendix B Treatment of the Viscous Fluxes

The general-relativistic large-eddy simulation (GRLES) approach requires the inclusion of the term \partial_{j}(\sqrt{-g}\tau^{ij}) in the fluid momentum equations. This numerical derivative needs to be carefully treated to avoid the odd-even decoupling instability. To this aim we introduce the left and right biased finite differencing operators

\partial_{i}^{\pm}u=\pm\frac{u(x\pm h\mathbf{e}_{i})-u(x)}{h}\,, (B1)

where \mathchoice{e_{\kern 0.0pt\hbox{$\scriptstyle i\hbox{}$}}\kern-8.924864pt\kern 0% .0pt\kern 0.0pt{}^{\hbox{$\scriptstyle{}$}\kern 0.0pt\hbox{$\kern 0.0pt% \scriptstyle\hbox{}j$}}}{e_{\kern 0.0pt\hbox{$\scriptstyle i\hbox{}$}}\kern-8.% 924864pt\kern 0.0pt\kern 0.0pt{}^{\raise 0.43pt\hbox{$\scriptstyle{}$}\kern 0.% 0pt\raise 0.43pt\hbox{$\kern 0.0pt\scriptstyle\hbox{}j$}}}{e_{\kern 0.0pt\hbox% {$\scriptscriptstyle i\hbox{}$}}\kern-5.624914pt\kern 0.0pt\kern 0.0pt{}^{% \raise 0.43pt\hbox{$\scriptscriptstyle{}$}\kern 0.0pt\raise 0.43pt\hbox{$\kern 0% .0pt\scriptscriptstyle\hbox{}j$}}}{e_{\kern 0.0pt\hbox{$\scriptscriptstyle i% \hbox{}$}}\kern-5.624914pt\kern 0.0pt\kern 0.0pt{}^{\raise 0.43pt\hbox{$% \scriptscriptstyle{}$}\kern 0.0pt\raise 0.43pt\hbox{$\kern 0.0pt% \scriptscriptstyle\hbox{}j$}}}=\mathchoice{\delta_{\kern 0.0pt\hbox{$% \scriptstyle i\hbox{}$}}\kern-8.924864pt\kern 0.0pt\kern 0.0pt{}^{\hbox{$% \scriptstyle{}$}\kern 0.0pt\hbox{$\kern 0.0pt\scriptstyle\hbox{}j$}}}{\delta_{% \kern 0.0pt\hbox{$\scriptstyle i\hbox{}$}}\kern-8.924864pt\kern 0.0pt\kern 0.0% pt{}^{\raise 0.43pt\hbox{$\scriptstyle{}$}\kern 0.0pt\raise 0.43pt\hbox{$\kern 0% .0pt\scriptstyle\hbox{}j$}}}{\delta_{\kern 0.0pt\hbox{$\scriptscriptstyle i% \hbox{}$}}\kern-5.624914pt\kern 0.0pt\kern 0.0pt{}^{\raise 0.43pt\hbox{$% \scriptscriptstyle{}$}\kern 0.0pt\raise 0.43pt\hbox{$\kern 0.0pt% \scriptscriptstyle\hbox{}j$}}}{\delta_{\kern 0.0pt\hbox{$\scriptscriptstyle i% \hbox{}$}}\kern-5.624914pt\kern 0.0pt\kern 0.0pt{}^{\raise 0.43pt\hbox{$% \scriptscriptstyle{}$}\kern 0.0pt\raise 0.43pt\hbox{$\kern 0.0pt% \scriptscriptstyle\hbox{}j$}}} are the coordinate basis vectors and h is the grid spacing. For the covariant derivatives we use the finite-differencing operators

D_{i}^{\pm}u^{j}=\partial_{i}^{\pm}u^{j}+\mathchoice{\Gamma_{\kern 0.0pt\hbox{% $\scriptstyle\hbox{}ik$}}\kern-14.174784pt\kern 0.0pt\kern 0.0pt{}^{\hbox{$% \scriptstyle{}$}\kern 0.0pt\hbox{$\kern 0.0pt\scriptstyle j\hbox{}$}}}{\Gamma_% {\kern 0.0pt\hbox{$\scriptstyle\hbox{}ik$}}\kern-14.174784pt\kern 0.0pt\kern 0% .0pt{}^{\raise 0.43pt\hbox{$\scriptstyle{}$}\kern 0.0pt\raise 0.43pt\hbox{$% \kern 0.0pt\scriptstyle j\hbox{}$}}}{\Gamma_{\kern 0.0pt\hbox{$% \scriptscriptstyle\hbox{}ik$}}\kern-9.374857pt\kern 0.0pt\kern 0.0pt{}^{\raise 0% .43pt\hbox{$\scriptscriptstyle{}$}\kern 0.0pt\raise 0.43pt\hbox{$\kern 0.0pt% \scriptscriptstyle j\hbox{}$}}}{\Gamma_{\kern 0.0pt\hbox{$\scriptscriptstyle% \hbox{}ik$}}\kern-9.374857pt\kern 0.0pt\kern 0.0pt{}^{\raise 0.43pt\hbox{$% \scriptscriptstyle{}$}\kern 0.0pt\raise 0.43pt\hbox{$\kern 0.0pt% \scriptscriptstyle j\hbox{}$}}}u^{k}\,, (B2)

where \mathchoice{\Gamma_{\kern 0.0pt\hbox{$\scriptstyle\hbox{}ik$}}\kern-14.174784% pt\kern 0.0pt\kern 0.0pt{}^{\hbox{$\scriptstyle{}$}\kern 0.0pt\hbox{$\kern 0.0% pt\scriptstyle j\hbox{}$}}}{\Gamma_{\kern 0.0pt\hbox{$\scriptstyle\hbox{}ik$}}% \kern-14.174784pt\kern 0.0pt\kern 0.0pt{}^{\raise 0.43pt\hbox{$\scriptstyle{}$% }\kern 0.0pt\raise 0.43pt\hbox{$\kern 0.0pt\scriptstyle j\hbox{}$}}}{\Gamma_{% \kern 0.0pt\hbox{$\scriptscriptstyle\hbox{}ik$}}\kern-9.374857pt\kern 0.0pt% \kern 0.0pt{}^{\raise 0.43pt\hbox{$\scriptscriptstyle{}$}\kern 0.0pt\raise 0.4% 3pt\hbox{$\kern 0.0pt\scriptscriptstyle j\hbox{}$}}}{\Gamma_{\kern 0.0pt\hbox{% $\scriptscriptstyle\hbox{}ik$}}\kern-9.374857pt\kern 0.0pt\kern 0.0pt{}^{% \raise 0.43pt\hbox{$\scriptscriptstyle{}$}\kern 0.0pt\raise 0.43pt\hbox{$\kern 0% .0pt\scriptscriptstyle j\hbox{}$}}} are the Christoffel symbols of the Levi-Civita connection associated with the spatial metric \gamma_{ij}. These are computed using a standard centered 2nd order finite-differencing operator.

Using the biased finite-differencing operators we define

\tau_{ij}^{\pm}=-2\,\nu_{{}_{T}}\,(\rho+p)\,W^{2}\left[\frac{1}{2}\big{(}D_{i}% ^{\pm}v_{j}+D_{j}^{\pm}v_{i}\big{)}-\frac{1}{3}D_{k}^{\pm}v^{k}\gamma_{ij}% \right]\,. (B3)

Then we discretize the divergence of \tau^{ij} as

\partial_{j}(\sqrt{-g}\tau^{ij})\simeq\frac{1}{2}\left\{\partial_{j}^{-}\big{[% }\sqrt{-g}\,(\tau^{+})^{ij}\big{]}+\partial_{j}^{+}\big{[}\sqrt{-g}\,(\tau^{-}% )^{ij}\big{]}\right\}\,. (B4)

It is easy to verify that this yields a flux-conservative scheme. Exact conservation is also ensured at refinement level boundaries using the flux correction algorithm of Berger & Colella (1989) as implemented in Reisswig et al. (2013a).

Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description