Properties of Neutrino-driven Ejecta from the Remnant of Binary Neutron Star Merger : Purely Radiation Hydrodynamics Case

TL;DR

This study uses general relativistic neutrino radiation hydrodynamics simulations to analyze neutrino-driven ejecta from binary neutron star merger remnants, highlighting the impact of pair-annihilation heating on ejecta dynamics and implications for gamma-ray burst outflows.

Contribution

It provides the first detailed analysis of neutrino-driven ejecta properties in a relativistic framework, emphasizing the role of pair-annihilation heating in merger remnants.

Findings

Pair-annihilation heating doubles the ejecta kinetic energy.

Relativistic outflows are suppressed due to baryon loading.

Conditions for launching gamma-ray burst jets are discussed.

Abstract

We performed general relativistic, long-term, axisymmetric neutrino radiation hydrodynamics simulations for the remnant formed after the binary neutron star merger, which consist of a massive neutron star and a torus surrounding it. As an initial condition, we employ the result derived in a three-dimensional, numerical relativity simulation for the binary neutron star merger. We investigate the properties of neutrino-driven ejecta. Due to the pair-annihilation heating, the dynamics of the neutrino-driven ejecta is significantly modified. The kinetic energy of the ejecta is about two times larger than that in the absence of the pair-annihilation heating. This suggests that the pair-annihilation heating plays an important role in the evolution of the merger remnants. The relativistic outflow, which is required for driving gamma-ray bursts, is not observed because the specific heating rate around the rotational axis is not sufficiently high due to the baryon loading caused by the neutrino-driven ejecta from the massive neutron star. We discuss the condition for launching the relativistic outflow and the nucleosynthesis in the ejecta.

Figures

Tables

Table 1. Table 1 List of the models. M NS subscript 𝑀 NS M_{\rm NS} is the gravitational mass of each NS in the merger simulation of Sekiguchi et al. ( 2015 ) . For Δ ​ x 0 Δ subscript 𝑥 0 \Delta x_{0} , δ 𝛿 \delta , R star subscript 𝑅 star R_{\rm star} , N 𝑁 N , and L 𝐿 L , see the text in Sec.  2.5 .

Model	EOS	$M_{\mathrm{NS}}$	Pair-annihilation	$\nu$-Absorption	$\Delta x_0$	$\delta$	$R_{\mathrm{star}}$	$N$	$L$
		$\left(M_{\odot }\right)$	heating		(m)		(km)		(km)
DD2-135135-On-H (fiducial,high)	DD2	1.35	on	Eq.(40)	150	0.0075	30	951	5440
DD2-135135-Off-H	DD2	1.35	off	Eq.(40)	150	0.0075	30	951	5440
DD2-135135-Iso-H	DD2	1.35	on, isotropic	Eq.(40)	150	0.0075	30	951	5440
DD2-135135-On-H-53	DD2	1.35	on	Eq.(41)	150	0.0075	30	951	5440
DD2-135135-On-M (medium)	DD2	1.35	on	Eq.(40)	200	0.0075	30	871	5790
DD2-135135-On-L (low)	DD2	1.35	on	Eq.(40)	250	0.0075	30	809	5690

Equations

$${\partial }_t\alpha$$

$${\partial }_t{\beta }^i$$

$${\nabla }_{\beta }{T}_{(\mathrm{fluid})}^{\alpha \beta }$$

$${\nabla }_{\beta }{T}_{({\nu }_i)}^{\alpha \beta }$$

$$T_{({\nu }_i)}^{\alpha \beta }=T_{({\nu }_i,\mathrm{T})}^{\alpha \beta }+T_{({\nu }_i,\mathrm{S})}^{\alpha \beta },$$

$$T_{(\mathrm{tot})}^{\alpha \beta }=T^{\alpha \beta }+\sum _i{T}_{({\nu }_i,\mathrm{S})}^{\alpha \beta },$$

$$T_{\alpha \beta }=\rho h{u}_{\alpha }{u}_{\beta }+P{g}_{\alpha \beta }.$$

$${\partial }_t({\rho }_{\ast }\hat{e})+\frac{1}{x}{\partial }_k[x({\rho }_{\ast }\hat{e}{v}^k+W^{-3}P(v^k+{\beta }^k))]$$

$${\partial }_t({\rho }_{\ast }{\hat{u}}_i)+\frac{1}{x}{\partial }_k[x({\rho }_{\ast }{\hat{u}}_i{v}^k+\alpha W^{-3}P{\delta }^k{}_i)]$$

$$\displaystyle\hskip 56.9055pt+\frac{1}{2}W^{-3}S_{kl}\partial_{i}\gamma^{kl}-W^{-3}Q_{\rm(leak)}^{\alpha}\gamma_{\alpha i}\biggr{]},$$

$${\nabla }_{\alpha }(\rho u^{\alpha })$$

$${\nabla }_{\alpha }(\rho Y_e{u}^{\alpha })$$

$${\nabla }_{\alpha }(\rho Y_{{\nu }_i}{u}^{\alpha })$$

$${\partial }_t{\rho }_{\ast }+\frac{1}{x}{\partial }_k(x{\rho }_{\ast }{v}^k)$$

$${\partial }_t({\rho }_{\ast }{Y}_e)+\frac{1}{x}{\partial }_k(x{\rho }_{\ast }{Y}_e{v}^k)$$

$${\partial }_t({\rho }_{\ast }{Y}_{{\nu }_i})+\frac{1}{x}{\partial }_k(x{\rho }_{\ast }{Y}_{{\nu }_i}{v}^k)$$

$$T_{({\nu }_i,\mathrm{S})}{}_{\alpha \beta }=E{n}_{\alpha }{n}_{\beta }+F_{\alpha }{n}_{\beta }+F_{\beta }{n}_{\alpha }+P_{\alpha \beta },$$

$${\partial }_t\tilde{E}+\frac{1}{x}{\partial }_k[x(\alpha {\tilde{F}}^k-{\beta }^k\tilde{E})]=$$

$${\partial }_t{\tilde{F}}_i+\frac{1}{x}{\partial }_k[x(\alpha {\tilde{P}}_i{}^k-{\beta }^k{\tilde{F}}_i)]=$$

$$P^{ij}$$

$${\left(P^{ij}\right)}_{\mathrm{thin}}$$

$${\left(P^{ij}\right)}_{\mathrm{thick}}$$

$$+\frac{1}{w}\left[F^i{V}^j+F^j{V}^i\right]$$

$$+\frac{2F^k{u}_k}{w(2w^2+1)}\left[-w^2{\gamma }^{ij}+V^i{V}^j\right],$$

$$\chi =\frac{3+4f^2}{5+2\sqrt{4-3f^2}}.$$

$$Q_{{\nu }_e}^{\alpha }=$$

$$-Q_{\mathrm{pair},{\nu }_e}^{(+)}{}^{\alpha }-Q_{\mathrm{abs},{\nu }_e}^{(+)}{}^{\alpha },$$

$$Q_{{\bar{\nu }}_e}^{\alpha }=$$

$$-Q_{\mathrm{pair},{\bar{\nu }}_e}^{(+)}{}^{\alpha }-Q_{\mathrm{abs},{\bar{\nu }}_e}^{(+)}{}^{\alpha },$$

$$Q_{{\nu }_x}^{\alpha }=$$

$$-Q_{\mathrm{pair},{\nu }_x}^{(+)}{}^{\alpha }.$$

$$\displaystyle Q^{(+)}_{{\rm pair},\nu_{i}}{}^{\alpha}=\$$

$$\displaystyle=\$$

$$Q_{pair/Iso,{\nu }_i}^{(+)}=$$

$${\mathcal{R}}_e$$

$${\mathcal{R}}_{{\nu }_e}$$

$${\mathcal{R}}_{{\bar{\nu }}_e}$$

$${\mathcal{R}}_{{\nu }_x}$$

$$Q_{\mathrm{abs},{\nu }_e}^{(+)}{}^{\alpha }$$

$$Q_{\mathrm{abs},{\bar{\nu }}_e}^{(+)}{}^{\alpha }$$

$${\mathcal{R}}_{\mathrm{abs},{\nu }_e}$$

Full text

1. Introduction

Binary neutron star (NS-NS) merger is one of the promising sources of gravitational waves for ground-based gravitational-wave detectors, such as advanced LIGO (Abadie et al. 2010), advanced Virgo (Accadia et al. 2011), and KAGRA (Kuroda & LCGT Collaboration 2010). In association with gravitational waves, the merger remnant, which is composed typically of a central compact object (NS or black hole) surrounded by a massive accretion torus, emits a huge amount of neutrinos. Neutrinos emitted could be a major source of various types of phenomena as follows.

First, a large amount of mass could be ejected as sub-relativistic neutrino-driven wind (Dessart et al. 2009; Metzger & Fernández 2014; Perego et al. 2014), in addition to the dynamical mass ejection, which would occur at NS-NS mergers. Then, a substantial amount of radioactive nuclei are likely to be synthesized via the $r$-process (Lattimer & Schramm 1974; Symbalisty & Schramm 1982; Eichler et al. 1989; Rosswog et al. 1999; Freiburghaus et al. 1999; Goriely et al. 2011; Wanajo et al. 2014). The ejecta is heated by the radioactive decay and fission of these heavy elements and would shine at optical and infra-red bands in $1-10$ days after the merger (Li & Paczyński 1998; Kasen et al. 2013; Tanaka & Hotokezaka 2013). Simultaneous detection of this so-called “kilonova” or “macronova” with gravitational waves will significantly improve positional accuracy of gravitational wave sources. The neutrino-driven ejecta could be an additional energy source to the dynamical ejecta.

Second, if a relativistic outflow is launched by neutrino heating, the merger remnant would drive short-duration gamma-ray bursts (SGRBs). The neutrino pair-annihilation could be the engine for driving the extremely relativistic ejecta (Eichler et al. 1989; Meszaros & Rees 1992; Narayan et al. 1992). This is because this process can deposit thermal energy into materials regardless of the baryon density, so that the ejecta could achieve large photon-to-baryon ratio if the baryon density is sufficiently low in the heating region. We should note that there are several other scenarios proposed to drive the relativistic ejecta (Nakar 2007), such as the electromagnetic energy extractions from rotating black holes (Blandford & Znajek 1977) and rapidly rotating pulsars of a strong magnetic field (Usov 1992).

Radiation hydrodynamics simulations of a merger-remnant torus surrounding a central compact object have been performed by many authors (Dessart et al. 2009; Fernandez et al. 2013; Metzger & Fernández 2014; Perego et al. 2014; Fernández et al. 2015; Just et al. 2015, 2016; Foucart et al. 2016). Recently, simulations of black hole-torus system approximately taking the pair-annihilation heating into account were performed (Just et al. 2016). It was concluded that the system is unlikely to drive SGRBs because the neutrino heating cannot overcome the baryon loading.

Recent numerical relativity simulations of the NS-NS mergers with finite-temperature equations of state (EOS) that can support the observed $\sim 2M_{\odot}$ NSs (Demorest et al. 2010; Antoniadis et al. 2013) have suggested that a massive neutron star (MNS) is likely to be formed as the remnant of the merger if the total binary mass is not extremely high ($\lesssim 2.8M_{\odot}$) (Shibata et al. 2005; Shibata & Taniguchi 2006; Kiuchi et al. 2009, 2010; Sekiguchi et al. 2011; Hotokezaka et al. 2011, 2013; Kaplan et al. 2014; Takami et al. 2015; Dietrich et al. 2015). The neutrino emissivity of the MNS-torus system is quite different from that of the black hole-torus system. The MNS itself emits a large amount of neutrinos. In addition, the torus matter does not accrete into the black hole, but stops at the surface of the MNS. Such an inner region of the torus could contribute significantly to the neutrino emission. Therefore, if the lifetime of MNSs is sufficiently long, a huge influence of neutrinos on the ejecta is expected (Metzger & Fernández 2014; Lippuner et al. 2017). Recent numerical relativity simulations for the NS-NS mergers (Sekiguchi et al. 2015, 2016) have indeed shown that the effect of neutrino irradiation is significant near the rotational axis. However, their simulation time is not so long as the accretion time of the torus ($\gtrsim 100$ ms). Hence, to explore the entire effects of neutrinos a longer-term simulation is required.

In Richers et al. (2015), the pair-annihilation heating rate in a MNS-torus system was calculated with a Monte-Carlo method using snapshots of the pseudo-Newtonian hydrodynamics simulations (Metzger & Fernández 2014) as their background configurations. Perego et al. (2017) also calculated the heating rate with a ray-tracing method as a post process of the Newtonian hydrodynamics simulations for the NS-NS merger (Perego et al. 2014). It was suggested that the total energy deposited by the pair-annihilation heating is not large enough to account for the typical energy of observed SGRB events, although the presence of the MNS does increase the pair-annihilation rate. Perego et al. (2017) concluded that much larger neutrino luminosity is required to explain most of the observed SGRBs. Compared with the Newtonian simulation, however, the merger becomes more violent in simulations with general relativity because of the stronger gravitational effects. Hence, the temperature of the merger remnant becomes higher than that with Newtonian gravity. Therefore, their Newtonian simulation may underestimate the neutrino luminosity. Moreover, Perego et al. (2017) employed a stiff EOS, referred to as TM1 (Shen et al. 1998a, b; Hempel et al. 2012). For the stiff EOS, the remnant MNS is less compact, and hence, the neutrino luminosity is lower than that with softer EOSs. The neutrino pair-annihilation heating could be more efficient in the general relativistic simulation with softer EOS.

Motivated by the above considerations, in this paper, we investigate neutrino-driven ejecta from the MNS-torus system. We perform long-term fully general relativistic simulations with an approximate neutrino transport in axial symmetry using the initial condition based on a three-dimensional numerical relativity merger simulation. We also investigate the effects of the neutrino-antineutrino pair-annihilation heating, which is not taken into account in our previous studies (Sekiguchi et al. 2015, 2016). The outline of this paper is as follows. A brief description of the simulation setting and the initial condition which we employ are given in Sec. 2. Then, the results of the present simulations are shown in Sec. 3. Finally, discussions and the conclusions are given in Secs. 4 and 5, respectively. Throughout this paper, we use the units of $c=1=G$, where $c$ and $G$ are the speed of light and the gravitational constant, respectively.

2. Method

2.1. Einstein’s Equation

In our simulation, Einstein’s equation is solved in a Baumgarte-Shapiro-Shibata-Nakamura (BSSN)-puncture formalism (Shibata & Nakamura 1995; Baumgarte & Shapiro 1999; Campanelli et al. 2006; Baker et al. 2006). In our formalism, we evolve the conformal three-metric $\tilde{\gamma}_{ij}=\gamma^{-1/3}\gamma_{ij}$, the conformal factor $W=\gamma^{-1/6}$ (Marronetti et al. 2008), the trace of the extrinsic curvature $K=\gamma^{ij}K_{ij}$, the trace-free part of the extrinsic curvature $\tilde{A}_{ij}=\gamma^{-1/3}(K_{ij}-\gamma_{ij}K/3)$, and the auxiliary variable $F_{i}=\delta^{jk}\partial_{j}\tilde{\gamma}_{ki}$. Here, $\gamma_{\alpha\beta}=g_{\alpha\beta}+n_{\alpha}n_{\beta}$ is the induced metric, where $g_{\alpha\beta}$ and $n_{\alpha}$ are the spacetime metric and the timelike unit vector normal to the time slice, respectively, $\gamma=\det(\gamma_{ij})$, and $\det(\tilde{\gamma}_{ij})$ is assumed to be unity. We adopt the so-called cartoon method (Alcubierre et al. 2001; Shibata 2003a) to impose axially symmetric conditions for the geometric quantities. We evaluate the spatial derivative by a fourth-order central finite differencing. We add sixth order Kreiss-Oliger-type dissipation terms (Brügmann et al. 2008) to the evolution equations of the geometrical variables for the stability of numerical evolution. For the gauge condition, we adopt dynamical lapse (Alcubierre et al. 2003) and shift (Shibata 2003b) conditions as

[TABLE]

where $\alpha$ and $\beta_{i}$ are the lapse function and the shift vector, respectively, and $\Delta t$ is the time-step interval.

2.2. Neutrino-radiation Hydrodynamics Equations

For exploring the neutrino-driven outflow from the NS-NS merger remnant, we need to solve hydrodynamics equations together with neutrino radiation transfer equations. In this work, we divide neutrinos into “streaming” and “trapped” components. For streaming-neutrinos, we employ the so-called M1-closure scheme based on Thorne’s moment formalism (Thorne 1981; Shibata et al. 2011). On the other hand, for trapped-neutrinos, we employ a leakage-based scheme (for the leakage scheme, see, e.g., Cooperstein 1988). In the following, we give a brief description of our formulation. In this work, we consider three species of neutrinos: electron neutrinos $\nu_{e}$, electron antineutrinos $\bar{\nu}_{e}$, and the other neutrino species $\nu_{x}$, which represents all of the heavy-lepton neutrinos $\nu_{\mu}$, $\bar{\nu}_{\mu}$, $\nu_{\tau}$, and $\bar{\nu}_{\tau}$.

The basic equation of the fluid is derived from the energy-momentum conservation equation. On the other hand, the basic equation of neutrinos is derived from the energy-integrated, first moment of Boltzmann’s equation. They are written, respectively, as

[TABLE]

where $T_{\rm(fluid)}^{\alpha\beta}$ and $T_{(\nu_{i})}^{\alpha\beta}$ are the energy-momentum tensors of the fluid and the $i$-th species of neutrinos, respectively, and $Q^{\alpha}_{\nu_{i}}$ is the energy-momentum source term of the $i$-th species of neutrinos, which is determined by weak-interaction processes. Here, the total energy-momentum tensor is $T^{\alpha\beta}_{\rm(tot)}=T_{\rm(fluid)}^{\alpha\beta}+\sum_{i}T_{(\nu_{i})}^{\alpha\beta}$, which satisfies $\nabla_{\beta}T^{\alpha\beta}_{\rm(tot)}=0$.

In our scheme, the energy-momentum tensor of neutrinos is phenomenologically decomposed as

[TABLE]

where $T_{(\nu_{i},{\rm T})}^{\alpha\beta}$ and $T_{(\nu_{i},{\rm S})}^{\alpha\beta}$ are the contributions of trapped- and streaming- neutrinos, respectively. We assume that a part of neutrinos produced at a rate $Q^{\alpha}_{\nu_{i}}$ becomes streaming-neutrinos at a leakage rate $Q_{({\rm leak})\nu_{i}}^{\alpha}$ (for the detail of this term, see Sekiguchi et al. 2012), and the other part becomes trapped-neutrinos at a rate $Q^{\alpha}_{\nu_{i}}-Q^{\alpha}_{({\rm leak})\nu_{i}}$.

We assume that trapped-neutrinos are tightly coupled with the fluid. Then the energy-momentum tensor of trapped-neutrinos can be written by that of the fluid of relativistic Fermi particles, and thus, we decompose the total energy-momentum tensor as

[TABLE]

where $T^{\alpha\beta}=T_{\rm(fluid)}^{\alpha\beta}+\sum_{i}T_{(\nu_{i},{\rm T})}^{\alpha\beta}$ is the energy-momentum tensor composed of the sum of the fluid and trapped-neutrinos. It is written in the form

[TABLE]

Here $\rho$ is the baryon rest-mass density, $u^{\alpha}$ is the four-velocity of the fluid, $h=1+\varepsilon+P/\rho$ is the specific enthalpy, $P$ is the pressure, and $\varepsilon$ is the specific internal energy. $P$ and $\varepsilon$ contain the contributions from baryons, electrons, positrons, and trapped-neutrinos.

We evolve the fluid variables in cylindrical coordinates $(R,\varphi,z)$. Hereafter, we regard the $y=0$ plane of Cartesian coordinates $(x,y,z)$ as a meridional ($\varphi=0$) plane of cylindrical coordinates $(R,\varphi,z)$. The Euler and energy equations can be written, respectively, as {widetext}

[TABLE]

where $\rho_{*}\equiv W^{-3}w\rho$, $\hat{u}_{i}\equiv hu_{i}$, $v^{i}\equiv u^{i}/u^{t}$, $\hat{e}\equiv hw-P/(w\rho)$, $S_{ij}\equiv T_{\alpha\beta}\gamma^{\alpha}{}_{i}\gamma^{\beta}{}_{j}=\rho hu_{i}u_{j}+P\gamma_{ij}$, $w\equiv-n_{\alpha}u^{\alpha}=\sqrt{1+\gamma^{ij}u_{i}u_{j}}$ is the Lorentz factor of the fluid, and $Q^{\alpha}_{\rm(leak)}$ is the sum of leakage source terms of all neutrino species.

The number density of baryons, electrons, and trapped-neutrinos is evolved by

[TABLE]

where ${\cal R}_{e}$ and ${\cal R}_{\nu_{i}}$ denote the source terms for the evolution of the lepton number density due to the weak interactions. Here we assumed that a part of neutrinos produced at a rate ${\cal R}_{\nu_{i}}$ becomes trapped-neutrinos at a rate ${\cal R}_{\nu_{i}}-{\cal R}_{({\rm leak})\nu_{i}}$, where ${\cal R}_{({\rm leak})\nu_{i}}$ is the emission rate of neutrinos of the flavor $i$ due to the diffusion (for the details, see Sekiguchi et al. 2012). Equations. (10)–(12) are written in the following forms, respectively;

[TABLE]

[TABLE]

For streaming-neutrinos, we decompose their energy-momentum tensor as

[TABLE]

where $E\equiv T_{({\rm S})}{}_{\mu\nu}n^{\mu}n^{\nu}$, $F_{\alpha}\equiv-T_{({\rm S})}{}_{\mu\nu}\gamma^{\mu}{}_{\alpha}n^{\nu}$, and $P_{\alpha\beta}\equiv T_{({\rm S})}{}_{\mu\nu}\gamma^{\mu}{}_{\alpha}\gamma^{\nu}{}_{\beta}$ are the energy density, energy flux, and spatial stress tensor of streaming-neutrinos, respectively. $F_{\alpha}$ and $P_{\alpha\beta}$ satisfy $F_{\alpha}n^{\alpha}=0=P_{\alpha\beta}n^{\alpha}$. In the above, we omitted the subscripts of the species of neutrinos. The evolution equations for $E$ and $F_{i}$ are written as {widetext}

[TABLE]

where $\tilde{E}\equiv W^{-3}E$, $\tilde{F}_{i}\equiv W^{-3}F_{i}$, and $\tilde{P}^{ij}\equiv W^{-3}P^{ij}$. In this formalism, we have to impose a condition to determine $P^{ij}$. For this, we use the so-called M1-closure relation (Levermore 1984; González et al. 2007). In this closure relation, the spatial stress tensor is written as

[TABLE]

where $\left(P^{ij}\right)_{\rm thin}$ and $\left(P^{ij}\right)_{\rm thick}$ are the spatial stress tensor in the optically thin and thick limits, respectively, and are written as

[TABLE]

where $V^{i}=\gamma^{ij}u_{j}$ (Shibata et al. 2011). $\chi$ is the so-called variable Eddington factor, which is a function of a normalized flux $f^{2}=h_{\alpha\beta}H^{\alpha}H^{\beta}/J^{2}$ and we choose the following form:

[TABLE]

Here we defined the energy density $J$, energy flux $H^{\alpha}$, and spatial stress-tensor $L^{\alpha\beta}$ of streaming-neutrinos in the fluid rest frame by $J=T_{({\rm S})}^{\alpha\beta}u_{\alpha}u_{\beta}$, $H^{\alpha}=-T_{({\rm S})}^{\gamma\beta}u_{\beta}h^{\alpha}{}_{\gamma}$, and $L^{\alpha\beta}=T_{({\rm S})}^{\gamma\delta}h^{\alpha}{}_{\gamma}h^{\beta}{}_{\delta}$, respectively, where $h_{\alpha\beta}=g_{\alpha\beta}+u_{\alpha}u_{\beta}$ is a projection operator.

2.3. Microphysics

2.3.1 Equation of State

In this paper, for the nuclear equation of state (EOS), we adopt a tabulated EOS referred to as DD2 (Banik et al. 2014), which is the same EOS as one of those used in three-dimensional simulations for the NS-NS merger (Sekiguchi et al. 2015). The remnant of the NS-NS merger of the typical total mass ($\sim 2.7M_{\odot}$) is a long-lived MNS, and hence, we can explore the long-term phenomena associated with the MNS surrounded by a torus.

The original table of this EOS covers only ranges of density $10^{3.2}-10^{16.2}\ {\rm g\ cm^{-3}}$ and temperature $10^{-1}-10^{2.2}$ MeV, but the density and the temperature of the ejecta become lower than the lowest values of the table. For treating such low density and temperature regions, we extend the table to ranges with $10^{-0.8}-10^{16.2}\ {\rm g\ cm^{-3}}$ and $10^{-3}-10^{2.2}$ MeV using an EOS by Timmes & Swesty (2000). This EOS includes contributions of nucleons, heavy nuclei, electrons, positrons and photons to the pressure and the internal energy.

2.3.2 Weak Interaction

The source terms of Eqs. (3) and (4) are

[TABLE]

Here, $Q_{{\rm brems},{\nu_{i}}}^{(-)}$, $Q_{{\rm pair},{\nu_{i}}}^{(-)}$, $Q_{{\rm plasm},{\nu_{i}}}^{(-)}$, $Q_{\rm ec}^{(-)}$, and $Q_{\rm pc}^{(-)}$ are the neutrino cooling (emission) rates due to the nucleon-nucleon Bremsstrahlung, the electron-positron pair-annihilation, the plasmon decay, the electron-capture, and the positron-capture processes, respectively. $Q_{{\rm pair},{\nu_{i}}}^{(+)}{}^{\alpha}$, $Q_{{\rm abs},\nu_{e}}^{(+)}{}^{\alpha}$, and $Q_{{\rm abs},\bar{\nu}_{e}}^{(+)}{}^{\alpha}$ are the matter-heating source terms due to the neutrino-antineutrino pair-annihilation, the electron-neutrino and electron-antineutrino absorption processes, respectively. For electron and positron capture processes, we use the rate in Fuller et al. (1985). For pair-production processes, we use the rates in Cooperstein et al. (1986) for electron-positron pair-annihilation, Ruffert et al. (1996) for plasmon-decay, and Burrows et al. (2006) for nucleon-nucleon Bremsstrahlung, respectively. Except for the neutrino pair-annihilation process, the rates adopted in the present simulations are the same as those adopted in the simulations for three-dimensional NS-NS mergers (Sekiguchi et al. 2015). The explicit forms of the cooling rates ($Q_{{\rm brems},{\nu_{i}}}^{(-)}$, $Q_{{\rm pair},{\nu_{i}}}^{(-)}$, $Q_{{\rm plasm},{\nu_{i}}}^{(-)}$, $Q_{\rm pc}^{(-)}$, and $Q_{\rm ec}^{(-)}$) are found in Sekiguchi et al. (2012).

In our moment formalism, the source term due to the neutrino pair-annihilation process is written (in $\hbar=1$ unit) as {widetext}

[TABLE]

where $\left<\omega\right>_{\rm pair}$ is the average energy of streaming-neutrinos that annihilate, and the quantities with bar imply those of antineutrinos. The exponential factor in the first line of Eq. (26) is introduced to suppress the heating rate in the optically thick region because the pair-annihilation heating balances with electron-positron pair-annihilation cooling in such a region. $G_{\rm F}$ is the Fermi coupling constant and $C_{\nu_{i}\bar{\nu}_{i}}^{\rm pair}$ is written by the Weinberg angle $\theta_{\rm W}$ as $C_{\nu_{i}\bar{\nu}_{i}}^{\rm pair}=1\pm 4\sin^{2}\theta_{\rm W}+8\sin^{4}\theta_{\rm W}$, where the plus and minus signs denote neutrinos of electron and heavy lepton types, respectively. We use the value $\sin^{2}\theta_{\rm W}\approx 0.2319$. The detailed derivation of Eq. (26) is described in Appendix A.

Because we assume an analytic closure relation, the evaluation of the pair-annihilation rate in our simulation could be significantly different from the true rate. In this paper, the uncertainty in the pair-annihilation heating rate is explored by performing a simulation in which the pair-annihilation heating rate is calculated assuming the isotropic momentum-space angular distribution for neutrinos. Under this assumption, the pair-annihilation source term $Q^{(+)}_{{\rm pair/Iso},\nu_{i}}\ (i=e,x)$ is written as

[TABLE]

In this model, the pair-annihilation heating rate may be incorporated in an optimistic manner (see Sec. 4.4 for discussion).

For the source terms of Eqs. (11) and (12), we employ the following equations:

[TABLE]

For these rates except for the electron (anti)neutrino absorption processes, we employ the rates in Fuller et al. (1985), Cooperstein et al. (1986), Ruffert et al. (1996), and Burrows et al. (2006), which are the same as in Sekiguchi et al. (2012).

The source terms due to the electron (anti)neutrino absorption processes are written as

[TABLE]

where $\left<\omega\right>_{\rm abs}$ is the average energy of streaming-neutrinos, $m_{u}$ is the atomic mass unit, $g_{\rm A}\approx 1.26$ is the axial vector coupling strength, and $X_{n}$ and $X_{p}$ are the mass fractions of neutrons and protons, respectively. For the same reason as in the pair-annihilation heating, the exponential factors in Eqs. (32)–(35) are introduced to suppress these rates in the optically thick region, in which these rates balance with those of their inverse processes. The detailed derivation for these rates is described in Appendix B.

2.3.3 Average Energy Estimation

For evaluating the average energy of streaming-neutrinos $\left<\omega\right>_{\rm pair}$ in our energy-integrated radiation transfer scheme, we assume that streaming-neutrinos have the Fermi-Dirac-type energy distribution, i.e.,

[TABLE]

Then the energy density of neutrinos in the fluid rest frame $J$, the “temperature” of the neutrino $T_{\nu}$, and the normalized “chemical potential” of streaming-neutrinos $\eta_{\nu}$ satisfy the relation

[TABLE]

where $F_{k}(\eta)$ is the relativistic Fermi integral of order $k$ defined by

[TABLE]

We assumed that the “temperature” of streaming-neutrinos $T_{\nu}$ is the same as the local temperature of the matter $T$. Then we obtain the normalized “chemical potential” of neutrinos, $\eta_{\nu}$, by solving Eq. (37). The average energy of streaming-neutrinos is then evaluated as

[TABLE]

This average energy would be smaller than the real one because the temperature of neutrinos is in reality comparable to the temperature of their emission region, which is usually higher than that of the free-streaming region of neutrinos. For this reason, we should keep in mind that the pair-annihilation heating rate would be evaluated by our scheme in a conservative manner.

In the same way as Eq. (39), we evaluate $\left<\omega\right>_{\rm abs}$ as

[TABLE]

where we simply derive the average energy of neutrinos assuming the Fermi-Dirac-type energy distribution.

We have to keep in mind that there are several uncertainties in the average energy estimation in our energy-integrated neutrino transport scheme. In particular, the neutrino absorption reactions have strong energy dependence, so that the change of the average energy would affect the dynamics of the neutrino-driven ejecta. For checking such uncertainties, we perform a simulation in which the average energy of neutrinos are estimated by

[TABLE]

instead of Eq. (40). This is derived when we take the energy dependence of the absorption reaction into account, and factor out the neutrino energy density of Eq. (37). We give a detail for deriving Eq. (41) in Appendix B.

2.4. Initial Condition

We mapped three-dimensional simulation data for a NS-NS merger (Sekiguchi et al. 2015) into two-dimensional data as the initial condition for the axisymmetric simulation. Although the merger remnant has a non-axisymmetric structure in an early phase after the merger, it gradually relaxes to a nearly axisymmetric structure. It also approaches a quasi-stationary state. In this work, we employ the remnant at about 50 ms after the onset of the merger as an initial condition. In such a phase, the merger remnant is approximately in an axisymmetric quasi-stationary state.

We generate axisymmetric data from three-dimensional data by taking average for $\rho_{*}$, $\hat{e}$, $\hat{u}_{i}$, $Y_{e}$, and $Y_{\nu_{i}}$ by

[TABLE]

where the quantities with superscript (3D) imply those taken from the corresponding three-dimensional data. After mapping, we solve the Hamiltonian and momentum constraint equations assuming the conformal flatness. Employing the conformal flatness is reasonable because the values of $\tilde{\gamma}_{ij}-\delta_{ij}$ are fairly small at $\sim$50 ms after the merger; specifically the absolute magnitude of all the components is smaller than 0.02.

The configurations of initial thermodynamical quantities are shown in Fig. 1. The MNS is surrounded by a massive torus of mass $\sim 0.2M_{\odot}$. Here we defined the mass of the torus as the baryonic mass in the region where the density is lower than $10^{13}\ {\rm g\ cm^{-3}}$. The maximum temperature is $\sim$30 MeV near the equatorial surface region of the MNS. The temperature of the inner region of the torus is $\sim 5-7$ MeV. The effective neutrino emissivity is highest on the polar surface of the MNS and in the inner region of the torus.

We summarize the models for our simulations in Table 1. We refer to the fiducial model as DD2-135135-On-H. For understanding the effects of the pair-annihilation heating, we also perform a simulation without the neutrino pair-annihilation (the model DD2-135135-Off-H). In addition to these models, we perform simulations of the models DD2-135135-Iso-H and DD2-135135-On-H-53 to check the uncertainty of the heating rates due to the neutrino pair-annihilation and the neutrino absorption processes. In the former model, the pair-annihilation heating rate is calculated by Eq. (27) assuming the isotropic momentum-space angular distribution for neutrinos. The latter is the simulation in which the average energy used in the neutrino absorption heating process is estimated by Eq. (41) instead of Eq. (40).

2.5. Grid Setting

We adopt a non-uniform grid in which the grid spacing is increased according to the rule

[TABLE]

where $\Delta x_{j}\equiv x_{j+1}-x_{j}$, $\Delta z_{l}\equiv z_{l+1}-z_{l}$, $\delta$ and $R_{\rm star}$ are constants, and $0\leq j\leq N$ ($0\leq l\leq N$). Here, $N+1$ is the total grid number for one direction. In this grid, a uniform grid with the grid spacing $\Delta x_{0}$ is adopted in the inner region $0\leq x\leq R_{\rm star}$ and $0\leq z\leq R_{\rm star}$ to resolve the MNS. The convergence of numerical results is examined by performing simulations with coarser grid resolutions (DD2-135135-On-M and DD2-135135-On-L, see Table 1). The values of the innermost grid spacing $\Delta x_{0}$, the constant $\delta$, the size of the region in which the uniform grid is adopted $R_{\rm star}$, the grid number $N$, and the size of computational domain $L$ are also tabulated in Table 1.

3. Result

3.1. Dynamics of the System

First, we briefly summarize the dynamics of the fiducial model DD2-135135-On-H. For an early phase of evolution ($t\lesssim 50$ ms, where $t$ is the time after mapping from the 3D simulation), a strong outflow is launched in the vicinity of the rotational axis. This is driven primarily by the neutrino pair-annihilation heating in the dilute matter. The neutrino absorption heating is subdominant in the fiducial model. The evolution of the rest-mass density, the velocity vector, and the specific heating rate is displayed in the left panel of Fig. 2. This clearly shows that the outflow is launched around the rotational axis. The maximum velocity of the ejecta is $\sim 0.5\ c$. We also display snapshots of the rest-mass density, the velocity vector, and the specific heating rate in the absence of the pair-annihilation heating (DD2-135135-Off-H) in the right panel of Fig. 2. For this case, an outflow is launched only by the neutrino absorption heating, and the ejecta velocity is $0.1-0.2\ c$. As found in the right halves of the both panels in Fig. 2, the specific heating rate is considerably different between the two simulations. The specific heating rate in the polar region of the model DD2-135135-On-H is much higher than that of the model DD2-135135-Off-H because the pair-annihilation heating dominates the other heating processes (i.e., the absorption of electron neutrinos and electron antineutrinos on nucleons). This illustrates that the high-speed ejecta is launched by the strong pair-annihilation heating.

At $t\sim$ 50 ms, the maximum velocity of the ejecta goes down to $\sim 0.2\ c$ even in the presence of the pair-annihilation heating. This is because the neutrino pair-annihilation heating rate decreases with time as found in the middle panel of Fig. 2 (and see also Fig. 3). The decrease of the pair-annihilation heating rate is caused by the decrease of the neutrino luminosity.

3.2. Neutrino Luminosity, Pair-annihilation Rate, and Efficiency

The top panel of Fig. 3 shows the luminosity curves of individual neutrinos as functions of time. At the beginning of the simulation, the luminosity of electron-type neutrinos is $\sim 10^{53}\ {\rm erg\ s^{-1}}$, and an order of magnitude larger than those by proto-neutron stars formed during the typical supernova explosions (Fischer et al. 2010). In 100 ms, the luminosity decreases rapidly, because the temperature of the torus, which is a strong emitter of electron-type neutrinos, decreases from $\sim 8$ MeV to $\sim 5$ MeV (see Fig. 4). For $t\gtrsim 300$ ms, the luminosity settles into nearly constant in time as $\sim 10^{52}\ {\rm erg\ s^{-1}}$. In this phase, the primary source of the neutrino emission is the MNS and the neutrino emissivity of the torus is much smaller than that of the MNS.

We define the total pair-annihilation heating rate outside the neutrinosphere by

[TABLE]

where $d\Sigma_{\alpha}=d^{3}x\sqrt{\gamma}n_{\alpha}$ is the three-dimensional volume element on spatial hyper-surfaces.

The bottom panel of Fig. 3 shows $Q_{\rm pair,tot}$ (solid curve) as a function of time. In the fiducial model, $Q_{\rm pair,tot}$ is $\sim 10^{50}\ {\rm erg\ s^{-1}}$ at $t\sim 50$ ms, and it decreases with time as the neutrino luminosity decreases. For $t\agt 300$ ms, it settles down to $\sim 10^{49}\ {\rm erg\ s^{-1}}$.

In the same panel, we also plot the efficiency of the pair-annihilation heating (dashed curves)

[TABLE]

as a function of time (Ruffert et al. 1997; Setiawan et al. 2004). Since $Q_{\rm pair,tot}\propto L_{\nu}^{2}$, the efficiency is proportional to the neutrino luminosity. Therefore, it decreases with time and eventually settles to $\sim 10^{-4}$ in the quasi-stationary phase.

This efficiency is consistent with those obtained with a Monte-Carlo method in the central NS case (Richers et al. 2015) and obtained with a multi-energy M1-scheme in the central black hole case (Just et al. 2016) when the individual neutrino luminosity is $\sim 10^{52}\ {\rm erg\ s^{-1}}$.

3.3. Properties of Neutrino-driven Ejecta

We define the baryonic mass, the total energy, and the internal energy of the ejecta by

[TABLE]

where we supposed that fluid elements with $|u_{t}|>1$ are gravitationally unbound. Using these quantities, we define the kinetic energy of the ejecta by

[TABLE]

We evaluate these quantities in the cylindrical region of $x\leq 2000$ km and $-2000$ km $\leq z\leq 2000$ km. We also take into account the mass and energy gone outside this inner region by integrating the fluxes at the boundary of this region and adding them to Eqs. (51)–(53).

We plot the results in Fig. 5. The top panel shows that the ejecta mass increases for $t\lesssim 100$ ms, but the growth rate decreases with time after that. The ejecta mass reaches about $8\times 10^{-4}M_{\odot}$ finally. The kinetic energy, shown in the middle panel, has the same trend as the ejecta mass, and its final value is about $1\times 10^{49}$ erg. These values are comparable to those of the dynamical ejecta for the models with the DD2 EOS, which are $\sim 10^{-3}M_{\odot}$ and $\sim 2\times 10^{49}$ erg (Sekiguchi et al. 2015).

Using these quantities, we estimate the average velocity of the ejecta $V_{\rm ej}$ by

[TABLE]

Note that we do not consider the contribution from a component that was already unbound at the beginning of the simulation, i.e., we only take into account the neutrino-driven ejecta. We show the evolution of $V_{\rm ej}$ in the bottom panel of Fig. 5. The average velocity of the ejecta is initially about $0.18\ c$, but decreases to $0.13\ c$ in the quasi-stationary phase. The initial enhancement for $V_{\rm ej}$ is induced by the strong pair-annihilation heating.

In Fig. 5, we also plot the results for the model DD2-135135-Iso-H, in which the pair-annihilation heating rate is taken into account assuming the isotropic angular distribution for neutrinos. The ejecta mass and the kinetic energy for the model DD2-135135-Iso-H are about 1.5 and 4 times larger than those for the fiducial model, respectively. Hence, the averaged velocity of the ejecta for the model DD2-135135-Iso-H is larger as $\sim 0.2\ c$ than that in the fiducial model. However, the relativistic outflow is not observed even in this model. The reason would be that the specific heating rate is not sufficiently high even in this model.

Next, we pay attention to the mass distributions of the electron fraction and the specific entropy of the ejecta. Figure 6 shows the mass histogram of the ejecta $Y_{e}$ and entropy at $t=300$ ms, at which the ejecta mass becomes approximately constant in time (see Fig. 5). We find from the left panel that the ejecta are mildly neutron-rich, and have relatively high electron fraction $Y_{e}>0.25$, with the typical value $\sim 0.4$. This property depends only weakly on the presence of the pair-annihilation heating. This is because the fluid velocity just above the MNS ($z\approx 10$ km) is very small $v\sim 0.01\ c$ even in the presence of the pair-annihilation heating (see the top panel of Fig. 7).

Hence the neutrino absorption timescale in the polar region, which is estimated using the average energy of the neutrinos $\left<\omega\right>$ (in $\hbar=1$ unit) by

[TABLE]

is shorter than the time for the fluid element needed to escape from the polar region

[TABLE]

That is, the electron fraction of the ejecta achieves an equilibrium value (e.g., see Eq. (77) in Qian & Woosley 1996) soon after the ejecta is launched from the MNS.

From the right panel of Fig. 6, we find that the typical specific entropy of the ejecta is $\sim 10\ k_{\rm B}$. We also find that a small fraction of the ejecta has very high specific entropy; the highest value is $500\ k_{\rm B}$. The pair-annihilation heating process generates the high-entropy ejecta because this process heats up the material regardless of the baryon density, and hence, a large amount of thermal energy can be injected in the low baryon density region. Implication of these results on the nucleosynthesis is discussed in Sec. 4.2.

3.4. Dependence of Ejecta Properties on Average-Energy Estimation of Neutrinos

Here we discuss the dependence of ejecta properties on the methods of estimating the average energy of neutrinos. The orange curves in Figs. 5 and 6 denote the mass, kinetic energy, average velocity, and mass histogram of the ejecta for the model DD2-135135-On-H-53, in which the average energy is estimated by Eq. (41). We find that the ejecta mass is enhanced, and that the high-entropy component of $s\gtrsim 200\ k_{\rm B}$ is absent for the model DD2-135135-On-H-53.

In the vicinity of the MNS and the torus, the neutrino absorption process dominates the heating rate due to the high density. Thus, the increase of the neutrino absorption heating rate enhances the mass ejection from the system. In general, $\sqrt{F_{5}(\eta_{\nu})/F_{3}(\eta_{\nu})}$ is larger than $F_{3}(\eta_{\nu})/F_{2}(\eta_{\nu})$, and hence, the neutrino absorption heating rate estimated by Eq. (41) are larger than those by Eq. (40) in the entire region. As a result, the ejecta mass is also larger for the model DD2-135135-On-H-53 as found in Fig. 5.

Because of the enhancement of the mass ejection from the MNS, the mass density around the rotational axis increases, and the specific heating rate decreases in the pole region where the pair-annihilation heating is efficient. As a result, the highest end of the specific entropy of the ejecta for the model DD2-135135-On-H-53 does not become as large as that for the model DD2-135135-On-H.

The average energy is thus a key quantity for launching the relativistic outflow because the baryon-loading would be more significant if the neutrino absorption heating is larger. We have to keep in mind that the condition for the relativistic outflow discussed in Sec 4.1 will depend on the estimation of the average energy of neutrinos.

To rigorously discuss the average energy of neutrinos and the neutrino absorption heating rates, a simulation with energy-dependent radiation transfer scheme is needed. Although such simulation is beyond scope of our present work, we have to perform it in the future in order to derive more convincing numerical results.

3.5. Convergence Test

In Figs. 3 and 5, we plot the results for the lower resolution models DD2-135135-On-M (solid) and DD2-135135-On-L (dashed) by the thin curves. For $t\lesssim 100$ ms, the time evolution of the neutrino luminosity depends weakly on the grid resolution. The disagreement in the results among the three different resolution models is $\lesssim$5 % , so that a convergence is reasonably achieved during this time. On the other hand, the convergence becomes poor for the late time. The difference in the luminosity of electron antineutrinos and heavy-lepton neutrinos between the highest- and lowest-resolution models is $\sim$40 % at $t=400$ ms, while the luminosity of electron neutrinos does not depend on the resolution significantly.

The resolution dependence of the average energy of neutrinos shows the same trend as that of the luminosity. The increases of the average energy of electron antineutrinos and heavy-lepton neutrinos in time become weaker in higher resolution models, and thus, we can only discuss the upper limit of the average energy of these neutrino species. The possible reason for this behavior is that the density gradient at the surface of the MNS, around which neutrinos are most significantly emitted, becomes steeper at that time, and hence, the diffusion of the neutrino emission is not accurately resolved with the low resolution. Taking this resolution dependence into account, we may conclude that $L_{\bar{\nu}_{e}}\lesssim 10^{52}\ {\rm erg\ s^{-1}}$, $L_{\nu_{x}}\lesssim 7\times 10^{51}\ {\rm erg\ s^{-1}}$, $\left<\omega_{\bar{\nu}_{e}}\right>_{\rm em}\lesssim 13\ {\rm MeV}$, and $\left<\omega_{\nu_{x}}\right>_{\rm em}\lesssim 20\ {\rm MeV}$ at $t=400$ ms.

For the ejected mass and its kinetic energy (see Fig. 5), the convergence is better achieved. The difference between the highest- and lowest-resolution models is $\sim$ 3 % and $\sim$ 2 % for the mass and the kinetic energy of the ejecta, respectively.

4. Discussion

4.1. Can Neutrinos Drive SGRBs?

As mentioned in Sec. 1, the neutrino pair-annihilation could be a driving force for launching a relativistic ejecta, which could be SGRBs. In the results of our fiducial model (DD2-135135-On-H), the kinetic energy of the ejecta ($\sim 10^{49}\ {\rm erg}$) is smaller than the typical energy (the sum of the energy of gamma-rays in prompt emission and the kinetic energy of the blast wave) of SGRBs of $\sim 10^{50}\ {\rm erg}$ (Fong et al. 2015). Moreover, no relativistic outflow is observed. The primary reason for this is the presence of high-density baryon in the polar region of the MNS. Even in the model in which the pair-annihilation heating is optimistically evaluated (DD2-135135-Iso-H), the kinetic energy of the ejecta ($\sim 4\times 10^{49}\ {\rm erg}$) is still lower than $\sim 10^{50}\ {\rm erg}$, and the Lorentz factor of the ejecta is also small ($\Gamma\lesssim 1.3$).

We can estimate the terminal Lorentz factor of the ejecta as

[TABLE]

where $Q$ is the heating rate, and $\tau_{\rm heat}\sim r/v$ is the heating timescale defined by the length scale of the heating region divided by the fluid velocity of the region. Here, as fiducial values, we take $Q/\rho\sim 10^{24}\ {\rm erg\ g^{-1}\ s^{-1}}$ and $\tau_{\rm heat}\sim 30\ {\rm km}/0.1\ c\approx 1$ ms from the results at $t=20$ ms shown in Fig. 7. Note that we assumed that the internal energy deposited in the material due to heating processes would be totally transformed into the kinetic energy. The small Lorentz factor in our simulation ($\Gamma\lesssim 1.3$) is consistent with the estimate of Eq. (58).

The Lorentz factor of $\sim 100$ could be obtained if the specific heating rate is $Q/\rho\sim 10^{26}\ {\rm erg\ g^{-1}\ s^{-1}}$ over a length scale of $\sim 30$ km. We plot the density structure along the rotational axis in the bottom panel of Fig. 7. This shows that the density decreases with time gradually, so that the heating rate required to achieve the outflow with the high Lorentz factor also decreases with time. For example, for $t\gtrsim$400 ms, the density decreases to $\lesssim 10^{5}\ {\rm g\ cm^{-3}}$. Then the heating rate required to launch the ultra-relativistic ejecta is $Q\sim 10^{31}\ {\rm erg\ cm^{-3}\ s^{-1}}$. Such a high heating rate is indeed achieved in the first $\sim 10$ ms in our simulation (see Fig. 8), at which the system has very large neutrino luminosity of $\sim 10^{53}\ {\rm erg\ s^{-1}}$. Therefore, the relativistic outflow could be launched for later times if such a high luminosity could be sustained.

4.2. Nucleosynthesis in the Neutrino-driven Ejecta

Most of the neutrino-driven ejecta found in this paper have the electron fraction between 0.3 and 0.5 and the specific entropy $\sim 10-20\ k_{\rm B}$. Therefore, the strong $r$-process may not occur in this type of the ejecta (Hoffman et al. 1997; Martin et al. 2015).

One interesting finding is that, due to the effect of the pair-annihilation heating, a small amount of very high specific entropy material of $s\sim 500\ k_{\rm B}$ is ejected (see Fig. 6). The electron fraction of such ejecta is $\sim 0.45-0.50$, and the expansion velocity is $\sim 0.2-0.5\ c$. In such high-entropy and fast-expanding material, heavy nuclei could be synthesized through the $r$-process (the condition for the $r$-process nucleosynthesis is found, e.g, see Hoffman et al. 1997). Even in a slightly proton-rich condition, if the entropy and expansion velocity of the ejecta are sufficiently high, a lot of alpha particles and nucleons remain in low-temperature environment due to alpha-rich freeze-out. Then the heavy nuclei could be produced through the nucleon capture process as described in Meyer (2002) and Fujibayashi et al. (2016). To explore these issues, we need a detailed nucleosynthesis calculation.

4.3. Effects of the Viscosity

In our simulations, the viscous effect, which is likely to be induced by the magnetohydrodynamical turbulent motion (Balbus & Hawley 1998; Hawley et al. 2013; Suzuki & Inutsuka 2014; Salvesen et al. 2016; Shi et al. 2016) is not taken into account. Recent general relativistic magnetohydrodynamics simulations for NS-NS mergers suggest that strongly magnetized MNS and torus are likely to be formed after the merger via the Kelvin-Helmholtz and the magnetorotational instabilities (Kiuchi et al. 2014, 2015a). The turbulent flow in the torus is likely to be sustained by the magnetorotational instability, and then, the flow will transport the angular momentum. Furthermore, the dissipation of the turbulent motion will heat up the torus. The dynamics of the system could be changed significantly by the effective viscosity.

First, the properties of the torus could be modified. For example, the geometrical thickness of the torus could be increased due to the viscous heating, and in addition, the torus could spread outward due to the angular momentum transport (Shibata & Kiuchi 2017). Furthermore, the viscous heating could rise the temperature inside the torus, and as a result, the neutrino luminosity of the torus would be increased. Then the heating due to neutrinos would be enhanced and will help driving a relativistic jet from the polar region of MNS as discussed in Sec. 4.1.

Second, in addition to the neutrino-driven wind, the viscosity-driven wind from the torus is expected. As suggested in the previous works (Fernández & Metzger 2013; Metzger & Fernández 2014; Fernández et al. 2015; Kiuchi et al. 2015b; Siegel & Metzger 2017), $\sim 10$ % of the torus mass may be ejected as the viscosity-driven wind, so that the properties of the mass ejection in the viscous time scale $\sim 100$ ms would be affected by the viscous effect. We plan to perform viscous radiation hydrodynamics simulations in our next work.

As the prelude, we roughly estimate the viscous effects using a snapshot of the present non-viscous simulations. In an axisymmetric, geometrically thin, and stationary accretion disk, the viscous heating rate $Q^{(+)}_{\rm vis}$ becomes

[TABLE]

where we employed a Shakura-Sunyaev parametrization of the kinetic viscosity coefficient (Shakura & Sunyaev 1973; Kato et al. 2008) as

[TABLE]

We show the viscous heating rate and neutrino emission cooling rate ($\sum_{i}Q_{({\rm leak})\nu_{i}}$) on the equatorial plane as functions of the radius in Fig. 9. Here we used the axisymmetric configuration at the beginning of the simulation. We assumed $\alpha\sim 10^{-2}$ as a fiducial value following the latest results by high-resolution magnetohydrodynamics simulations (Hawley et al. 2013; Suzuki & Inutsuka 2014; Salvesen et al. 2016; Shi et al. 2016). We find that the viscous heating rate is comparable to the cooling rate by the neutrino emission. This suggests that, due to the viscous heating, the energy loss by the neutrino cooling may be compensated, and the temperature of the torus would decrease more slowly than that in the simulation without the viscosity. Therefore, in reality, the initially large neutrino luminosity could be sustained for a longer timescale and the mass and the kinetic energy of the neutrino-driven ejecta would be larger in the presence of the viscous effect.

4.4. Uncertainties of Pair-annihilation Heating Rate

It is known that the M1-closure scheme has a problem in the optically thin region (e.g., Jiang et al. 2014; Ohsuga & Takahashi 2016; Takahashi et al. 2016; Frank et al. 2006). As a result, the pair-annihilation heating rate evaluated in our present work can be significantly different from more realistic one. This can be understood, for example, from the fact that the cross section of the pair-annihilation process has strong angle-dependence ($\propto(1-\cos\Theta)^{2}$ with $\Theta$ being the scattering angle between two neutrinos).

In Just et al. (2015), radiation fields obtained with their M1 scheme are compared to those obtained with a ray-tracing method. It is shown that the radiation energy density around the rotational axis is unphysically enhanced in the M1 scheme. Although they compared only the energy density and flux of neutrinos, their result implies that our pair-annihilation heating rate could be overestimated because of the enhancement of energy density. On the other hand, our estimation of the average energy deposited per reaction (see Eq. (39)) could be underestimated because we use the local matter temperature for its estimation.

Here, we estimate the heating rates in terms of a leakage-based volume-integration method described in Ruffert et al. (1997) using each snapshot of the simulation and compare the heating rate with that of our M1-scheme. According to Ruffert et al. (1997), the heating rate due to the pair-annihilation in this scheme can be written (in $\hbar=1$ unit) as {widetext}

[TABLE]

where $dV=d^{3}x\sqrt{\gamma}$ is the physical volume element, $\cos\Theta=(x-x_{1})\cdot(x-x_{2})/|x-x_{1}||x-x_{2}|$ is the angle between the momenta of colliding neutrinos, and $Q^{(-)}_{\rm(leak)}{}_{\nu_{i}}(x_{1})$ is the leakage source term for neutrinos of the flavor $i$ except for the heating processes. Here we do not consider general relativistic effects such as the bending of neutrino trajectories and gravitational redshift. For the average energy of neutrinos $\left<\omega\right>_{\rm pair}$, we employ the value used in the radiation hydrodynamics simulation (i.e., Eq. (39)) and focus only on the difference in the scattering angle.

Figure 10 shows the pair-annihilation heating rate using the volume-integration method along the rotational ($z$-) axis. It is found that the heating rate based on the moment formalism is by a factor of $\gtrsim 10$ smaller than the results in the volume-integration framework.

On the other hand, the pair-annihilation heating rate based on the assumption of isotropic angular distribution of neutrinos could provide a result that agrees with that in the volume-integration framework in a better manner than that in the moment formalism for the region shown in Fig. 10. This is because, near the MNS, neutrinos propagate in various directions due to the neutrino emission from the torus and the MNS. This suggests that, for the region at which the pair-annihilation heating is efficient, we need to take into account the angular distribution of neutrinos more carefully.

5. Summary

We performed a fully general relativistic, axisymmetric numerical relativity simulation for a MNS surrounded by a torus, which is the typical remnant of the NS-NS merger. We took into account neutrino transport using the truncated moment formalism with a M1-closure and relevant weak interaction reactions between neutrinos and the fluid material including the neutrino pair-annihilation in an approximate manner. For the initial condition of this simulation, we used a configuration obtained in a three-dimensional, numerical relativity simulation for the NS-NS merger. Our purpose is to investigate the amount and the properties of the material ejected due to the neutrino heating in the framework of purely radiation-hydrodynamics simulation.

In this setup, a quasi-stationary neutrino-driven outflow is launched for $\sim$300 ms from the beginning of the simulation. The effect of the neutrino pair-annihilation heating is quite large because of the very high neutrino luminosity of the MNS and the torus. Furthermore, due to the existence of the dense and hot torus, the structure of the heating rate density is quite different from the isolated NS usually considered as normal core-collapse supernova remnants. For the DD2 EOS, the mass and the kinetic energy of the neutrino-driven ejecta are comparable to those of the dynamical ejecta (Sekiguchi et al. 2015). We expect that this would be also the case for stiff EOS or for the merger of low-mass NS. Therefore, the neutrino-driven ejecta would contribute to the mass and kinetic energy of the entire ejecta of the NS-NS merger for this EOS.

The relativistic outflow required for SGRBs is not found in our present simulation because the specific heating rate around the rotational axis is not sufficiently high for driving it. The reasons for this might be the small pair-annihilation heating rate and the baryon pollution due to the neutrino-driven wind from the MNS. Moreover, the kinetic energy of the ejecta is smaller than the typical value of the SGRBs. Our results suggest that, in this purely-radiation hydrodynamics case, the neutrino pair-annihilation process in MNS-torus systems cannot account for the majority of observed SGRB events. However, we used an approximate neutrino transport scheme at this point, and hence the final conclusion should be drawn by performing simulations fully solving Boltzmann’s equation for neutrino transport. Moreover, we do not consider the effects of magnetohydrodynamic turbulence in our simulations although the MNS is likely to be strongly magnetized and the resulting turbulent motion could play an important role for its evolution (Kiuchi et al. 2015a). Hence the neutrino emissivity could be underestimated because the viscous heating associated with the turbulence motion would enhance the neutrino luminosity. For this reason, in the future, we plan to perform more sophisticated simulations considering these missing elements.

In order to obtain the relativistic outflow from this MNS-torus system, the high neutrino luminosity of $\sim 10^{53}\ {\rm erg\ s^{-1}}$ would be needed at later time, for which the density in the polar region decreases to $\lesssim 10^{5}\ {\rm g\ cm^{-3}}$. Such high neutrino luminosity may be achieved due to the viscous heating in the accretion torus.

S. F. thanks Kunihito Ioka, Yudai Suwa, and Takahiro Tanaka for fruitful discussions, Kohsuke Sumiyoshi for discussions about the radiation transfer and the emission processes of neutrinos, and the anonymous referee for helping to improve the quality of this article. This work was initiated during S.F.’s stay in Toho University, and S. F. also thanks Tetsu Kitayama for willingly accepting S.F.’s visit to the university. S. F. is supported by Research Fellowship of Japan Society for the Promotion of Science (JSPS) for Young Scientists (No.26-1329). This work was supported by JSPS Grant-in-Aid for Scientific Research (16H02183, 15H00836, 15K05077, 15H00783, 16K17706, 16H06341, 15H00782), by the HPCI Strategic Program of Japanese MEXT, and by a MEXT program “Priority Issue 9 to be tackled by using Post K computer”. Numerical simulations were carried out on Cray XC40 in Yukawa Institute for Theoretical Physics, Kyoto University and FX10 in Information Technology Center, the University of Tokyo. Numerical analyses were carried out on computers at Center for Computational Astrophysics, National Astronomical Observatory of Japan.

Appendix A A. Prescription of Neutrino-antineutrino Pair-annihilation Heating Rate in Energy-integrated M1-scheme

In the fluid rest frame, the energy deposition rate for a neutrino species $\nu_{i}$ due to the pair-annihilation process is written (in $\hbar=1$ unit) as

[TABLE]

where $f_{\nu_{i}}(k)$ and $f_{\bar{\nu}_{i}}(k^{\prime})$ are the distribution functions of the $i$-th species of neutrinos and antineutrinos, which are functions of momenta $k$ and $k^{\prime}$, respectively, and $\omega=-u^{\alpha}k_{\alpha}$ is the energy of neutrinos in the fluid rest frame. The cross section of the pair-annihilation process ($\nu_{i}+\bar{\nu}_{i}\rightarrow e^{-}+e^{+}$), $\sigma_{{\rm pair},i}(k,k^{\prime})$, is written (see, e.g., Salmonson & Wilson 1999) as

[TABLE]

where we ignored the phase space blocking for electrons and positrons, and the electron mass ($\approx 0.511$ MeV) because the energy of neutrinos considered in this paper ($\agt 5$ MeV) is much higher. Here we defined a spatial unit four-vector $\ell^{\alpha}$ orthogonal to $u^{\alpha}$ as

[TABLE]

In our moment formalism (Shibata et al. 2011), we define the energy-integrated zeroth-, first-, and second-rank moments of neutrinos by

[TABLE]

where $\Omega$ denotes the solid angle in the momentum space of neutrinos in the fluid rest frame. In order to evaluate the heating rate using the above moments, we approximate the energy of neutrinos, $\omega$, by its average value $\left<\omega\right>_{\rm pair}$. Then Eq. (A1) becomes

[TABLE]

which appeared in Eq. (26).

Appendix B B. Prescription of Neutrino Absorption Processes in Energy-integrated M1-scheme

In the same way as Eq. (A1), the source terms due to the absorption of electron (anti)neutrinos are written (in $\hbar=1$ unit) as

[TABLE]

The cross section of the absorption processes is written as

[TABLE]

where we ignored the electron mass and the difference between proton and neutron masses ($m_{n}-m_{p}\approx 1.293$ MeV). As in Eq. (A7), the energy deposition rates are evaluated in our scheme as

[TABLE]

The average energy estimated by Eq. (40) is the simple spectral average assuming the Fermi-Dirac-type energy distribution for neutrinos. On the other hand, the average energy estimated by Eq. (41) is derived when we take the energy dependence of the absorption reaction into account. We can directly perform the energy integral in Eq. (B1) as

[TABLE]

where we used Eq. (37). Thus, the average energy of this reaction is the square root of the square bracket, i.e., Eq. (41).

Using these rates, the source terms for electrons (Eq.(11)) are written as

[TABLE]