Unequal Mass Binary Neutron Star Simulations with Neutrino Transport: Ejecta and Neutrino Emission

TL;DR

This study presents detailed simulations of unequal mass neutron star mergers, analyzing ejecta and neutrino emissions, revealing how binary parameters influence matter ejection, composition, and neutrino characteristics.

Contribution

The paper introduces twelve new simulations with advanced neutrino transport, exploring the effects of different equations of state and mass ratios on ejecta and neutrino emissions.

Findings

Dynamical ejecta mass exceeds 0.01 solar masses mainly for SFHo EOS.

Ejecta electron fraction distribution is broad, with mean around 0.2.

Average ejecta velocity ranges from 0.2c to 0.3c, decreasing with asymmetry.

Abstract

We present twelve new simulations of unequal mass neutron star mergers. The simulations were preformed with the SpEC code, and utilize nuclear-theory based equations of state and a two-moment gray neutrino transport scheme with an improved energy estimate based on evolving the number density. We model the neutron stars with the SFHo, LS220 and DD2 equations of state (EOS) and we study the neutrino and matter emission of all twelve models to search for robust trends between binary parameters and emission characteristics. We find that the total mass of the dynamical ejecta exceeds $0.01M_\odot$ only for SFHo with weak dependence on the mass-ratio across all models. We find that the ejecta have a broad electron fraction ($Y_e$) distribution ($\approx 0.06-0.48$), with mean $0.2$. $Y_e$ increases with neutrino irradiation over time, but decreases with increasing binary asymmetry. We also find that the models have ejecta with a broad asymptotic velocity distribution ($\approx 0.05-0.7c$). The average velocity lies in the range $0.2c - 0.3c$ and decreases with binary asymmetry. Furthermore, we find that disk mass increases with binary asymmetry and stiffness of the EOS. The $Y_e$ of the disk increases with softness of the EOS. The strongest neutrino emission occurs for the models with soft EOS. For (anti) electron neutrinos we find no significant dependence of the magnitude or angular distribution or neutrino luminosity with mass-ratio. The heavier neutrino species have a luminosity dependence on mass-ratio but an angular distribution which does not change with mass-ratio.

Tables

Table 1. Table 1: Parameters for the initial models presented in this paper. M 1 subscript 𝑀 1 M_{1} and M 2 subscript 𝑀 2 M_{2} are the gravitational masses of the two neutron stars, C 1 subscript 𝐶 1 C_{1} and C 2 subscript 𝐶 2 C_{2} are the compactness of the respective neutron stars and d 0 subscript 𝑑 0 d_{0} is the initial coordinate distance between the centers of stars. The last column gives the post-merger remnant as of 7.5-ms post-merger. If the binary collapsed before 7.5-ms a rough estimate for the time of collapse is shown in parenthesis. We have codenamed each simulation by joining the first letter of the EOS name with the digits of the gravitational masses for that system, e.g. the DD2 1.2 ​ M ⊙ + 1.44 ​ M ⊙ 1.2 subscript 𝑀 direct-product 1.44 subscript 𝑀 direct-product 1.2M_{\odot}+1.44M_{\odot} run is abbreviated D12144 throughout the paper.

Model	EOS	$M_1$	$M_2$	$q$	$C_1$	$C_2$	$d_0\left(km\right)$	Collapse?
D144144	DD2	1.44	1.44	1.0	.161	.161	48.7	No
D12132	DD2	1.2	1.32	.91	.134	.147	48.7	No
D12144	DD2	1.2	1.44	.83	.134	.161	48.7	No
D12156	DD2	1.2	1.56	.77	.134	.173	48.7	No
L144144	LS220	1.44	1.44	1.0	.175	.175	44.3	$\left(\sim 2ms\right)$
L12132	LS220	1.2	1.32	.91	.146	.161	44.3	No
L12144	LS220	1.2	1.44	.83	.146	.175	44.1	No
L12156	LS220	1.2	1.56	.77	.146	.191	44.3	$\left(\sim 4.5ms\right)$
S144144	SFHo	1.44	1.44	1.0	.179	.179	44.3	$\left(\sim .5ms\right)$
S12132	SFHo	1.2	1.32	.91	.148	.163	44.3	No
S12144	SFHo	1.2	1.44	.83	.148	.179	44.3	No
S12156	SFHo	1.2	1.56	.77	.148	.195	44.3	$\left(\sim 1.4ms\right)$

Table 2. Table 2: Finite difference grid sizes for the initial models. All simulations were set so that they would have a 300m resolution during merger. For reference, the DD2 1.44 ​ M ⊙ 1.44 subscript 𝑀 direct-product 1.44M_{\odot} star has a ∼ 10.48 ​ k ​ m similar-to absent 10.48 𝑘 𝑚 \sim 10.48km radius in our grid coordinates.

Name	$d{x}_{\text{ins}}$ (m)	$d{x}_{\text{mer}}$ (m)
D12132	279	300
D12144	280	300
D12156	279	300
D144144	273	300
S12132	251	300
S12144	250	300
S12156	238	300
S144144	238	300
L12132	252	300
L12144	253	300
L12156	252	300
L144144	245	300

Equations

$$d{s}^2$$

$$=-{\alpha }^{2}d{t}^2+{\gamma }_{ij}(d{x}^i+{\beta }^i)(d{x}^j+{\beta }^j)$$

$${\gamma }_{\alpha \beta }=g_{\alpha \beta }+n_{\alpha }{n}_{\beta },$$

$$R_{\alpha \beta }=8\pi \left(T_{\alpha \beta }-\frac{1}{2}{g}_{\alpha \beta }T\right),$$

$$g_{\alpha \beta }{\nabla }^{\zeta }{\nabla }_{\zeta }{x}^{\beta }=H_{\alpha }\left(x^{ϵ}\right),$$

$$\begin{array}{rl}{g}^{\delta \gamma }{\partial }_{\gamma }{\partial }_{\delta }{g}_{\alpha \beta }+{\partial }_{\beta }{g}^{\gamma \delta }{\partial }_{\gamma }{g}_{\alpha \delta }& +{\partial }_{\alpha }{g}^{\gamma \delta }{\partial }_{\gamma }{g}_{\beta \delta }+2{\partial }_{(\beta ,}{H}_{\alpha )}\\ -2H_{\delta }{\Gamma }_{\alpha \beta }^{\delta }+2{\Gamma }_{\delta \beta }^{\gamma }{\Gamma }_{\gamma \alpha }^{\delta }& =-8\pi (2T_{\alpha \beta }-g_{\alpha \beta }T),\end{array}$$

$${\Gamma }_{\alpha \beta }^{\gamma }=\frac{1}{2}{g}^{\gamma ϵ}\left({\partial }_{\beta }{g}_{\alpha ϵ}+{\partial }_{\alpha }{g}_{\beta ϵ}-{\partial }_{ϵ}{g}_{\alpha \beta }\right).$$

$$H_{\alpha }\left(x_c^i,t\right)=H_{\alpha }\left(x_c^i,0\right)\exp \left(-\frac{t^2}{{\tau }^2}\right),$$

$$H_{\alpha }={\left(\log \frac{\sqrt{\gamma }}{\alpha }\right)}^2\left(\frac{\sqrt{\gamma }}{\alpha }{t}_{\alpha }-{\gamma }_{\alpha i}\frac{{\beta }^i}{\alpha }\right),$$

$$H_{\alpha }(t)=H_{\alpha }\left(1-\exp \frac{-(t-t_{DH}{)}^2}{{\tau }_m^2}\right),$$

$${\partial }_t{g}_{\alpha \beta }$$

$${\partial }_t{\Pi }_{\alpha \beta }$$

$$=2\alpha g^{\zeta \delta }({\gamma }^{ij}{\Phi }_{i\zeta \alpha }{\Phi }_{j\delta \beta }-{\Pi }_{\zeta \alpha }{\Pi }_{\delta \beta }-g^{ϵ\sigma }{\Gamma }_{\alpha \zeta ϵ}{\Gamma }_{\beta \delta \sigma })$$

$$-2\alpha {\nabla }_{(\alpha }{H}_{\beta )}-\frac{1}{2}\alpha t^{\zeta }{t}^{\delta }{\Pi }_{\zeta \delta }{\Pi }_{\alpha \beta }-\alpha t^{\zeta }{\Pi }_{\zeta i}{\gamma }^{ij}{\Phi }_{j\alpha \beta }$$

$$+\alpha {\gamma }_0[2{\delta }_{(\alpha }^c{t}_{\beta )}-g_{\alpha \beta }{t}^{\zeta }](H_{\zeta }+{\Gamma }_{\zeta })-{\gamma }_1{\gamma }_2{\beta }^i{\Phi }_{i\alpha \beta }$$

$$-2\alpha (T_{\alpha \beta }-\frac{1}{2}{g}_{\alpha \beta }{T}^{\zeta \delta }{g}_{\zeta \delta }),$$

$${\partial }_t{\varphi }_{i\alpha \beta }$$

$$=\frac{1}{2}\alpha t^{\zeta }{t}^{\delta }{\Phi }_{i\zeta \delta }{\Pi }_{\alpha \beta }+\alpha {\gamma }^{jk}{t}^{\zeta }{\Phi }_{ij\zeta }{\Phi }_{k\alpha \beta }-\alpha {\gamma }_2{\Phi }_{i\alpha \beta },$$

$${\gamma }_0$$

$$+\frac{.2}{M_2}f(r_2^i,2.5R_2)+\frac{.075}{M_2}f(r_c^i,2.5d),$$

$${\gamma }_1$$

$${\gamma }_2$$

$$+\frac{3}{M_2}f(r_2^i,2.5R_2)+\frac{.075}{M_2}f(r_c^i,2.5d),$$

$$f(r^i$$

$$p^{\alpha }\left[\frac{\partial f_{(\nu )}}{\partial x^{\alpha }}-{\Gamma }_{\alpha \gamma }^{\beta }{p}^{\gamma }\frac{\partial f_{(\nu )}}{\partial p^{\beta }}\right]=C\left[f_{(\nu )}\right],$$

$$T_{rad}^{\mu \nu }=J{u}^{\mu }{u}^{\nu }+H^{\mu }{u}^{\nu }+H^{\nu }{u}^{\mu }+S^{\mu \nu },$$

$$p^{\alpha }=\nu (u^{\alpha }+l^{\alpha }),$$

$$J$$

$$H^{\mu }$$

$$S^{\mu \nu }$$

$$T_{rad}^{\mu \nu }=E{n}^{\mu }{n}^{\nu }+F^{\mu }{n}^{\nu }+F^{\nu }{n}^{\mu }+P^{\mu \nu },$$

$$h_{\alpha \beta }=g_{\alpha \beta }+u_{\alpha }{u}_{\beta }.$$

$$E$$

$$F_{\mu }$$

$$+W{v}_{\mu }{v}_{\nu }{H}^{\nu }+(g_{\mu \nu }-n_{\mu }{v}_{\nu })H^{\nu }+W{v}_{\mu }{v}_{\nu }{H}^{\nu }$$

$$+(g_{\mu \nu }+n_{\mu }{v}_{\nu })v_{\rho }{S}^{\nu \rho },$$

$$P_{\mu \nu }$$

$$+(g_{\mu \rho }-n_{\mu }{v}_{\rho })(g_{\nu \kappa }-n_{\nu }{v}_{\kappa })S^{\rho \kappa }$$

$$+W(g_{\rho \nu }-n_{\rho }{v}_{\nu })v_{\mu }{H}^{\rho },$$

$$u^{\mu }=W(n^{\mu }+v^{\mu }),$$

$${\partial }_t\tilde{E}$$

Full text

I Introduction

The binary NS merger GW170817 was a landmark event, combining the first detection of gravitational waves with an observation of a short gamma ray burst and a kilonova Abbott et al. (2017a, b, 2017); Ajello et al. (2018). Following the event, studies emerged analyzing many aspects of GW170817, from the internal structure of neutron stars The LIGO Scientific Collaboration et al. (2018a), to the production of short gamma-ray bursts Mochkovitch et al. (1993); Lee and Kluzniak (1999); Janka et al. (1999); Abbott et al. (2017b, 2017); Mooley et al. (2018) and the synthesis of r-process elements Li and Paczynski (1998); Lattimer and Schramm (1976); Rosswog et al. (1999); Kulkarni (2005); Metzger et al. (2010); Metzger (2017a); Kasliwal et al. (2017); Coulter et al. (2017); Chornock et al. (2017); Coulter et al. (2017); Cowperthwaite et al. (2017); Kasen et al. (2017); Kilpatrick et al. (2017); McCully et al. (2017); Nicholl et al. (2017); Pian et al. (2017); Smartt et al. (2017); Soares-Santos et al. (2017); Tanvir et al. (2017); Evans et al. (2017). With the increasing sensitivity of advanced gravitational wave interferometers many binary neutron star (BNS) detections are expected in the next decade Collaboration et al. (2018).

Numerical simulations of mergers play a crucial role in efforts to model the gravitational wave signal, predict the properties of its electromagnetic counterparts, and estimate the production of r-process elements from the merger. In this work, we focus on the matter and neutrino emissions from BNS mergers. In particular we look at the ejecta and neutrino emission from a new set of twelve BNS simulations which extends a previous set of four equal mass BNS simulations Foucart et al. (2016) to include mass-ratios not equal to one. While asymmetric mass-ratio BNS simulations have been studied before in the context of general-relativistic-radiation hydrodynamics, e.g. Sekiguchi et al. (2016a), Lehner et al. (2016), Radice et al. (2016), these previous studies use either a simpler neutrino scheme and/or different mass-ratios and EOS. Thus this paper adds to the ongoing effort to simulate BNS systems and characterize their observables.

Neutrino interactions were first included in general relativistic simulations of neutron star mergers through a simple leakage scheme Sekiguchi et al. (2011), based on approximate methods developed for Newtonian simulations Ruffert et al. (1996); Rosswog and Liebendörfer (2003). A leakage scheme uses the local properties of the fluid and an estimate of the neutrino optical depth to determine the amount of energy lost locally to neutrino-matter interactions, and the associated change in the composition of the fluid. Leakage schemes provide an order-of-magnitude accurate estimate of neutrino cooling in the post-merger remnant, and have thus been used to capture the first-order effect of neutrino-matter interactions in general relativistic simulations of compact binary mergers Sekiguchi et al. (2011); Wanajo et al. (2014); Lehner et al. (2016); Radice et al. (2016); Palenzuela et al. (2015); Deaton et al. (2013); Foucart et al. (2014). The inclusion of neutrino-matter effects with the leakage scheme, while crude, significantly affected the composition, morphology and total mass of the outflows, with some studies showing a factor 2 difference in total ejecta mass Radice et al. (2016). However, most implementations of leakage do not account for irradiation of low density regions by neutrinos emitted from hot, dense regions. This potentially leads to large errors in the composition of the outflows, mostly by underestimating the number of protons Foucart et al. (2016a, 2015). Accordingly, the simplest leakage schemes are very inaccurate when attempting to predict the properties of post-merger electromagnetic signals. The only general relativistic simulations going beyond leakage use a moment formalism with an analytic closure to approximate the Boltzmann equation Thorne (1981); Shibata et al. (2011). In particular, neutron star merger simulations have been performed with a gray two-moment scheme Foucart et al. (2016a, 2015); Sekiguchi et al. (2015, 2016a), in which the energy density and flux density of each neutrino species are evolved. In BNS mergers, the use of this moment formalism showed that a range of compositions and thus of nucleosynthesis outcomes, exists in the material ejected by the merger Wanajo et al. (2014).

Previous studies of asymmetric mass-ratio BNS systems with fully general-relativistic radiation-hydrodynamics include Sekiguchi et al. (2016a); Lehner et al. (2016); Radice et al. (2018). Sekiguchi et al. Sekiguchi et al. (2016a) studied two EOS, SFHo and DD2 with mass ratios between .86 and 1.0 and a fixed total mass of $2.7M_{\odot}$ to around 30-ms post-merger using a two-moment neutrino transport scheme. They found that for SFHo the ejecta mass depended weakly on mass-ratio, but the average electron number per baryon decreased with mass ratio. For DD2 these trends were reversed. They also found that only the soft EOS, SFHo, produced ejecta mass above $0.01M_{\odot}$. Lehner et al. Lehner et al. (2016) studied three EOS, NL3, SFHo and DD2 with mass ratios between .76 and 1.0 and a fixed total mass of $2.7M_{\odot}$ at 3-ms post-merger using a neutrino leakage scheme. They found that there was a greater ejecta mass with increasing binary asymmetry. Unlike Sekiguchi et al., Lehner et. al found that none of the EOS produced ejecta above $0.01M_{\odot}$. Finally, Radice et al. Radice et al. (2018) studied four EOS, BHB$\Lambda\phi$, SFHo, DD2 and LS220, with mass ratios between .85 and 1.0 and a fixed total mass of $2.7M_{\odot}$ to around 20-ms post-merger using a neutrino leakage scheme and a viscous hydrodynamics scheme. Radice et al./ found that none of their models produced ejecta over $0.01M_{\odot}$ and their numbers agreed with Lehner et al. Lehner et al. (2016). The discrepancy between these results could however be due different choices for the definition of the unbound material in these studies. This paper adds onto previous works in the following ways. First, we look at a different set of parameters not found in the above studies. We use the SFHo, LS220 and DD2 EOS to study the effects of EOS on the merger emissions. We use mass-ratios ranging from $q\sim.76-1$ to study the effects of mass asymmetry on emissions. Unlike the previous studies, we do not fix the total mass and allow it to vary from $\sim 2.5-2.9M_{\odot}$. On top of this, we use a new two-moment neutrino transport scheme which evolves the number density, allowing for consistent lepton number evolution Foucart et al. (2016b).

We organize the paper as follows. In Section II, we discuss the numerical implementation we use and the equations we solve, including the gray two moment scheme for neutrino transport. In the following sections of the paper, we discuss the matter and neutrino emission from a new set of twelve binary neutron star merger simulations, ranging in mass ratio and equation of state. Finally we conclude with ideas for future work. We use a system of units such that $c=G=M_{\odot}=1$, where c is the speed of light in vacuum, G is the gravitational constant, and $M_{\odot}$ is the mass of the Sun. We use Einstein’s convention of summation over repeated indices. Latin indices run over 1, 2, 3, while Greek indices run over 0, 1, 2, 3. The spacetime metric signature is $(-,+,+,+)$.

II Numerical Implementation

II.1 General Overview

We evolve Einstein’s equations and the general relativistic equations of ideal radiation-hydrodynamics using the Spectral Einstein Code (SpEC)spe . SpEC evolves those equations on two separate grids: a pseudospectral grid for Einstein’s equations, written in the generalized harmonic formulation Lindblom et al. (2006), and a finite volume grid for the general relativistic equations of neutrino-hydrodynamics, written in conservative form. The latter makes use of an approximate Riemann solver (HLL Harten et al. (1983)) and high-order shock capturing methods (fifth order WENO scheme Liu et al. (1994); Jiang and Shu (1996)), resulting in a second-order accurate evolution scheme. For the time evolution, we use a third-order Runge-Kutta algorithm. Finally, after each time step, the two grids communicate the required source terms, using a third-order accurate spatial interpolation scheme. Those source terms are the metric and its derivatives (from the pseudospectral grid to the finite volume grid) and the fluid variables, which are rest-mass density, pressure, spatial velocity, the Lorentz factor and enthalpy. The following sections will give more detail on individual segments of this numerical method.

In the following sections, we make use the 3+1 decomposition of the metric

[TABLE]

where $\alpha$ is the lapse, $\beta^{i}$ the shift, and $\gamma_{ij}$ the 3-metric on a slice of constant coordinate t. The extension of $\gamma_{ij}$ to the full 4-dimensional space is the projection operator:

[TABLE]

with $n_{\mu}=(-\alpha,0,0,0)$ the unit normal to a $t={\rm constant}$ slice.

II.2 Initial Data

Initial data for this simulation was produced by an BNS initial data solver based on the work of Foucart et al. for BHNS systems Foucart et al. (2008), which was built upon the elliptic solver Spells Pfeiffer et al. (2003) and further improved for BNS systems in Tacik et al. (2015) and Haas et al. (2016). We start by considering systems in quasi-equilibrium, where time derivatives vanish in a corotating frame. We take the metric to be conformally flat and solve for the lapse, shift, and conformal factor using the extended conformal thin sandwich (XCTS) equations Pfeiffer and York (2003). The matter in the stars is modeled as a cold perfect fluid with an irrotational velocity profile.

Low eccentricity can be achieved through an iterative procedure requiring the evolution of the system for 2-3 orbits in each iteration Pfeiffer et al. (2007). All of the simulations in this paper use that algorithm to achieve estimated eccentricities of approximately $e\sim 0.001$.

Neutrinos are initialized in thermal equilibrium with the fluid. As the initial data is composed of two cold neutron stars ($T=0.1{\rm MeV}$), the initial energy density of neutrinos is negligible compared to its value at any later time in the evolution.

II.3 Spacetime Evolution

The Einstein field equations can be written in the form

[TABLE]

where $R_{\alpha\beta}$ is the Ricci tensor, $g_{\alpha\beta}$ is the metric tensor and $T_{\alpha\beta}$ is the stress energy tensor with trace T. SpEC uses the generalized harmonic decomposition Lindblom et al. (2006) to write the Einstein field equations in a form that allows stable numerical computation.

In the generalized harmonic formalism, the evolution of the coordinates follows the wave equation

[TABLE]

where $g_{ab}$ is the spacetime metric, and $H_{\alpha}(x^{c})$ a set of four arbitrary functions. Using Eq. (5) we can rewrite Eq. (4) as Pretorius (2005)

[TABLE]

where the Christoffel symbols $\Gamma^{\delta}_{\alpha\beta}$ are defined by

[TABLE]

Eq. (6) introduces four independent gauge functions $H_{\alpha}$, which need to be chosen. At t = 0, we set

[TABLE]

with $\tau=\sqrt{d^{3}_{0}/M_{\infty}}$, $d_{0}$ the initial separation, $M_{\infty}$ is the total mass of the binary at infinite separation and $x^{i}_{c}$ are comoving spatial coordinates which follow the rotation and inspiral of the binary. The initial data is constructed in a gauge $H^{initial}_{a}$ that assumes the time derivatives in the comoving frame are zero. At the beginning of the simulation, we set $H_{\alpha}\left(x^{i}_{c},0\right)=\hat{H}_{\alpha}$, where $\hat{H}_{\alpha}$ is a tensor that agrees with $H^{initial}_{a}$ in a frame comoving with the grid and is constant in time. The gauge will thus evolve into the harmonic condition $H_{\alpha}=0$. During merger, we find that transitioning to the “damped harmonic” gauge condition leads to more accurate numerical evolution Szilágyi (2014). Defining

[TABLE]

we transition according to

[TABLE]

with $t_{DH}$ the time at which we turn on the damped harmonic gauge and $\tau_{m}=100M$.

With these gauge choices we solve a first order representation of the generalized harmonic system (Eq. 6) in which the fundamental variables are the spacetime metric $g_{ab}$, its spatial first derivatives $\phi_{iab}$ and its first derivatives in the direction normal to the $t=$ constant slice $\Pi_{ab}=t^{c}\partial_{c}g_{ab}$. The generalized harmonic first order system is then:

[TABLE]

This amounts to a symmetric hyperbolic system of 50 coupled nonlinear equations. The constraint damping parameters $\{\gamma_{0},\gamma_{1},\gamma_{2}\}$ are additional free parameters which dampen any constraint violating modes which may grow due to small numerical errors in the evolution. To set the $\gamma$ parameters we use a trial-and-error process which concludes when we find that the constraint-violating modes do not grow significantly over time. In practice, we set the parameters to

[TABLE]

where $r^{i}_{1,2,c}$ correspond to the coordinate locations of the first neutron star, second neutron star and center of mass respectively, $R_{1}/M_{1}$ is the radius and mass of the 1st neutron star, $R_{2}/M_{2}$ is the radius and mass of the second neutron star, $M=M_{1}+M_{2}$ and $d$ is the separation of the two neutron star centers.

II.4 Neutrino Evolution

We use a gray two moment scheme for neutrino transport introduced in Foucart et al. (2016a). By evolving the number density of neutrinos in addition to their energy and flux densities, this new scheme guarantees exact conservation of the total lepton number. It also provides is with a local estimate of the average neutrino energy, an important quantity given the strong dependence of neutrino-matter cross-sections on the energy of neutrinos. In this section we will give an overview of this algorithm.

For each species of (anti)neutrino $\nu_{i}$ we can describe the neutrinos by their distribution function $f_{\nu}(x^{\mu},p^{\mu})$ where $x^{\mu}=(t,x^{i})$ gives the time and the position of the neutrinos and $p^{\alpha}$ is the 4-momentum of the neutrinos. The distribution function $f_{(\nu)}$ evolves according to the Boltzmann equation:

[TABLE]

where the right-hand side includes all collisional processes (emissions, absorptions, scatterings). In general, this is a 7-dimensional problem which is extremely expensive to solve numerically. Approximations to the Boltzmann equation have thus been developed for numerical applications. In this work, we consider the moment formalism developed by Thorne Thorne (1981) and Shibata et al. Shibata et al. (2011), in which only the lowest moments of the distribution function in momentum space are evolved. We limit ourselves to the use of this formalism in the gray approximation, that is we only consider energy-integrated moments. We will consider three independent neutrino species: the electron neutrino $\nu_{e}$, the electron antineutrinos $\bar{\nu}_{e}$, and the heavy lepton neutrinos $\nu_{x}$. The latter represents the 4 species $(\nu_{\mu},\bar{\nu}_{\mu},\nu_{\tau},\bar{\nu}_{\tau})$. This merging is justified because the temperatures and neutrino energies reached in our merger calculations are low enough to suppress the formation of the corresponding heavy leptons. The presence of heavy leptons would then require including the charged current neutrino interactions that differentiate between these individual species.

In the gray approximation with the first two moments of the distribution function, we evolve for each species projections of the stress-energy tensor of the neutrino radiation $T^{\mu\nu}_{rad}$. We write

[TABLE]

with $H^{\mu}u_{\mu}=S^{\mu\nu}u_{\mu}=0$ and $u^{\mu}$ the 4-velocity of the fluid. We can decompose the momentum as follows

[TABLE]

with $l^{\alpha}u_{\alpha}=0$ and $l^{\alpha}l_{\alpha}=1$. With this decomposition, the energy $J$, flux $H^{\mu}$ and stress tensor $S^{\mu\nu}$ of the neutrino radiation as observed by an observer comoving with the fluid are related to the neutrino distribution function by

[TABLE]

where $\nu$ is the neutrino energy in the fluid frame, and $\int\mathrm{d}\Omega$ denotes integrals over solid angle on a unit sphere in momentum space. We also utilize the decomposition of $T^{\mu\nu}_{rad}$ in terms of the energy, flux and stress tensor observed by an inertial observer (i.e. an observer whose worldline is tangent to the normal vector $n^{\mu}$),

[TABLE]

with $F^{\mu}n_{\nu}=P^{\mu\nu}n_{\mu}=F^{t}=P^{t\nu}=0$.

We define a projection operator on the reference frame of an observer comoving with the fluid:

[TABLE]

We can then write equations relating the fluid frame moments to the inertial frame moments

[TABLE]

using the decomposition of the 4-velocity

[TABLE]

where $v^{\mu}n_{\mu}=0$ and $W=\sqrt{1+v_{i}v^{i}}$.

By taking moments of the Boltzmann equation, the evolution equations for $\tilde{E}=\sqrt{\gamma}E$ and $\tilde{F}^{i}=\sqrt{\gamma}F^{i}$ can then be written in conservative form

[TABLE]

where $\gamma$ is the determinant of $\gamma_{ij}$, $\tilde{P}_{ij}=\sqrt{\gamma}P_{ij}$, and $\tilde{S}_{rad}^{\alpha}=\sqrt{\gamma}S_{rad}^{\alpha}$ includes all collisional source terms. Additionally, we consider the number current density for each species of neutrino:

[TABLE]

where N is the number density of neutrinos, and $\mathcal{F}^{\mu}$ the number flux density. We can define the number current in the fluid frame by

[TABLE]

where $\langle\nu\rangle$ is the average neutrino energy and $\langle\nu^{F}\rangle$ is the flux-weighted average neutrino energy. Using Eq. 33 we can get an estimate for the energy

[TABLE]

where $W$ is the Lorentz factor introduced in Eq. 29. The evolution equation for $\tilde{N}:=\sqrt{\gamma}N$ is

[TABLE]

where the source term accounts for neutrino-matter interactions.

To close this system of equations we need three additional ingredients: a prescription for the computation of $P^{ij}(E,F_{i})$ which is called the closure relation, a prescription for the computation of the number flux $\mathcal{F}^{j}$ (specific to te evolution of the number density N in this paper) and the collisional source terms $\tilde{S}^{\alpha}_{\rm rad},C_{(0)}$.

For $P^{ij}(E,F_{i})$ we interpolate between optically thick and thin limits:

[TABLE]

Here the parameter p is known as the variable Eddington factor and our choice for the functional form of p in terms of the lower moments $H,J$ is known as the Minerbo closure Minerbo (1978). Our choices for $P^{ij}_{\rm thin},P^{ij}_{\rm thick}$ and $p$ are discussed in the Appendix of Foucart et al. (2015).

For $\mathcal{F}^{i}$, by definition we have

[TABLE]

where the flux-weighted average neutrino energy, $\langle\nu^{F}\rangle$ is computed in such a way to take the effects of a finite optical depth on the spectrum into account, see Foucart et al. (2016) for details.

For the source terms $\tilde{S}^{\alpha}$, we assume that the fluid has an energy-integrated emissitivity $\bar{\eta}$ due to the charged-current reactions

[TABLE]

as well as electron-positron pair annihilation

[TABLE]

plasmon decay

[TABLE]

and nucleon-nucleon Bremsstrahlung

[TABLE]

The inverse reactions are responsible for an energy-averaged absorption opacity $\bar{\kappa}_{a}$. We also consider an energy-averaged scattering opacity $\bar{\kappa}_{s}$ due to elastic scattering of neutrinos on nucleons and heavy nuclei. Neglecting other reactions (e.g. inelastic scatterings and $\nu\bar{\nu}$ annihilation) the source terms $S^{\alpha}$ can then be shown to be Shibata et al. (2011)

[TABLE]

The collisional source term for the number density $\bar{N}$ is given by

[TABLE]

Thus we need the energy integrated emissivities $(\bar{\eta},\bar{\eta}_{N})$ and the energy-averaged opacities $(\bar{\kappa}_{A},\bar{\kappa}_{S},\bar{\kappa}_{N})$ to compute the source terms. We begin by computing the energy-integrated energy and number emissivities $(\bar{\eta},\bar{\eta}_{N})$ and energy-averaged equilibrium opacities $(\bar{\kappa}^{eq}_{A},\bar{\kappa}^{eq}_{S})$, which assume a thermal distribution of neutrinos in equilibrium with the fluid. The emissivities are computed following Ruffert et al. Ruffert et al. (1996), except for nucleon-nucleon Bremsstrahlung for which the emissivity is computed following Burrows et al. Burrows et al. (2006). The equilibrium opacities are computed using Kirchoff’s Law:

[TABLE]

where $B_{(\nu)}$ is the blackbody spectrum in equilibrium with the fluid. This ensures that, in regions where neutrinos are in thermal equilibrium with the fluid, they have the correct energy density $\bar{\eta}/\bar{\kappa}^{eq}$. To account for the fact that the average energy of neutrinos may be very different from its equilibrium value, we then correct the opacities by making use of the fact that processes computed here have cross-sections scaling as $T_{\nu}^{2}$:

[TABLE]

where $T_{\nu}$ is the neutrino temperature computed from the neutrino energy and number density, assuming a blackbody spectrum (see Foucart et al. (2016a) for more details). Lastly we set

[TABLE]

where $\eta_{\nu}=\mu_{\nu}/T$, $\mu_{\nu}$ is the chemical potential of neutrinos in equilibrium with the fluid, and $F_{k}(\eta_{\nu})$ is the Fermi integral

[TABLE]

This choice ensures that the average energy of neutrinos in equilibrium with the fluid takes its expected value,

[TABLE]

II.5 Fluid Evolution

The neutron stars are described by ideal fluids with stress tensor

[TABLE]

where $\rho_{0}$ is the rest mass density, h the specific enthalpy, P the pressure and $u^{\mu}$ the 4-velocity. The general relativistic equations of hydrodynamics are evolved in conservative form, using the conservative variables

[TABLE]

The 3 + 1 stress-energy conservation equations for radiation-hydrodynamics $\nabla_{\nu}T^{\mu\nu}=-S^{\mu}_{\text{rad}}$ and the baryon conservation equations $\nabla_{\mu}(\rho_{0}u^{\mu})=0$ then become Shibata et al. (2011)

[TABLE]

where $v_{t}^{i}=u^{i}/u^{0}$ is the transport velocity.

We evolve these equations using standard high-order shock-capturing finite volume method. We compute values of the physical quantities $(\rho_{0},T,u_{i})$ on cell faces from their cell-average values using fifth order WENO5 reconstruction Borges et al. (2008), and compute numerical fluxes on faces using the approximate HLL Riemann solver Harten et al. (1983).

II.6 Composition Evolution

The fluid composition in our simulations is described by the electron fraction:

[TABLE]

where $n_{p}$ and $n_{n}$ are the proton and neutron number densities, respectively (the net electron number density $n_{e^{-}}-n_{e^{+}}=n_{p}$, due to charge neutrality in the fluid). From lepton number conservation, we have:

[TABLE]

where $\sum_{\nu}$ sums over neutrino species with $\textrm{sign}(\nu)$ set to 1 for $\nu_{e}$, -1 for $\bar{\nu}_{e}$, and 0 for $\nu_{x}$. Importantly for the electron fraction, evolving the neutrino number density $N$ frees us from having to guess at the average neutrino energy when computing the coupling to the fluid. It also guarantees that the source term for the evolution of the electron fraction of the fluid is fully consistent with the evolution of the neutrino number density, thus conserving the total lepton number of the system Foucart et al. (2016a).

II.7 Equation of State (EOS)

We use three finite-temperature, composition dependent nuclear-theory based equations of state. Two of them are based on relativistic mean field (RMF) models Walecka (1974) and one based on the single nucleus approximation for heavy nuclei Lattimer and Swesty (1991). They are:

DD2 Hempel et al. (2012): This EOS is based on nuclear statistical equilibrium with a finite volume correction coupled to a relativistic mean field theory for treating high-density nuclear matter. DD2 contains neutrons, protons, light nuclei such as deuterons, helions, tritons and alpha particles and heavy nuclei. DD2 does not satisfy the so-called flow-constraint Hempel et al. (2017). 2. 2.

LS220 Lattimer and Swesty (1991): is based on the single nucleus approximation for heavy nuclei where the thermal distribution of different nuclear species is replaced by a single representative heavy nucleus. LS220 contains neutrons, protons, alpha particles and heavy nuclei. LS220 does not satisfy the constraints from Chiral effective field theory Hempel et al. (2017). 3. 3.

SFHo Steiner et al. (2013): This EOS, like DD2, also uses an RMF model, containing neutrons, protons, light nuclei such as deuterons, helions, tritons and alpha particles and heavy nuclei. However, SFHo uses a different RMF parameterization, specifically designed to match neutron star properties as inferred by observations. SFHo shows some minor deviations from Chiral effective field theory calculations Hempel et al. (2017).

While each EOS shows deviations from theoretical calculations, they all have a radius, maximum mass and tidal deformalibity that are compatible with observations of isolated neutron stars Demorest et al. (2010); Hempel et al. (2017). The LS220 and SFHo EOSs are also compatible with all constraints derived from the gravitational wave signal emitted by GW170817 The LIGO Scientific Collaboration et al. (2018a, b), while the DD2 equation of state has a tidal deformability for $1.4M_{\odot}$ neutron stars just outside of the strictest bound derived using parametrized equation of state models The LIGO Scientific Collaboration et al. (2018b). For all three equations of state, we show mass-radius curves in Figure 1. These curves were computed by integrating the TOV equations for the EOS by using the table values at $T=0.1$ MeV in beta equilibrium. From Figure 1, we can clearly see that for each EOS there is a maximum mass for non-rotating isolated neutron stars. Furthermore we can associate with each EOS an average density at this maximum mass as $\langle\rho\rangle:=3M_{max}/4\pi R^{3}_{max}$ and the ratio of this average density and the central density $\rho_{c}$ of the star is an indication of how “stiff” the EOS is. Higher values of $\langle\rho\rangle/\rho_{c}$ correspond to stiffer EOS, while lower values correspond to softer EOS Bauswein et al. (2013). In our case, from stiffest to softest, the EOS are ranked DD2, LS220 and SFHo. Softer EOS tend to make smaller stars at a fixed mass. For a fiducial mass of $1.4M_{\odot}$, the radii are 13.2 km, 12.7 km, 11.9 km for DD2, LS220 and SFHo respectively. We can also introduce a useful quantity called the compactness, which is defined by $C=M/R$. Softer EOS tend to make more compact stars. The maximum mass however is not a good indicator of whether the post-merger remnant will promptly collapse to a black hole (BH) after merger because the remnant will most likely be differentially rotating and held up by centrifugal and thermal forces against collapse. A more accurate estimate for the maximum mass at which the merger promptly collapses to a BH is given by Bauswein et al. Bauswein et al. (2013) who ran simulations of BNS mergers with a smoothed particle hydrodynamics code that employs the conformal flatness approximation of the Einstein field equations and includes a GW backreaction scheme to determine a thershold mass for prompt collapse across 12 different realistic EOS. Bauswein et al. found the threshold masses for SFHo, LS220 and DD2 were $2.95M_{\odot}$, $3.05M_{\odot}$ and $3.35M_{\odot}$ respectively. Thus the SFHo and LS220 models are more likely to undergo prompt collapse for the higher end of the total masses considered in this paper ($\sim 2.9M_{\odot}$, see Sec. II.8).

II.8 Initial Models

We extend our previous work Foucart et al. (2016) and study the merger of unequal mass neutron star binaries with the neutrino transport scheme introduced in Foucart et al. (2016a) and discussed above. We study mass ratios between .76 and 1 and masses between 1.2$M_{\odot}$ and 1.56$M_{\odot}$ which are both within the ranges of current binary neutron star observations (see e.g. Lattimer (2012)). We focus here on the late stages of the coalescence, comprising the last $3-5$ orbits (depending on EoS) before merger and we stop evolving any systems once they collapse. The parameters for the initial models and the type of post-merger remnants they create (before or as of 7.5-ms post-merger) are shown in Table 1. The specific grid setup for each model will be discussed in the following section.

II.9 Grid Setup

Before the two neutron stars enter into contact, the pseudospectral grid on which we evolve Einstein’s equations takes advantage of the approximate spherical symmetry of the neighborhood of each star, and in the far-field region. The evolved spatial slice is decomposed into two small balls around the center of each neutron star, sets of spherical shells around each star, spherical shells in the wave-zone region far from the stars and distorted cubes to connect the three spherically symmetric regions. The inner ball is expanded into Zernike polynomials, the shells into Chebyshev polynomials (in radius) and spherical harmonics (in angle) and the distorted cubes in Chebyshev polynomials. The grid follows the centers of the neutron stars defined as the center of mass of the matter in the $x<0$ and $x>0$ half planes, through a simple rotation and scaling of the grid coordinates.

We maintain this grid decomposition for the evolution of Einstein’s equations up to the point at which the maximum density on the grid increases beyond the low-level oscillations observed during the inspiral. This rise in the density signifies the transition from two well-separated neutron star cores to a single, more massive object. At that point we switch to a grid which is fully centered on the coordinate center of mass of the system. This grid is made of a ball at the origin of the coordinate system, surrounded by 59 spherical shells extending to the outer edge of the computational domain. Both before and after merger, the outer boundary is located at $40d_{0}$, with $d_{0}$ the initial separation of the binary, provided in Table 1.

The finite volume grid on which we evolve the general relativistic equations of hydrodynamics is very simple. Before the two neutron stars get into contact, it is composed of two cubes, each centered on a neutron star and composed of $96^{3}$ cells. In the coordinate system comoving with the neutron star centers, the neutron stars expand as the binary inspirals. To avoid losing matter to the outer boundary of the finite volume grid, we expand the grid by $4.5\%$ every time the flux of matter across the outer boundary exceeds $0.015M_{\odot}s^{-1}$. As the inspiral lasts less than 10 ms, this implies a mass loss well below $10^{-4}M_{\odot}$ before merger. As the two neutron stars approach each other, the two finite volume boxes will eventually intersect. During merger, we would like to follow the forming massive neutron star remnant, the tidal tails, the accretion disk, and any ejected material. We switch to a finite difference grid centered on the forming remnant with 3 levels of refinement. Each level has, at our standard resolution $200^{2}\times 100$ cells, with the finest grid spacing listed in Table 2 and each coarser level increasing the grid spacing by a factor of 2. The lower number of cells in the vertical direction reflects the fact that the remnant is less extended in that direction, and thus that we do not need the finest grid to extend as far vertically as horizontally.

II.10 Ejecta Analysis

In a BNS merger, matter expelled at high velocity may ultimately become unbound from the central gravitational potential. Two indicators have been used to label matter unbound:

$\mathbf{u_{t}<-1}$: For a stationary spacetime, $u_{t}$ (the projection of the 4-velocity along the timelike Killing vector field) is a constant of motion for geodesics. Assuming the space is also asymptotically flat, the Lorentz factor $W$ satisfies $W=-u_{t}$ at infinity. Therefore we may flag matter as unbound using the condition $u_{t}<-1$. This assumes that the outflow is made of isolated particles following geodesics, thus neglecting the impact of pressure gradients and r-process heating on the fluid. It also assumes that the metric is time-independent. Neither assumption is entirely correct in the dynamical spacetime of a post-merger remnant. 2. 2.

$\mathbf{hu_{t}<-1}$: For a stationary relativistic fluid flow, the relativistic Bernoulli equation Rezzolla and Zanotti (2013) implies $hu_{t}$ is constant along fluid worldlines. In an asymptotically flat spacetime, we would expect $W=-u_{t}$ (if the fluid particles follow geodesics). Since the relativistic enthalpy $h$ is only defined up to a constant factor which can be set such that $h\leftarrow 1$ at spatial infinity. Therefore, we may flag matter as unbound using the condition $hu_{t}<-1$. The main difference with the previous criteria is that all of the thermal energy of the fluid is now assumed to be transformed into kinetic energy as the fluid decompresses. Energy deposition due to r-process heating is also treated differently. We implicitly assume that the difference in binding energy between particles in nuclear statistical equilibrium (NSE) at the current density, temperature, and composition, and their binding energy at low density and temperature but for the same $Y_{e}$, is entirely deposited/removed from the fluid’s kinetic energy. This neglects significant out-of-NSE evolution during r-process nucleosynthesis (see Foucart et al. (2017)).

We have found previously Foucart et al. (2016) that the second indicator $hu_{t}<-1$, the Bernoulli criterion, produces qualitatively more accurate results in SpEC, therefore all material labelled unbound in this paper uses this indicator. It is important to note that no indicator is exact. In fact, the Bernoulli criterion has been shown to result in as much as twice the ejected matter as the $u_{t}$ condition Kastaun and Galeazzi (2015).

With the Bernoulli criterion, we flag matter that is still on the grid if it is at least $50M_{\odot}$ away from the center of the remnant and compute the ejected mass as follows

[TABLE]

where $\mathcal{H}(\cdot)$ is the Heaviside function. While the $50M_{\odot}$ threshold is arbitrary, we have found that it comfortably excludes any ejecta very close to the remnant which may not ultimately become unbound.

We also measure unbound material that leaves the computational grid (for each run, the grid is roughly $200M_{\odot}\times 200M_{\odot}\times 100M_{\odot}$) using

[TABLE]

Here $v^{i}$ is the fluid 3-velocity, $W$ is the Lorentz factor, defined by $W:=(1-v^{i}v_{i})^{-1/2}$, $S$ is the grid boundary and $\rho_{0}\mathcal{H}(-hu_{t}-1)W(\alpha v^{i}-\beta^{i})n_{i}$ represents the flux of unbound material leaving the boundary. At a time $t$, we estimate the full ejecta as

[TABLE]

Finally, to compute average quantities for the ejected matter, we use the following definition of the mass-weighted average for a quantity X.

[TABLE]

II.11 Errors

Since we performed simulations only at one resolution, it is difficult to derive error estimates. However, Hotokezaka et al. Hotokezaka et al. (2013) find a $\sim 10\%$ relative error for their mass ejecta properties for unequal mass binaries and we have independently confirmed their results with our code using the same piecewise polytropes EOS they use at similar resolution. However in this work, we use tabulated EOS and have a slightly coarser resolution. Furthermore, in equal mass runs (see Foucart et al. (2016)) we have found that relative errors can be as high as $\sim 50\%$. Hotokezaka et al. Hotokezaka et al. (2013) find similar variations in the mass ejected by equal mass binaries, with some systems showing little sign of convergence even at higher resolution. We therefore conclude that our errors in the ejected mass may be as high as $\sim 50\%$ and conservatively set our error as

[TABLE]

with the lower bound $10^{-4}M_{\odot}$ coming from the fact that we ignore outflows of this size during the regridding in the inspiral stage. Practically speaking, this is likely to be a significant overestimate of the error for the average simulation in our dataset and is more representative of the error for the worst simulations presented here. Our numerical setup also ignores magnetic fields. Over the short time scales considered here, magnetohydrodynamics (MHD) effects are not expected to affect the evolution of neutron star remnant, but could drive additional outflows from the disk Kiuchi et al. (2014); Neilsen et al. (2014). Over longer time scales, magnetic fields would be critical to the spin evolution of the remnant neutron star, angular momentum transport, heating in the disk, and possibly the formation of relativistic jets and magnetically-driven outflows. General relativistic MHD simulations of postmerger disks show that up to $\sim 40\%$ of the accretion disk can be ejected over 9 seconds Fernandez et al. (2019). We discuss the difference magnetic fields and longterm neutrino-winds could make on ejecta estimates in Section III.2.3.

III Numerical Results

III.1 General Overview

Prior to merger, which is defined as the peak of the gravitational-wave amplitude, the compact objects inspiral around each other for 3-5 orbits, with the actual number varying between EOS. This stage is shown in first panel of Fig. 2 for the DD2 $1.44M_{\odot}+1.44M_{\odot}$ model and takes approximately eighteen milliseconds. After the stars start merging, which is shown in the second panel of Fig. 2, the properties of the system are then largely determined by the compactness and the mass-ratio of the pre-merger neutron stars. These parameters control two important mechanisms. The first mechanism is tidal forces, which rip off matter as the binary gets closer, ejecting matter at angles close to the orbital plane. The second mechanism is the contraction and recoil of the cores (called core-bounce), which produces shock waves that eject matter quasi-spherically. From these two mechanisms, we can make predictions on expected trends. For example, as the mass-ratio increases, the less massive star gets more tidally disrupted and thus on merger, the effect of the shock on ejection of matter is lessened, but the effect of tidal torque propelling ejecta off the grid is increased. Similarly, as the EOS softens, the compactness of the star increases and the effect of core-bounces on the ejection of matter increases. In general we might expect that changing the mass-ratio might have a small effect on the emission of SFHo and LS220, but a larger effect on the ejecta of DD2 (see e.g Sekiguchi et al. (2016a)). As we see in the third panel of Fig 2, the DD2 EOS showcases well defined tidal tails larger than eighty kilometers in length around three milliseconds post-merger. This is to be expected from a very stiff EOS, which creates large stars that are torn apart by tidal forces in the late inspiral. A few milliseconds after this, at around 7.5 milliseconds post-merger (as shown in the fourth panel of Fig. 2), the tidal tails have left the grid and we are left with a large, circularized, 125 km radii disk surrounding the remnant with matter outflows continuously pouring off and leaving the grid.

Differences in the EOS can have a large impact on the evolution and emission of the system. With the SFHo equation of state, i.e. for the most compact neutron stars, a compact core forms rapidly. In the higher total mass models, the SFHo star collapses promptly to a BH, within a few ms. For the other equations of state (LS220, DD2), more strongly developed tidal features appear at merger. Figure 3 shows each of the $1.2M_{\odot}+1.44M_{\odot}$ models across the three EOS at 3 ms post-merger.

Finally, at around 7.5 ms after merger, the tidal tail has left the computational grid and the remnant is fully formed with an accretion torus. Figure 4 shows a snapshot of the density profile, temperature profile and electron fraction profile of the binaries at 7.5 ms for the mass-ratio $(1:1.2)$ initial models. As we can see the SFHo EOS has a more dense and hot core than the other two EOS. This is expected as the core bounce for SFHo is much more violent and it occurs deeper in the gravitational potential of the system due to the softness of the EOS. The outflow from the SFHo remnant is also more symmetric because outflows powered by core bounce tend to have a more quasi-spherical morphology than outflows from tidal tails. Stiff EOSs like DD2 undergo less shock heating, but tend to have larger tidal tails. As cold tidal tails are typically neutron-rich, the ejecta in the $z=0$ plane also has a lower $Y_{e}$. The electron fraction is also noticeably higher in the disk region for SFHo. This could be due to the fact the SFHo remnant is much hotter and therefore neutrino irradiation effects cause a higher $Y_{e}$.

III.2 Matter outflow

III.2.1 General properties

We have compiled a table of the general properties of the ejecta up to either collapse or 7.5 ms post-merger in Table 3. To help analyze the ejecta properties, we distinguish between two different regions in our grid: polar and equatorial. The polar ejecta is defined to be any matter within an angle $\theta<30^{\circ}$ of the z-axis (the direction of the total angular momentum of the system), while the equatorial ejecta is the rest of the ejected matter. From Table 3 we first notice that for each EOS, the majority of the ejecta is equatorial. This is not surprising because tidal disruption tends to produce outflows that remain close to the orbital plane, while the main contribution to the dynamical polar ejecta is from quasi-spherical matter ejection during core-bounce and subsequent oscillations of the post-merger remnant. Given that for spherical mass ejection only $\sim 13\%$ of the matter has $\theta<30^{\circ}$, even that quasi-spherical ejecta contributes more to the mass of equatorial ejecta than to the polar ejecta. A more accurate analysis of this table is thus that simulations with the DD2 EOS and mass ratio $(1:1,1:1.1)$, as well as the $(1:1.1)$ simulation with the LS220 EOS have $(10-12)\%$ of their ejecta at $\theta<30^{\circ}$ and are thus nearly compatible with spherically symmetric mass ejection. The other configurations have a clear excess of equatorial ejecta, particularly the more asymmetric binaries with the stiffest (DD2) EOS, which have the strongest tidal features.

Secondly, Table 3 shows that the LS220 and SFHo $1.44M_{\odot}+1.44M_{\odot}$ models, which collapse promptly after the merger, have little to no ejecta. Third, the softest EOS, SFHo, has the most ejecta if the remnant avoids prompt collapse ($\sim 0.02M_{\odot}$), while the LS220 models eject, on average, less mass than the DD2 models. Since LS220 is considered to be softer than DD2, it may be a bit surprising that it ejects less matter. However, as the LS220 EOS has a lower maximum mass than the DD2 EOS, two of the LS220 models lead to relatively rapid collapse of the post-merger remnant to a black hole, which is known to suppress mass ejection. Not all results fit a clear and understandable narrative, however: none of these arguments explain the low mass ejection of the non-collapsing model L12132, or the extremely high mass ejection of the collapsing model S12156. We also note that the dependence of the ejected mass in the mass ratio is clearly non-monotonous. This was also observed by Sekiguchi et al. Sekiguchi et al. (2016a) for their SFHo models (but not their DD2 models). Our results thus add to the growing evidence that while overall trends in mass ejection can be associated with the compactness of neutron stars, the mass asymmetry of the binary, the total mass of the system, and the maximum mass of isolated neutron stars, there is also a significant scatter around ’best fit’ models that cannot be easily understood at this point (see also Dietrich and Ujevic (2017); Radice et al. (2018)).

Turning to the electron fraction and velocity, we see that $\left<Y_{e}\right>$ and $\left<v_{\infty}\right>$ are much higher in the polar regions and decrease roughly with mass ratio. $\left<Y_{e}\right>$ tends to hover around $\approx 0.2$ while $\left<v_{\infty}\right>$ varies between $0.2c$ and $0.3c$. These trends will be touched on in the coming section, where we look at the $Y_{e}$ and $v_{\infty}$ distributions.

Table 3 quantifies the unbound matter both on the grid (greater than $50M_{\odot}$ from the remnant) and off the grid. The material on the grid that is still labelled unbound at 7.5 ms is most likely due to neutrino-driven wind. To illustrate that matter ejection continues late in our simulations, we show in Figure 5 the $1.2M_{\odot}+1.32M_{\odot}$ models for all EOSs. The black contours separate the bound and unbound matter, while the white lines are isocontours of the density. Lastly, the colormap denotes the electron fraction, $Y_{e}$. We note that there is still a small amount of matter labelled unbound, particularly in the SFHo models. Also we again note that the accretion disk is more proton-rich for the SFHo EOS than for the other two systems. Both effects are natural consequences of the hotter post-merger remnant for the SFHo equations of state, which leads to higher neutrino luminosities, stronger neutrino-driven winds, and higher neutrino irradiation of the disk.

III.2.2 $Y_{e}$ and $v_{\infty}$ distributions

The electron fraction and ejecta velocity are important in determining the attributes of the kilonova and the outcome of r-process nucleosynthesis in the outflows. Indeed, the ejecta velocity has a significant impact on the evolution timescale and brightness of the electromagnetic emission Barnes and Kasen (2013), while the electron fraction is maybe the most important parameter in determining which elements are produced by the r-process. Roughly speaking, for $Y_{e}\lesssim 0.25$ nucleosynthesis results in the production of heavy r-process elements (mass number $A\gtrsim 120$), while less neutron rich ejecta produces lower mass elements Lippuner and Roberts (2015). This has consequences beyond astrophysical nucleosynthesis: in neutron-rich environments, r-process nucleosynthesis produces high-opacity lanthanides and actinides. A higher opacity means that the thermal emission from the ejecta only becomes observable at later times, when the ejecta is cooler. Accordingly, a neutron-rich ejecta typically produces an infrared, week-long kilonova, while a neutron-poor ejecta leads to an optical, day long kilonova Kasen et al. (2013); Barnes and Kasen (2013).

Figure 6 shows the inferred velocity of the ejecta at infinity across EOS for both the equatorial and polar regions at 7.5 ms post-merger. First we notice that the polar ejecta generally have a larger velocity than the equatorial ejecta. This is most likely due to the polar ejecta originating primarily from core-bounce and the collision of the two neutron stars, which produce faster moving ejecta than the more equatorial-plane bounded tidal ejecta. Furthermore, the faster-moving core-bounce ejecta traveling in the equatorial-plane may collide with tidal ejecta and be slowed down. We also notice that overall, most of the ejecta is at the lower end of the velocity range, $v_{\infty}\sim 0.2c$ with very little ejecta at $v_{\infty}>0.8c$ for any of the models pictured here. LS220 appears to eject matter with smaller velocities, although not very conclusively considering the small amount of matter ejected by that simulation. Thus, across EOSs, we cannot establish any robust differences in velocity distribution.

Figure 7 shows the electron fraction of ejected material for different EOSs in both the equatorial and polar regions. Qualitatively, the shapes of the distributions appear mostly EOS-independent. The equatorial ejecta have a broad range of electron fraction values whereas the polar ejecta, which make a much smaller fraction of the ejected mass, take on much higher $Y_{e}$ values on average. The production of dynamical ejecta with a broad $Y_{e}$ distribution and the presence of higher-$Y_{e}$ outflows in the polar regions are both robust results observed in all simulations of binary neutron star mergers using advanced neutrino transport methods Wanajo et al. (2014); Sekiguchi et al. (2016a); Foucart et al. (2016, 2016a).

Figure  6 and  7 investigated EOS-dependence at fixed component-masses. We now consider varying the component masses, while keeping the EOS fixed. Figure 8 shows the distribution of $v_{\infty}$ for the SFHo EOS at different mass ratios. As the system becomes more asymmetric, the average velocity of the outflows generally decreases, and the tail of the distribution at high velocities is suppressed. The lack of high-velocity ejecta for more asymmetric binaries could affect the ability of these systems to power fast, blue/UV emission Metzger (2017a). The decrease in the average velocity of the ejecta as the binary becomes more asymmetric is notable, as that trends is exactly the opposite of what is observed in black hole-neutron star mergers Kawaguchi et al. (2016); Foucart et al. (2016b). We interpret this difference as being due to the impact of the core-bounce ejecta, which is absent in black hole-neutron star mergers: while the velocity of the tidal ejecta may increase with the asymmetry of the system, the relative importance of the (faster) core-bounce ejecta decreases for more asymmetric systems, thus reducing the average velocity of the ejecta for binary neutron star systems.

Finally, Figure 9 shows the distribution of electron fraction for the SFHo EOS at different mass ratios. We see that as mass asymmetry increases, the average $Y_{e}$ of the ejecta decreases. The distributions remain broad, with higher values in the polar regions.

III.2.3 Long-term Emission

Analysis of the kilonova emission associated with GW170817 showed that the ejecta mass was on the order of $\sim 0.05M_{\odot}$ Metzger (2017b); Shibata et al. (2017). This is above the estimates from dynamical ejecta seen in Table 3, or for that matter above the mass of dynamical ejecta in any low-eccentricity, low-spin neutron star merger simulation performed to date (see e.g. Dietrich and Ujevic (2017) for a summary). This is not really a problem, however: our analysis ignores other forms of ejecta that become more important over longer evolution time scales. After merger, magnetically-driven Siegel and Metzger (2017); Fernandez et al. (2019) and neutrino-driven Just et al. (2015) winds as well as outflows due to viscous angular momentum transport and $\alpha$-particle recombination Fernández and Metzger (2013) should combine to unbind a large fraction of the remnant disk ($\sim 40\%$ for disks around black holes in Fernandez et al. (2019)). Based on these results for the long-term evolution of post-merger disk, we can provide order-of-magnitude estimates of the evolution of the remnant produced at the end of our simulations. We consider the following definition of the disk mass:

[TABLE]

which has been used in previous works, see Radice et al. (2018); Shibata et al. (2017). This definition is used for two reasons. Firstly, based on visualizations of the remnant, such as Fig. 5 one can clearly see that the $10^{13}g/cm^{3}$ density contour encapsulates only the core of the remnant. Secondly, it was shown that for densities approximately below $10^{13}g/cm^{3}$, the matter is rotationally bound with roughly a $r^{-3/2}$ fall-offHanauske et al. (2017). Unfortunately $M_{disk}$ defined by Eq. (65) drifts to larger values over time and indeed it has not converged yet at 7.5 ms. Nonetheless, we note similar evolution of this quantity over time for all EOS and mass-ratios, so using estimates at 7.5 ms will give us a crude estimate on the lower bound of the disk mass. We compute the disk masses at 7.5-ms post-merger except for the models that collapse (e.g. S12156, S144144, L12156, L144144) where we compute the disk mass right before collapse. The masses are summarized in Fig. 10.

We first note that the stiffer EOS produce larger disks. Moreover, disk mass increases as the degree of binary asymmetry increases. Most importantly, all but one simulation (L144144) can potentially eject enough mass to explain the kilonova observed after GW170917, with required efficiency for the conversion of disk mass to disk outflows ranging from $\sim 25\%$ to $\sim 50\%$. This is well within the (very uncertain) range of efficiencies allowed by current post-merger evolutions.

III.3 Neutrino Emission

Finally we turn to the neutrino emission of our models. Neutrinos play an important role in the post-merger evolution, not only cooling the central massive regions, but also transporting energy into the surface layers of the disk, providing a long-term outflow. On top of this, neutrinos deposit net lepton number through weak interactions which can drastically change the electron fraction of outflows and thus affect the kilonova emission. Furthermore, they can deposit large amounts of thermal energy into the polar regions through $\nu\bar{\nu}$ annihilation, potentially contributing to the creation of baryon free regions above the remnant and helping the production of jet-like emissions.

As seen in the bottom portion of Fig. 11, which shows the electron anti-neutrino energy density for the D144144 model at 3 ms post-merger, the main emission regions are the hot, dense parts of the remnant; the central core, and the shocked tidal arms. In these dense parts of the remant the neutrinos are advected with the flow. The main escape route for the neutrinos trapped in the dense part of the remnant is via the polar regions, which show a greater neutrino energy density in the top portion of Fig. 11 than the equatorial parts. Due to the more compact and hotter remnants, we would expect that the SFHo models might show greater neutrino luminosity than the other EOS and this is confirmed in Figure 12 for the electron neutrino luminosity. Although we do not plot them, the other species show exactly the same type of behavior, with SFHo having a much higher luminosity in all three species evolved in the transport scheme.

As the mass ratio changes, we find only mild variations in the neutrino luminosity as shown in Fig. 13. The heavy neutrino luminosity has the strongest dependance on the mass ratio, showing a sizable increase as the binary becomes more asymmetric. This is in contrast to Sekiguchi et al. Sekiguchi et al. (2016b) which finds a decrease in neutrino luminosity as the binary asymmetry increases. However, Sekiguchi et al. fix the total mass at $2.7M_{\odot}$. Since we only fix the lower mass neutron star in the binary, a larger mass asymmetry implies a higher total mass, more compact remnants, higher temperatures and higher neutrino emission. In our simulations, the impact of the total mass of the system is cancelled out by the change in mass ratio, which, as Sekiguchi showed, tends to decrease the neutrino luminosity of the post-merger remnant.

Next we take a look at the angular distribution of the neutrino emission. Figure 14 shows the neutrino flux density as a function of its angle with respect to the equatorial plane at 7.5 ms after merger. From this figure, one can clearly see that most of the neutrinos are emitted in the polar directions. Once a disk forms, the neutrinos are mostly confined within a 40 degrees cone around the poles, with an amplitude peak at 30 - 40 degrees from the poles. As argues in Foucart et al. (2016a), this peak is probably due to neutrinos beamed from the shocked tidal arms, which become less optically thick as time passes. The confinement of the neutrinos to the polar directions stems from the fact that neutrinos escape through the low-density regions above and below the disk and are confined by the optically thick accretion disk. As the EOS changes, the distributions stay qualitatively the same. The binary mass ratio as an equally negligible impact on the angular distribution of neutrinos, and is thus not depicted here.

While neutrino irradiation of the disk corona can drive outflows, the total mass of these outflows is negligible, especially over the fairly short time scales evolved here. We do however see some effects of neutrino irradiation in the electron fraction of the outflows. Figure 15 shows the electron fraction distributions for all outflows escaping the grid before 7.5-ms and 10-ms post-merger. As we can see, the disk outflows observed at the end of the simulation have a particularly high electron fraction. The reason for this is that in the presence of strong electron neutrino luminosity, the neutrino capture processes are activated, increasing the overall electron fraction.

IV Conclusion and Future Work

We presented a first set of SpEC simulations of unequal mass neutron star mergers using 3 nuclear-theory based equations of state (SFHo, LS220, DD2), and a two-moment grey neutrino transport scheme with an improved energy estimate based on evolving the number density. This set of simulations varied the equation of state of the neutron stars and the total mass and mass-ratio of the binary in order to determine any robust trends between these parameters and emission characteristics. We found the following

• •

For the models that promptly collapsed, namely L144144 and S144144, there was little to no ejecta and the disk masses were much smaller. This is consistent with other simulations of promptly-collapsing equal mass neutron star-neutron star binaries.

• •

The softest models, those of the SFHo EOS, had the most ejecta and SFHo was the only EOS with $M_{ej}\sim 10^{-2}M_{\odot}$.

• •

In a majority of the simulations, we found that matter was preferentially ejected more than $30^{\circ}$ away from the poles. This result is however not universal: for near equal-mass systems and stiff EOSs, we find that the outflows are dominated by the more spherically symmetric matter ejected at core-bounce and during subsequent oscillations of the neutron stars.

• •

The $Y_{e}$ distributions were very broad $(\sim 0.06-0.48)$ and quite similar across EOS and mass-ratio. $\langle Y_{e}\rangle$ is around $\sim 0.2$, but increases over time as the ejected matter is irradiated by neutrinos. We see a decrease in the average $Y_{e}$ with increasing binary asymmetry, but there is some scatter, so the relationship might be more complicated.

• •

The asymptotic velocity distributions were also very broad $(\sim 0.05c-0.7c)$, with only a small amount of matter in the high-velocity tail of these distributions and average velocities within the range $(0.2-0.3)c$. We see a decrease in the average asymptotic velocity for more asymmetric mass ratios.

• •

The disk masses of the models at 7.5 ms post-merger appear to increase with mass-ratio and stiffness of the EOS. The softest EOS, SFHo, has a disk with a much higher $Y_{e}$ than the other EOS, presumably because softest EOSs are associated with more compact stars, hotter post-merger remnants, and thus stronger neutrino irradiation of both the disk and the matter outflows.

• •

Neutrino emission was accordingly largest for the SFHo EOS. We found no dependence of the neutrino luminosity on the mass-ratio, however, except in the case of heavy-lepton neutrinos. Most of the neutrino emission was in the polar region and we did not find any dependence of its morphology in either the mass-ratio or the EOS.

The results presented here are limited by two important assumptions in our simulations. Firstly, the absence of magnetic fields. Over the short time scales considered here, magnetohydrodynamics effects are not expected to affect the evolution of the neutron star remnant, but could drive additional outflows from the disk Kiuchi et al. (2014); Neilsen et al. (2014). Over longer time scales, magnetic fields would be critical to the spin evolution of the remnant neutron star, angular momentum transport, heating in the disk, and possibly the formation of relativistic jets and magnetically-driven outflows. The second limitation is the relatively small numerical grid on which we evolve the equations of general relativistic hydrodynamics. The grid limits our ability to measure the mass and properties of the outflows accurately. We expect to address these issues in the future. However, we do not expect these assumptions to significantly affect the main results of this work, i.e. the properties of the ejecta and the trends across EOS and mass-ratio in the first 7.5 ms after merger.

Acknowledgements.

We would like to acknowledge helpful discussions with Tim Dietrich, Wyatt Brege and Sergei Ossokine. F.F. gratefully acknowledges support from NASA through grant 80NSSC18K0565, and from the NSF through grant PHY-1806278. H.P. gratefully acknowledges support from the NSERC Canada. L.K. acknowledges support from NSF grant PHY-1606654. M.S. acknowledges support from NSF Grants PHY-170212 and PHY-1708213. L.K. and M.S. also thank the Sherman Fairchild Foundation for their support. Computations were performed on the supercomputer Briaree from the Universite de Montreal, managed by Calcul Quebec and Compute Canada. The operation of these supercomputers is funded by the Canada Foundation for Innovation (CFI), NanoQuebec, RMGA and the Fonds de recherche du Quebec - Nature et Technologie (FRQ-NT). Computations were also performed on the Minerva cluster at the Max-Planck-Institute for Gravitational Physics, and the GPC and Niagara supercomputers at the SciNet HPC Consortium Loken et al. (2010). SciNet is funded by: the Canada Foundation for Innovation; the Government of Ontario; Ontario Research Fund - Research Excellence; and the University of Toronto.