Mass Ejection from the Remnant of a Binary Neutron Star Merger: Viscous-radiation Hydrodynamics Study

In our code, Einstein's equation is solved using a version of the Baumgarte–Shapiro–Shibata–Nakamura puncture formalism (Shibata & Nakamura 1995; Baumgarte & Shapiro 1999; Baker et al. 2006; Campanelli et al. 2006). The quantities we evolve in our formalism (in the Cartesian coordinate basis) are listed in Table 1. Here, from the spacetime metric g_αβ and the time-like unit vector normal to spatial hypersurfaces n_α, we define the induced spatial metric as γ_αβ ≡ g_αβ + n_αn_β. In addition, we define the determinant of the induced metric as γ = det(γ_ij) and the extrinsic curvature as ${K}_{\alpha \beta }\,=-(1/2){{ \mathcal L }}_{n}{\gamma }_{\alpha \beta }$, where ${{ \mathcal L }}_{n}$ is the Lie derivative with respect to n^α. For the gauge conditions, we employ dynamical lapse and shift gauge conditions that are the same as Equations (1) and (2) in Fujibayashi et al. (2017). We adopt the so-called cartoon method (Alcubierre et al. 2001; Shibata 2003) to impose axisymmetric conditions for the geometric quantities.

Table 1.  List of the Quantities that We Evolve in Our Baumgarte–Shapiro–Shibata–Nakamura Puncture Formulation

Notation	Definition
${\tilde{\gamma }}_{{ij}}={\gamma }^{-1/3}{\gamma }_{{ij}}$	Conformal three-metric
W = γ−1/6	Conformal factor
K = γijKij	Trace of the extrinsic curvature Kij
${\tilde{A}}_{{ij}}={\gamma }^{-1/3}({K}_{{ij}}-{\gamma }_{{ij}}K/3)$	Trace-free part of Kij
${F}_{i}={\delta }^{{jk}}{\partial }_{j}{\tilde{\gamma }}_{{ki}}$	Auxiliary variable

Download table as:  ASCIITypeset image

As in Fujibayashi et al. (2017), we decompose neutrinos into "streaming" and "trapped" components and write the total energy-momentum tensor of the matter (fluid and neutrinos) as

$$\begin{eqnarray}&&{T}_{(\mathrm{tot})}^{\alpha \beta }={T}^{\alpha \beta }+\displaystyle \sum _{i}{T}_{({\nu }_{i},{\rm{S}})}^{\alpha \beta },\end{eqnarray} \tag{ 1 }$$

where ${T}^{\alpha \beta }={T}_{(\mathrm{fluid})}^{\alpha \beta }+{\sum }_{i}{T}_{({\nu }_{i},T)}^{\alpha \beta }$ is the energy-momentum tensor composed of the sum of the fluid and trapped neutrinos, and ${T}_{({\nu }_{i},S)}^{\alpha \beta }$ denotes the energy-momentum tensor for free-streaming neutrinos. These energy-momentum tensors obey the following equations:

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{\beta }{T}^{\alpha \beta }=-{Q}_{(\mathrm{leak})}^{\alpha }=-\displaystyle \sum _{i}{Q}_{(\mathrm{leak}){\nu }_{i}}^{\alpha },\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{\beta }{T}_{({\nu }_{i},{\rm{S}})}^{\alpha \beta }={Q}_{(\mathrm{leak}){\nu }_{i}}^{\alpha },\end{eqnarray} \tag{ 3 }$$

where ${Q}_{(\mathrm{leak}){\nu }_{i}}^{\alpha }$ denotes the leakage rate of ith species of neutrinos. For solving the evolution equations for free-streaming neutrinos, we employ the so-called M1 scheme with a closure relation (Shibata et al. 2011), and, for trapped neutrinos, we employ a leakage-based scheme developed by Sekiguchi (2010). The detailed description of these schemes is found in Sekiguchi (2010) and Fujibayashi et al. (2017).

Then, we describe the basic equations for viscous hydrodynamics. The formulation of our general relativistic viscous hydrodynamics is described in Shibata et al. (2017b) and Shibata & Kiuchi (2017). In the following, we briefly outline our formulation.

The energy-momentum tensor of the viscous fluid with trapped neutrinos is given as

$$\begin{eqnarray}&&{T}_{\alpha \beta }=\rho {{hu}}_{\alpha }{u}_{\beta }+{{Pg}}_{\alpha \beta }-\rho h\nu {\tau }^{0}{}_{\alpha \beta },\end{eqnarray} \tag{ 4 }$$

where ρ is the baryon rest-mass density, h = 1 +  + P/ρ is the specific enthalpy, is the specific internal energy, u^α is the fluid four-velocity, P is the pressure, ν is the viscosity coefficient, and τ⁰_αβ is the viscous stress tensor. Here τ⁰_αβ is a symmetric tensor that satisfies the relation τ⁰_αβu^α = 0, and, in the formalism of Israel & Stewart (1979), it obeys the evolution equation

$$\begin{eqnarray}&&{{ \mathcal L }}_{u}{\tau }^{0}{}_{\alpha \beta }=-\zeta ({\tau }^{0}{}_{\alpha \beta }-{\sigma }_{\alpha \beta }),\end{eqnarray} \tag{ 5 }$$

where σ_αβ is the shear tenor defined by

$$\begin{eqnarray}&&{\sigma }_{\alpha \beta }={h}_{\alpha }{}^{\mu }{h}_{\beta }{}^{\nu }({{\rm{\nabla }}}_{\mu }{u}_{\nu }+{{\rm{\nabla }}}_{\nu }{u}_{\mu })={{ \mathcal L }}_{u}{h}_{\alpha \beta },\end{eqnarray} \tag{ 6 }$$

and h_αβ = g_αβ + u_αu_β. Here ζ is a nonzero constant of (time)⁻¹ dimension, and it has to be chosen in an appropriate manner for τ⁰_αβ to approach σ_αβ in a short timescale because it is reasonable to suppose that ${\tau }_{\alpha \beta }^{0}$ should approach σ_αβ in a microscopic timescale.

We can rewrite Equation (5) as

$$\begin{eqnarray}&&{{ \mathcal L }}_{u}{\tau }_{\alpha \beta }=-\zeta {\tau }^{0}{}_{\alpha \beta },\end{eqnarray} \tag{ 7 }$$

where ${\tau }_{\alpha \beta }\equiv {\tau }^{0}{}_{\alpha \beta }-\zeta {h}_{\alpha \beta }$. Thus, in addition to hydrodynamics equations, we solve Equation (7) as a basic equation that describes the evolution of τ_αβ.

The energy-momentum tensor of the fluid is rewritten as

$$\begin{eqnarray}{T}_{\alpha \beta } & = & \rho h(1-\nu \zeta ){u}_{\alpha }{u}_{\beta }+(P-\rho h\nu \zeta ){g}_{\alpha \beta }-\rho h\nu {\tau }_{\alpha \beta }\\ & = & \rho w\hat{e}{n}_{\alpha }{n}_{\beta }+{J}_{\alpha }{n}_{\beta }+{J}_{\beta }{n}_{\alpha }+{S}_{\alpha \beta },\end{eqnarray} \tag{ 8 }$$

where we defined

$$\begin{eqnarray} & & \hat{e}\equiv {T}_{\alpha \beta }{n}^{\alpha }{n}^{\beta }/\rho {w}\\ & = & {hw}(1-\nu \zeta )-\displaystyle \frac{P-\rho h\nu \zeta }{\rho w}-h\nu {w}^{-3}{\tau }_{{ij}}{V}^{i}{V}^{j},\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}{J}_{i} & \equiv & -{T}_{\alpha \beta }{n}^{\alpha }{\gamma }_{i}{}^{\beta }\\ & = & \rho {{hwu}}_{i}(1-\nu \zeta )-\rho h\nu {w}^{-1}{\tau }_{{ij}}{V}^{j},\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}{S}_{{ij}} & \equiv & {T}_{\alpha \beta }{\gamma }^{\alpha }{}_{i}{\gamma }^{\beta }{}_{j}\\ & = & \rho h(1-\nu \zeta ){u}_{i}{u}_{j}+(P-\rho h\nu \zeta ){\gamma }_{{ij}}-\rho h\nu {\tau }_{{ij}}.\end{eqnarray} \tag{ 11 }$$

Here w ≡ −n_αu^α (the Lorentz factor of the fluid), and ${V}^{i}\,\equiv {\gamma }^{{ij}}{u}_{j}$. The general relativistic Navier–Stokes and energy equations are derived, respectively, from the space-like and time-like components of Equation (2), i.e., ${\gamma }_{i\alpha }{{\rm{\nabla }}}_{\beta }{T}^{\alpha \beta }\,=-{\gamma }_{i\alpha }{Q}_{(\mathrm{leak})}^{\alpha }$ and ${n}_{\alpha }{{\rm{\nabla }}}_{\beta }{T}^{\alpha \beta }=-{n}_{\alpha }{Q}_{(\mathrm{leak})}^{\alpha }$, as

$$\begin{eqnarray}&&{\partial }_{t}(\sqrt{\gamma }{J}_{i})+{\partial }_{k}[\sqrt{\gamma }(\alpha {S}^{k}{}_{i}-{\beta }^{k}{J}_{i})]\\ &&\quad =\sqrt{\gamma }\left(-{\rho }_{* }\hat{e}{\partial }_{i}\alpha +{J}_{k}{\partial }_{i}{\beta }^{k}-\displaystyle \frac{\alpha }{2}{S}_{{kl}}{\partial }_{i}{\gamma }^{{kl}}-{\gamma }_{i\alpha }{Q}_{(\mathrm{leak})}^{\alpha }\right),\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{\partial }_{t}({\rho }_{* }\hat{e})+{\partial }_{k}[\sqrt{\gamma }(\alpha {J}^{k}-{\beta }^{k}{\rho }_{* }\hat{e})]\\ &&\quad =\sqrt{\gamma }({S}_{{kl}}{K}^{{kl}}-{\gamma }^{{kl}}{J}_{k}{\partial }_{l}\alpha +{n}_{\alpha }{Q}_{(\mathrm{leak})}^{\alpha }),\end{eqnarray} \tag{ 13 }$$

where ${\rho }_{* }\equiv \rho w\sqrt{\gamma }$, and α and β^k denote the lapse function and shift vector, respectively. On the other hand, the evolution equation of τ_ij is derived from Equation (7) as

$$\begin{eqnarray}&&{\partial }_{t}({\rho }_{* }{\tau }_{{ij}})+{\partial }_{k}({\rho }_{* }{\tau }_{{ij}}{v}^{k})\\ &&\quad =-{\rho }_{* }({\tau }_{{ik}}{\partial }_{j}{v}^{k}+{\tau }_{{jk}}{\partial }_{i}{v}^{k})-{\rho }_{* }\alpha {w}^{-1}\zeta {\tau }^{0}{}_{{ij}},\end{eqnarray} \tag{ 14 }$$

where we used the equation of continuity ${\partial }_{t}{\rho }_{* }+{\partial }_{k}({\rho }_{* }{v}^{k})=0$ and v^k ≡ u^k/u^t. We solve these equations in cylindrical coordinates as in Shibata et al. (2017b).

For the neutron star EOS, we employ a tabulated EOS referred to as DD2 (Banik et al. 2014). We extend this EOS to low-temperature and low-density ranges down to 10⁻³ MeV and ≈0.1 g cm⁻³ using an EOS by Timmes & Swesty (2000) as in our previous work (Fujibayashi et al. 2017). This extended DD2 EOS includes contributions of nucleons, heavy nuclei, electrons, positrons, and photons to the pressure and internal energy.

In our formalism (Sekiguchi 2010; Fujibayashi et al. 2017), we solve the equation for the energy-momentum density of streaming neutrinos and number density of electrons and trapped neutrinos. For the source terms of these, we take into account the relevant processes due to the weak interaction. For (anti)neutrino production processes, we adopt the electron and positron capture of nucleons and heavy nuclei, electron–positron pair annihilation, nucleon–nucleon bremsstrahlung, and plasmon decay based on Fuller et al. (1985), Cooperstein et al. (1986), Burrows et al. (2006), and Ruffert et al. (1996), respectively. On the other hand, for (anti)neutrino absorption processes, we adopt the absorption of neutrinos and antineutrinos on nucleons and heavy nuclei and neutrino–antineutrino pair annihilation to electron–positron pairs. The detailed description of calculating the rates is also found in Fujibayashi et al. (2017).

Since we use the energy-integrated formalism for neutrino radiation transfer, we have to assume the average energy of neutrinos for calculating the rates of neutrino absorption and neutrino pair-annihilation processes, which depend strongly on the energy of neutrinos. In Fujibayashi et al. (2017), we found that the results of the simulation, such as the ejecta mass and kinetic energy, depend weakly on the assumption for the average neutrino energy. In this work, we adopt the average energy estimated by Equation (41) in Fujibayashi et al. (2017) because it can be derived directly from the energy integration of the energy-dependent neutrino absorption heating rate.

In this work, we model the dynamical shear viscosity coefficient following Shakura & Sunyaev (1973) as

$$\begin{eqnarray}&&\nu ={\alpha }_{\mathrm{vis}}{c}_{s}{H}_{\mathrm{tur}},\end{eqnarray} \tag{ 15 }$$

where α_vis is the so-called α-viscosity parameter, c_s is the sound speed, and H_tur is the possible largest scale of eddies in a turbulence state. In this work, we set ${H}_{\mathrm{tur}}=10$ km, because it should be approximately equal to the size of the MNS (for the central region of the system).

As already mentioned in the previous subsection, ζ⁻¹ has to be a timescale short enough that ${\tau }_{\alpha \beta }^{0}$ quickly approaches σ_αβ. In this work, we set ζ = 3 × 10⁴ s⁻¹. This value is about 4 times larger than the maximum angular velocity of the system (${\rm{\Omega }}\sim 8000\,\mathrm{rad}\,{{\rm{s}}}^{-1}$); thus, the requirement for ζ is reasonably satisfied.

We note that for the outer part of the torus, the viscosity coefficient would be underestimated because the scale height of the torus increases with the radius in general. Thus, in our prescription, we conservatively take into account the effect of the viscosity for large radii.

The initial condition in this work is given by the same method as that in Fujibayashi et al. (2017): we first perform a three-dimensional, numerical-relativity simulation for the merger of equal-mass binary neutron stars with the total mass 2.7 M_⊙ (Sekiguchi et al. 2015) using DD2 EOS (Banik et al. 2014). For these choices of the total mass and EOS, the remnant formed after the merger is a long-lived MNS surrounded by a torus. At a few tens of ms after the formation of this system, the MNS and torus relax approximately to an axisymmetric configuration; hence, we employ such a system as the initial condition. The total baryon and Arnowitt--Deser--Misner (ADM) mass of the initial condition are $\approx 2.95\,{M}_{\odot }$ and $2.65\,{M}_{\odot }$, respectively. We note that $\approx 0.05\,{M}_{\odot }$ is carried by gravitational radiation during the inspiral and merger. The baryon mass of the torus of the initial condition is $\approx 0.2\,{M}_{\odot }$. Here we define the torus mass by

$$\begin{eqnarray}&&{M}_{{\rm{b}},\mathrm{torus}}={\int }_{\rho \lt {10}^{13}{\rm{g}}{\mathrm{cm}}^{-3}}{d}^{3}x\ {\rho }_{* }.\end{eqnarray} \tag{ 16 }$$

We note that the mass of the merger-remnant torus also depends on the total mass, mass ratio of the binary, and neutron star EOS. Therefore, the result in this work would be quantitatively different for different parameters and EOSs.

In addition to the dynamical variables employed in Fujibayashi et al. (2017), we need to prepare an initial condition for the spatial components of the viscous tensor τ_ij. For simplicity, we assume that the viscous tensor ${\tau }^{0}{}_{\alpha \beta }$ vanishes in the initial state. That is, as the initial condition for τ_ij, we set τ_ij = −ζ h_ij.

For evolving radiation-viscous hydrodynamics equations in cylindrical coordinates (x, z), we employ the same nonuniform grid as that in Fujibayashi et al. (2017), in which the grid structure is determined by the uniform grid spacing of the inner region, Δx₀; the range of the inner region, R_star; the increase rate of the grid spacing in the outer region, 1 + δ; and the grid number, N. For all the models in this paper, we set R_star = 30 km and δ = 0.0075. We list Δx₀ and N, together with the size of the computational domain, L, in Table 2. We performed simulations with three different viscosity parameters, α_vis =  0.01, 0.02, and 0.04, among which we refer to the model with α_vis =  0.02 as our fiducial model. Our choice of the viscosity parameter α_vis = O(0.01) is justified, at least for the outer part of the MNS (ρ ≲ 10¹⁴ g cm⁻³), by the recent high-resolution global MHD simulations for a binary neutron star merger (Kiuchi et al. 2017). For a higher-density region, such as an inner region of the MNS, the value of the viscosity parameter is still uncertain, but for simplicity, we use the same value of the viscosity parameter for the entire region, supposing that a turbulence state would be achieved entirely in the MNS.⁴

To investigate the long-term evolution of the system, the simulations are performed for ∼3 s for our fiducial and α_vis = 0.04 models, while the simulations of the other models are performed for a shorter time (≲1 s).

In addition to these models, to confirm that numerical results with different grid resolutions reasonably agree, we performed simulations for α_vis = 0.01 with two lower resolutions as Δx₀ = 200 m (M) and 250 m (L). The grid structure for these cases is also listed in Table 2. The result of this check is described in the Appendix.

From the component of Equation (2) along u^α, we obtain

$$\begin{eqnarray}-{Q}_{(\mathrm{leak})}^{\alpha }{u}_{\alpha } & = & -\rho {u}^{\alpha }{{\rm{\nabla }}}_{\alpha }h+{u}^{\alpha }{{\rm{\nabla }}}_{\alpha }P-{u}^{\alpha }{{\rm{\nabla }}}^{\beta }(\rho h\nu {\tau }^{0}{}_{\alpha \beta })\\ & = & -\rho {{Tu}}^{\alpha }{{\rm{\nabla }}}_{\alpha }s+\rho h\nu {\tau }^{0}{}_{\alpha \beta }{{\rm{\nabla }}}^{\alpha }{u}^{\beta },\end{eqnarray} \tag{ 17 }$$

where we used the relation ${\tau }^{0}{}_{\alpha \beta }{u}^{\alpha }=0$ and the first law of thermodynamics dh = dP/ρ + Tds, with T and s being the temperature and specific entropy. From the above relation, the viscous heating rate in the fluid rest frame in our formalism can be written approximately as

$$\begin{eqnarray}&&{Q}_{\mathrm{vis}}^{(+)}=\rho h\nu {\tau }^{0}{}^{\alpha \beta }{{\rm{\nabla }}}_{\alpha }{u}_{\beta }=\displaystyle \frac{1}{2}\rho h\nu {\tau }^{0}{}^{\alpha \beta }{\tau }^{0}{}_{\alpha \beta },\end{eqnarray} \tag{ 18 }$$

where we assumed that ${\tau }^{0}{}_{\alpha \beta }$ quickly approaches σ_αβ.

Using the y component of Equation (12), the angular momentum transport flux due to the viscous effect is described by

$$\begin{eqnarray}&&-\rho h\nu {\bar{\tau }}^{0}{}^{i}{}_{y}=-\rho h\nu {\gamma }^{{ik}}{\tau }^{0}{}_{{ky}}.\end{eqnarray} \tag{ 19 }$$

From the surface integral of this, the local angular momentum transport rate is derived.

Before describing the numerical results, we estimate the timescales for the evolution of the MNS-torus system in the presence of viscosity. In this problem, the viscous effect plays two crucial roles.

In the first ≲20 ms, the viscous effect plays an important role in the MNS. The MNS formed after the merger initially has a differentially rotating velocity profile. By the viscous effect, the angular momentum in a high angular velocity region is transported to a low angular velocity region; thus, the MNS approaches a rigidly rotating state. The timescale for the angular momentum transport in the MNS is estimated by

$$\begin{eqnarray}\displaystyle \frac{{R}_{\mathrm{eq}}{}^{2}}{\nu } & \approx & \,10\,\mathrm{ms}{\left(\displaystyle \frac{{\alpha }_{\mathrm{vis}}}{0.02}\right)}^{-1}{\left(\displaystyle \frac{{c}_{s}}{c/3}\right)}^{-1}\\ & & \times {\left(\displaystyle \frac{{R}_{\mathrm{eq}}}{15\mathrm{km}}\right)}^{2}{\left(\displaystyle \frac{{H}_{\mathrm{tur}}}{10\mathrm{km}}\right)}^{-1},\end{eqnarray} \tag{ 20 }$$

where R_eq denotes the equatorial radius of the MNS. Thus, the differential rotation profile of the MNS is expected to be erased within ∼10 (α_vis/0.02)⁻¹ ms by the viscous effect.

For longer timescales, the viscous effect plays a primary role in determining the evolution of the torus surrounding the MNS. Assuming that the torus evolution could be described by a standard accretion disk theory (Shakura & Sunyaev 1973), the viscous accretion timescale is estimated by

$$\begin{eqnarray}{t}_{\mathrm{vis}} & \sim & \ {\alpha }_{\mathrm{vis}}{}^{-1}{\left(\displaystyle \frac{{R}_{\mathrm{torus}}{}^{3}}{{G{M}}_{\mathrm{MNS}}}\right)}^{1/2}{\left(\displaystyle \frac{H}{{R}_{\mathrm{torus}}}\right)}^{-2}\\ & \approx & \ 0.55\,{\rm{s}}\ {\left(\displaystyle \frac{{\alpha }_{\mathrm{vis}}}{0.01}\right)}^{-1}{\left(\displaystyle \frac{{M}_{\mathrm{MNS}}}{2.5{M}_{\odot }}\right)}^{-1/2}\\ & & \times {\left(\displaystyle \frac{{R}_{\mathrm{torus}}}{50\mathrm{km}}\right)}^{3/2}{\left(\displaystyle \frac{H/{R}_{\mathrm{torus}}}{1/3}\right)}^{-2},\end{eqnarray} \tag{ 21 }$$

where M_MNS is the gravitational mass of the MNS, ${R}_{\mathrm{torus}}$ is the typical radius of the torus, and H/R_torus is the aspect ratio of the torus. Then, the mass accretion rate is estimated as

$$\begin{eqnarray}{\dot{M}}_{\mathrm{MNS}} & \sim & \displaystyle \frac{{M}_{\mathrm{torus}}}{{t}_{\mathrm{vis}}}\\ & \approx & \ 0.36\ {M}_{\odot }\,{{\rm{s}}}^{-1}\ \left(\displaystyle \frac{{\alpha }_{\mathrm{vis}}}{0.01}\right){\left(\displaystyle \frac{{M}_{\mathrm{MNS}}}{2.5{M}_{\odot }}\right)}^{1/2}\\ & & \times {\left(\displaystyle \frac{{R}_{\mathrm{torus}}}{50\mathrm{km}}\right)}^{-3/2}{\left(\displaystyle \frac{H/{R}_{\mathrm{torus}}}{1/3}\right)}^{2}\left(\displaystyle \frac{{M}_{\mathrm{torus}}}{0.2\ {M}_{\odot }}\right).\end{eqnarray} \tag{ 22 }$$

We note that in these estimates, we used the standard accretion disk in a stationary state as a model. In reality, the MNS-torus system is evolved by the viscous timescale (i.e., the values of R_torus and M_torus change with time); hence, we should check whether the stationary disk model is valid or not by performing the simulation of this system.

First, we briefly describe the evolution process of the system for the fiducial model DD2-135135-0.02-H.

In the early stage of the evolution, the angular momentum transport occurs in the MNS. Figure 1 shows the angular velocity profiles in the equatorial plane at different time slices for the fiducial model. This shows that initial differential rotation of the MNS within the radius of ∼12 km disappears and the MNS settles into a rigid rotation state in ∼10 ms, as estimated by Equation (20).

Zoom In Zoom Out Reset image size

Then, a sound wave is formed in the vicinity of the MNS. This is due to the variation of the MNS density profile caused by the angular momentum redistribution.

The variation of the rotational kinetic energy in the redistribution of the rotational profile is approximately estimated by

$$\begin{eqnarray}\displaystyle \frac{1}{2}{\rm{\Delta }}({v}^{2})M & \sim & \displaystyle \frac{1}{2}[{\rm{\Delta }}({{\rm{\Omega }}}^{2}){R}^{2}]M\\ & \approx & 2.3\times {10}^{52}\,\mathrm{erg}\ \left(\displaystyle \frac{{\rm{\Delta }}({{\rm{\Omega }}}^{2})}{9\times {10}^{6}\,{\mathrm{rad}}^{2}\,{{\rm{s}}}^{-2}}\right)\\ & & \times {\left(\displaystyle \frac{R}{10\mathrm{km}}\right)}^{2}\left(\displaystyle \frac{M}{2.5{M}_{\odot }}\right),\end{eqnarray} \tag{ 23 }$$

where Δ(v²) and Δ(Ω²) are the variations of the rotational velocity squared and angular velocity squared, and M is the mass suffered from this variation (approximately equal to ${M}_{\mathrm{MNS}}$). This shows that the energy of ∼10⁵² erg could be redistributed in the inner region of the MNS in its viscous timescale, which is proportional to ${\alpha }_{\mathrm{vis}}{}^{-1}$. This implies that the sound wave becomes stronger for the models with higher viscosity parameters.

During its outward propagation, the sound wave becomes a shock wave (see panel (1) of Figure 2). It sweeps the material surrounding the torus and then induces mass ejection of the material for the first ∼50 ms (see panels (1) and (2) of Figure 2). We refer to this mass ejection as "early viscosity-driven mass ejection."

Zoom In Zoom Out Reset image size

The material swept up by the shock wave for small radii (r ≲ 500 km) is still gravitationally bound while it continues to expand for ∼0.1 s, as shown in panel (3) of Figure 2. After ∼0.1 s, the material begins to turn over and fall again (see panel (4) of Figure 2). We note that for the model with α_vis = 0.04, this turnover occurs only weakly and more mass is ejected from the system (see the discussion of Section 3.4.2).

After the early viscosity-driven mass ejection, the neutrino- and viscosity-driven mass ejection is developed in the polar region. It is difficult to separate the contributions of the heating due to the neutrino irradiation and viscosity, but for t ≳ 0.2 s, the viscous heating becomes important rather than the neutrino heating. This is because the neutrino pair-annihilation heating, which is the primary neutrino heating source, plays a minor role due to the decrease of the neutrino luminosity in this phase for the inviscid model (Fujibayashi et al. 2017).

The material in the torus moves in various directions for the first tens of ms (see panel (1) of Figure 2), but the flow structure gradually relaxes to a laminal state due to the viscosity-driven redistribution (see panels (3) and (4) of Figure 2). After the relaxation, the MNS-torus system evolves in a quasi-stationary manner. The torus material around the equatorial plane expands outward, while the torus material near the torus surface accretes onto the central MNS (see also Figure 5 in Section 3.2). Figure 3 shows the density profiles on the equatorial plane at different time slices. As found in this figure, the torus gradually expands with time; hence, the torus density decreases. This behavior of the torus is determined by the viscous angular momentum transport (see Section 3.2 for details).

Zoom In Zoom Out Reset image size

The top and middle panels of Figure 4 show the time evolution of the baryon mass and the angular momentum of the torus defined by

$$\begin{eqnarray}&&{J}_{\mathrm{torus}}={\int }_{\rho \lt {10}^{13}{\rm{g}}{\mathrm{cm}}^{-3}}{d}^{3}x\ {\rho }_{* }{{hu}}_{y}x.\end{eqnarray} \tag{ 24 }$$

For all of the models, they decrease with time due to mass accretion onto the central MNS and outflow.⁵

Zoom In Zoom Out Reset image size

In the phase of the early viscosity-driven mass ejection, in which the density profiles of the MNS and torus vary in a short timescale, the decrease timescale is slightly shorter than that estimated by Equation (21). However, after the torus relaxes to a quasi-stationary state, the timescale agrees approximately with that by Equation (21).

The decrease rates of M_b,torus and J_torus approach zero for t ≳ 1 s. This implies that the accretion onto the MNS becomes inefficient. This point will be revisited in Section 3.3.

Because the angular momentum profile of the torus is approximately described by the Keplerian profile around the MNS, we define a typical torus radius by

$$\begin{eqnarray}{R}_{\mathrm{torus}} & \equiv & \displaystyle \frac{{J}_{\mathrm{torus}}^{2}}{{G{M}}_{{\rm{b}},\mathrm{torus}}^{2}}\displaystyle \frac{1}{{M}_{\mathrm{MNS}}}\\ & \approx & 70\,\mathrm{km}{\left(\displaystyle \frac{{J}_{\mathrm{torus}}}{{10}^{49}\mathrm{erg}{\rm{s}}}\right)}^{2}{\left(\displaystyle \frac{{M}_{{\rm{b}},\mathrm{torus}}}{0.1{M}_{\odot }}\right)}^{-2}{\left(\displaystyle \frac{{M}_{\mathrm{MNS}}}{2.6{M}_{\odot }}\right)}^{-1}.\end{eqnarray} \tag{ 25 }$$

Here we assumed that J_torus would be approximately given by ${M}_{{\rm{b}},\mathrm{torus}}\sqrt{{G{M}}_{\mathrm{MNS}}{R}_{\mathrm{torus}}}$. The bottom panel of Figure 4 shows the radii for all the models. Here we used ${M}_{\mathrm{MNS}}=2.6\ \,{M}_{\odot }$ for all of the models, which is approximately equal to the ADM mass of the system. The typical radius is 30–40 km at the beginning of the simulation. For the viscous models, it increases with time, while the radius is approximately constant in time for the inviscid model. This shows that, by viscous angular momentum transport, the torus expands along the equatorial direction. This effect eventually induces viscosity-driven mass ejection from the torus. We will describe this process in Section 3.4.

When the central remnant composed of a rigidly rotating MNS and torus relaxes to a quasi-stationary state, the flow structure also reaches a quasi-stationary state approximately. This state is preserved until the matter temperature of the torus decreases sufficiently. The top right panel of Figure 5 shows the radial mass flux in the meridian plane at t = 0.2 s for the fiducial model. Broadly speaking, the structure of the flow is divided into three parts: the outflow near the polar region (Region I; red), the expanding flow along the equatorial plane (Region II; red), and inflow between the two regions (Region III; blue).

Zoom In Zoom Out Reset image size

The top left panel of Figure 5 displays the profile of the source of the angular momentum transport due to the viscosity (i.e., ${x}^{-2}{\partial }_{j}({x}^{2}\rho h\nu {\bar{\tau }}^{0j}{}_{y})$). This panel also shows the direction of the angular momentum flux due to the viscosity (i.e., $-\rho h\nu ({\bar{\tau }}^{x}{}_{y},{\bar{\tau }}^{z}{}_{y})$). This clearly illustrates that the angular momentum is transported from Region III to Region II. Therefore, the expanding part in Region II is driven by the angular momentum transport from the inflow in Region III.

The top left panel of Figure 6 shows the vertically and azimuthally integrated heating rate (solid red curve)

$$\begin{eqnarray}&&{ \mathcal H }=2{\int }_{0}^{L}{dz}\ 2\pi x\ ({Q}_{\mathrm{vis}}^{(+)}+{Q}_{\nu }^{(+)})\end{eqnarray} \tag{ 26 }$$

and cooling rate (solid blue curve)

$$\begin{eqnarray}&&{ \mathcal C }=2{\int }_{0}^{L}{dz}\ 2\pi x\ {Q}_{(\mathrm{leak})}^{(-)}.\end{eqnarray} \tag{ 27 }$$

Here ${Q}_{\nu }^{(+)}$ is the sum of the matter-heating source terms due to the neutrino absorption and pair-annihilation processes (see Fujibayashi et al. 2017 for the description of the individual heating rates) and ${Q}_{(\mathrm{leak})}^{(-)}$ is the sum of the leakage source terms for all flavors of neutrinos. This figure shows that the viscous heating rate balances approximately with the net neutrino cooling rate.

Zoom In Zoom Out Reset image size

In the same panel, the heating and cooling rates in the outgoing and ingoing regions,

$$\begin{eqnarray}&&{{ \mathcal H }}_{\left\{\begin{array}{c}\lt \\ \gt \end{array}\right\}}=2{\int }_{0,\left\{\begin{array}{c}{v}^{x}\lt 0\\ {v}^{x}\gt 0\end{array}\right\}}^{L}{dz}\ 2\pi x\ ({Q}_{\mathrm{vis}}^{(+)}+{Q}_{\nu }^{(+)}),\end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}&&{{ \mathcal C }}_{\left\{\begin{array}{c}\lt \\ \gt \end{array}\right\}}=2{\int }_{0,\left\{\begin{array}{c}{v}^{x}\lt 0\\ {v}^{x}\gt 0\end{array}\right\}}^{L}{dz}\ 2\pi x\ {Q}_{(\mathrm{leak})}^{(-)},\end{eqnarray} \tag{ 29 }$$

are plotted in the dashed and dotted curves, respectively. These clearly show that the balances between the viscous heating and the net neutrino cooling are approximately achieved both in Region II (expanding region) and Region III (inflowing region). This indicates that the neutrino-dominated accretion flow (NDAF) is achieved in Region III.

The middle left panel of Figure 6 shows the mass inflow and outflow rates (${\dot{M}}_{\mathrm{in}}$ and ${\dot{M}}_{\mathrm{out}}$) along the equatorial plane, together with the net rates ($\dot{M}$), defined, respectively, by

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{in}}=-2{\int }_{0,{v}^{x}\lt 0}^{L}{dz}\ 2\pi x{\rho }_{* }{v}^{x},\end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{out}}=2{\int }_{0,{v}^{x}\gt 0}^{L}{dz}\ 2\pi x{\rho }_{* }{v}^{x},\end{eqnarray} \tag{ 31 }$$

and

$$\begin{eqnarray}&&\dot{M}={\dot{M}}_{\mathrm{in}}-{\dot{M}}_{\mathrm{out}}.\end{eqnarray} \tag{ 32 }$$

In the same panel, we plot the mass accretion rate in the "standard" disk picture, i.e., a stationary thin accretion disk,

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{sd}}\approx 3\pi \nu {\rm{\Sigma }}.\end{eqnarray} \tag{ 33 }$$

Here we neglected the terms from the inner boundary condition and defined the column density of the torus by ${\rm{\Sigma }}\,=2\rho {c}_{s}/{\rm{\Omega }}{| }_{z=0}$, which is evaluated using the values on the equatorial plane. The net mass accretion rate $\dot{M}$ agrees approximately with the accretion rate in the standard disk picture ${\dot{M}}_{\mathrm{sd}}$ for the innermost region of the torus x ≲ 50 km.

In the bottom left panel of Figure 6, we compare the scale height of the torus assuming the vertically hydrostatic structure ${H}_{\mathrm{sd}}={c}_{s}/{\rm{\Omega }}{| }_{z=0}$ and the scale height calculated from the simulation data H_num to check the validity of the approximation of the flow as the standard disk. Here H_num is defined as

$$\begin{eqnarray}&&\rho (z={H}_{\mathrm{num}})=\rho (z=0)/e,\end{eqnarray} \tag{ 34 }$$

where e is the base of natural logarithms. In the innermost region (x ≲ 100 km) of the torus, the two scale heights agree well.

The flow structure, heating and cooling rates, mass accretion rate, and scale height at the late time t = 1.6 s are shown in the bottom panels of Figure 5 and the right panels of Figure 6, respectively. We find in Figure 5 that the flow structure at t = 1.6 s is qualitatively the same as that at t = 0.2 s, while the mass flux becomes smaller than that in the early phase (see the middle panel of Figure 6). At t = 1.6 s, the region for which the net mass accretion rate agrees with the standard picture accretion rate is wider than that at t = 0.2 s.

To conclude, the global structure of the flow can be described approximately by the standard NDAF for a long timescale 0.1 s ≲ t ≲ 2 s.

The top panel of Figure 7 shows the accretion rate of the torus material onto the MNS, which is defined by the radial mass inflow rate (i.e., Equation (32)) at x = 20 km. The mass accretion rate is $\sim 0.5\,{M}_{\odot }\,{{\rm{s}}}^{-1}$ at ∼0.1 s, which is consistent with the estimation in Equation (22). The mass accretion rate monotonically decreases, and it becomes $\sim 0.004\,{M}_{\odot }\,{{\rm{s}}}^{-1}$ at t = 2 s for the fiducial model. This decrease stems partly from the decrease of the torus mass, but the major reason is the radial expansion of the torus due to the outward angular momentum transport, as described in the previous subsection. Unlike the estimation of Equation (22), the mass accretion rate for t ≲ 1 s depends only weakly on α_vis (compare the results for α_vis = 0.01, 0.02, and 0.04). This is because the typical radius of the torus becomes large more rapidly for the model with higher viscosity parameters (see the bottom panel of Figure 4). The mass accretion rate for t ≳ 1 s depends on α_vis; the rate for the model with α_vis = 0.04 decreases more rapidly than that for the fiducial model. This is because the neutrino cooling in the torus becomes inefficient at the earlier time; hence, the torus material starts being ejected rather than accreted onto the MNS (see Section 3.4 for details).

Zoom In Zoom Out Reset image size

The middle panel of Figure 7 shows the time evolution of the neutrino emission rates. Because of the presence of the viscous heating, the neutrino emission rate for ${\alpha }_{\mathrm{vis}}\ne 0$ is higher than that of the inviscid model. For the viscous models, the high emission rate of $\sim {10}^{53}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ is sustained for ∼0.1 s. However, for the late time t ≳ 1 s, the emission rate decreases to $\sim {10}^{52}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ as the accretion rate, ${\dot{M}}_{{\rm{b}},\mathrm{MNS}}$, decreases. Since the mass accretion rate depends weakly on α_vis for t ≲ 1 s, the neutrino emission rate also depends weakly on the viscosity parameter.

The bottom panel of Figure 7 shows the baryon mass of the MNS, defined by

$$\begin{eqnarray}&&{M}_{{\rm{b}},\mathrm{MNS}}={\int }_{x\lt 20\mathrm{km}}{d}^{3}x\ {\rho }_{* }.\end{eqnarray} \tag{ 35 }$$

As found from the top panel of Figure 7, the MNS mass increases with time, but eventually, it is saturated because the mass accretion from the torus decreases significantly. The final relaxed value, $\approx 2.9\,{M}_{\odot }$, depends only weakly on α_vis.

For the DD2 EOS, the MNS with ${M}_{{\rm{b}},\mathrm{MNS}}=2.9\,\,{M}_{\odot }$ does not collapse to a black hole, at least in a few seconds after the onset of the merger, because the maximum gravitational mass for cold spherical neutron stars for this EOS is quite high, $2.42\,\,{M}_{\odot }$. However, for softer EOSs, in which the maximum gravitational mass is smaller, say $2.1\,\,{M}_{\odot }$, the remnant MNS of the baryon mass of $\sim 2.9\,\,{M}_{\odot }$ may collapse to a black hole due to the mass accretion.

3.4.1. Early Viscosity-driven Mass Ejection

Panel (a) of Figure 8 shows the evolution of the baryon mass of the ejecta, which is calculated by integrating the flux of the unbound material at an outer surface as

$$\begin{eqnarray}{M}_{{\rm{b}},\mathrm{ej}}(t) & = & {\displaystyle \int }_{0}^{t}{dt}^{\prime} \,2\displaystyle \int {{dS}}_{k}{\rho }_{* }{v}^{k}{\rm{\Theta }}(| {{hu}}_{t}| -1)\\ & = & {\displaystyle \int }_{0}^{t}{dt}^{\prime} \,2{\displaystyle \int }_{0}^{{L}_{\mathrm{bnd}}}{dx}\ 2\pi x{\rho }_{* }{v}^{z}{\rm{\Theta }}(| {{hu}}_{t}| -1)\\ & & +{\displaystyle \int }_{0}^{t}{dt}^{\prime} \,2{\displaystyle \int }_{0}^{{L}_{\mathrm{bnd}}}{dz}\ 2\pi x{\rho }_{* }{v}^{x}{\rm{\Theta }}(| {{hu}}_{t}| -1),\end{eqnarray} \tag{ 36 }$$

where dS_k is the two-dimensional area element on a cylinder of both radius and height L_bnd. Here we set L_bnd = 4000 km and suppose that fluid elements with $| {{hu}}_{t}| \gt 1$ are gravitationally unbound.⁶

Zoom In Zoom Out Reset image size

This figure shows that the mass ejection rate in the early phase t ≲ 0.2 s is rather high, 0.03–$0.1\,{M}_{\odot }\,{{\rm{s}}}^{-1}$. This is achieved by the early viscosity-driven ejecta (see Section 3.1).

We also calculate the evolution of the total ejecta energy by

$$\begin{eqnarray}&&{E}_{\mathrm{tot},\mathrm{ej}}(t)={\int }_{0}^{t}{dt}^{\prime} \int {{dS}}_{k}{\rho }_{* }{\hat{e}}_{0}{v}^{k}{\rm{\Theta }}(| {{hu}}_{t}| -1).\end{eqnarray} \tag{ 37 }$$

Using E_tot,ej and M_b,ej, the kinetic energy of the ejecta is defined by

$$\begin{eqnarray}&&{E}_{\mathrm{kin},\mathrm{ej}}={E}_{\mathrm{tot},\mathrm{ej}}-{M}_{{\rm{b}},\mathrm{ej}},\end{eqnarray} \tag{ 38 }$$

which is shown in panel (b) of Figure 8. Here we assumed that the internal energy of the ejecta would be totally transformed into the kinetic energy during the subsequent expansion. Unlike the mass accretion rate and the neutrino emission rate shown in Section 3.3, the mass and kinetic energy of the ejecta increase monotonically with α_vis for $t\lesssim 0.2\,{\rm{s}}$.

Panel (c) shows the evolution of the average velocity of the ejecta, which is defined by ${V}_{\mathrm{ej}}=\sqrt{2{E}_{{\rm{k}},\mathrm{ej}}/{M}_{{\rm{b}},\mathrm{ej}}}$. For the early viscosity-driven ejecta, which is ejected for t ≲ 0.2 s, the average velocity is in the range 0.15–0.2 c and monotonically increases with α_vis as the mass and kinetic energy of the ejecta. This is because the shock driven by the early viscous effect on the MNS is stronger for the higher viscosity parameter model.

In our simulations, we suppose that the viscosity suddenly arises after the merger remnant settles to a quasi-stationary state. We should note that the variation of the quasi-equilibrium configuration occurs only if the viscous effect is enhanced after the dynamical phase of the merger ends, e.g., by MHD turbulence induced by the significantly amplified magnetic fields at the onset of the merger (Kiuchi et al. 2014, 2015b, 2017). If the viscosity in the MNS arises slowly enough, the variation of the quasi-stationary states occurs smoothly, so that the density wave would not appear strongly at the MNS surface and the mass ejection due to the shock would not occur significantly.

3.4.2. Late-time Viscosity-driven Mass Ejection

The mass ejection for t ≳ 0.2 s is driven primarily by the viscous effects in the torus. The ejecta is launched mainly toward the polar region for the relatively early time with $t\lesssim 1\,{\rm{s}}$, as found in panels (4)–(6) of Figure 2. However, the major mass ejection mechanism is subsequently changed as discussed below. In the following, we focus only on the models for α_vis = 0.02 and 0.04 because long-term simulations are performed only for these models.

Figure 9 shows the evolution of the mass, kinetic energy, and average velocity of the ejecta that becomes unbound after t = 0.5 s for the components toward the polar direction of 0° ≤ θ ≤ 30° and the other (equatorial) direction of 30° < θ ≤ 90°, respectively. From this figure, it is found that the material ejected toward the polar direction is smaller than or approximately as large as that ejected toward the equatorial direction. On the other hand, the polar ejecta has larger kinetic energy than the equatorial ejecta because the average velocity of the polar ejecta is larger (∼0.15 c) than that of the equatorial ejecta (∼0.05 c). This clearly shows that there are two components for the ejecta in the late time. We note that the two components are distinct because they are launched from completely different regions in the system. The polar ejecta is launched from the vicinity of the MNS, while the equatorial ejecta is launched from the outer region of the torus.

Zoom In Zoom Out Reset image size

As found in the top panel of Figure 9, for α_vis = 0.04, the mass ejection in this late phase (for t ≳ 0.5 s) is primarily toward the equatorial direction, while for the fiducial model (α_vis = 0.02), the equatorial mass ejection is as weak as the polar one up to t ∼ 2.7 s. For α_vis = 0.04, the early viscosity-driven mass ejection continues for a later time than that for the fiducial (α_vis = 0.02) model because the shock driven by the early viscous effect on the MNS is stronger for the larger viscosity parameter models. This is the reason for the larger mass ejection rate toward the equatorial direction for α_vis = 0.04 than that for the fiducial model for t ≲ 0.8 s. On the other hand, the reason for the rapid increase of the mass ejection rate for t ≳ 0.8 s is the viscous effect on the torus. The torus expands due to the viscous angular momentum transport; hence, its density and temperature decrease. As a result, the neutrino cooling becomes inefficient in the outer region of the torus. Then, the materials in the outer region can be ejected toward the equatorial region by the viscous heating and angular momentum transport without suffering from the neutrino cooling.

This late-time viscosity-driven mass ejection from the torus is first suggested in Fernández & Metzger (2013) and Metzger & Fernández (2014) and is also found in Just et al. (2015). The velocity of the ejecta in our work, ∼0.05 c, as well as the mass ejection mechanism, is totally consistent with their results. We confirm their results in a more realistic setup with an initial condition derived from three-dimensional simulation by a fully general relativistic simulation that self-consistently solves the evolution of both the MNS and the torus surrounding the MNS.

After the neutrino cooling becomes subdominant in the torus, the torus expands in the viscous timescale defined by Equation (21); hence, this late-time mass ejection rate can be estimated by dividing the torus mass by the viscous timescale as

$$\begin{eqnarray}\dot{M} & \sim & \displaystyle \frac{{M}_{{\rm{b}},\mathrm{torus}}}{{t}_{\mathrm{vis}}}\\ & \sim & 0.012\,{M}_{\odot }\,{{\rm{s}}}^{-1}\left(\displaystyle \frac{{\alpha }_{\mathrm{vis}}}{0.04}\right){\left(\displaystyle \frac{{H}_{\mathrm{tur}}}{10\mathrm{km}}\right)}^{2}\\ & & \times \left(\displaystyle \frac{{M}_{{\rm{b}},\mathrm{torus}}}{0.05\,{M}_{\odot }}\right){\left(\displaystyle \frac{{R}_{\mathrm{torus}}}{100\mathrm{km}}\right)}^{-7/2}{\left(\displaystyle \frac{{M}_{\mathrm{MNS}}}{2.6{M}_{\odot }}\right)}^{1/2},\end{eqnarray} \tag{ 39 }$$

where we used the torus mass for $t\gtrsim 1\,{\rm{s}}$ (see the top panel of Figure 4). This mass ejection rate agrees reasonably well with our results (see panel (d) of Figure 8).

For the fiducial model, this late-time equatorial viscosity-driven mass ejection also sets in at t ∼ 2.7 s. As found in the last panel of Figure 2, the velocity in the dense region (ρ ∼ 10⁷ g cm⁻³) turns outward. In addition, the top panel of Figure 9 indicates that the equatorial mass ejection rate increases for t ≳ 2.7 s. The time delay for the onset of the late-time viscosity-driven mass ejection toward the equatorial plane for the smaller value of α_vis is simply due to the less viscous power.

The mass ejection efficiency, defined by ${\dot{M}}_{{\rm{b}},\mathrm{ej}}/{\dot{M}}_{{\rm{b}},\mathrm{MNS}}$, is shown in panel (e) of Figure 8. We find that this ratio becomes larger than unity for $t\gtrsim 1\,$ and $t\gtrsim 2.3\,{\rm{s}}$ for α_vis = 0.04 and 0.02, respectively, and it eventually exceeds 10 because the mass accretion onto the MNS approaches zero. Therefore, a large fraction of the torus material, which is $\sim 0.05\ \,{M}_{\odot }$ at t ≳ 1 s (see the top panel of Figure 4), is likely to be ejected from the system eventually. This speculation is supported by the result for the model with α_vis = 0.04, for which the mass of the late-time equatorial viscosity-driven ejecta is already $\sim 0.05\,{M}_{\odot }$ at t ∼ 2.7 s.

For the α = 0.02 model, a large fraction of the material that would eventually be ejected is still in the region 0 ≤ x < L_bnd and 0 ≤ z < L_bnd; hence, they are still not considered an ejecta at the end of the simulations. The main neutrino emitter is the MNS for t ≳ 0.5 s, and the neutrino emission rate is ∼10⁵² erg s⁻¹. For this emission rate, the timescale for the neutrino absorption

$$\begin{eqnarray}&&\displaystyle \frac{4\pi {r}^{2}\langle \omega \rangle }{{L}_{\nu }}\displaystyle \frac{1}{{G}_{{\rm{F}}}^{2}{\langle \omega \rangle }^{2}}\\ &&\quad \sim 0.4\,{\rm{s}}{\left(\displaystyle \frac{\langle \omega \rangle }{10\mathrm{MeV}}\right)}^{-1}{\left(\displaystyle \frac{{L}_{\nu }}{{10}^{52}\mathrm{erg}{{\rm{s}}}^{-1}}\right)}^{-1}{\left(\displaystyle \frac{r}{100\mathrm{km}}\right)}^{2}\end{eqnarray} \tag{ 40 }$$

is much longer than the time for the expansion timescale

$$\begin{eqnarray}&&r/v\sim 7\times {10}^{-3}\,{\rm{s}}\left(\displaystyle \frac{r}{100\,\mathrm{km}}\right){\left(\displaystyle \frac{v}{0.05c}\right)}^{-1}.\end{eqnarray} \tag{ 41 }$$

Thus, the electron fraction of the material in the outer part of the torus, which is in the range 0.3–0.4, is frozen out; hence, the electron fraction of the material that would eventually be ejected is also in the range 0.3–0.4.

Figure 10 shows the density profiles on the rotational axis (left panel) and equatorial plane (right panel) for the fiducial model. We find that the density decreases approximately in proportion to z⁻² along the rotational axis. This behavior does not depend strongly on the time, since the polar mass ejection develops from the early phase of the evolution of the system. On the other hand, along the equatorial plane, the density for the large radius (x ≳ 500 km) decreases in proportion approximately to x⁻³ after the equatorial viscosity-driven mass ejection sets in (see the blue curve in the right panel). This radius dependence does not change significantly in time for the late phase (t ≳ 2.7 s); hence, we expect that the density structure of the equatorial ejecta is also proportional approximately to x⁻³. This radius dependence of the equatorial ejecta is similar to that of the dynamical ejecta with a velocity ≲0.4 c (Nagakura et al. 2014).

Zoom In Zoom Out Reset image size

To summarize this subsection, we find two mass ejection mechanisms for the system composed of the long-lived MNS and torus. One is the early viscosity-driven mass ejection, in which the differential rotation of the MNS is the engine, and the other is the late-time viscosity-driven mass ejection from the torus, in which the differential rotation of the torus is the engine. Specifically, there are two different components for the late-time viscosity-driven ejecta from the torus, one of which is the viscosity-driven ejecta toward the polar direction in assistance with the neutrino irradiation, and the other is launched primarily toward the equatorial direction in the later phase for which the neutrino cooling becomes inefficient because the temperature of the torus decreases sufficiently (see Table 3 for a summary).

Table 3.  Properties of the Ejecta for the DD2 EOS

Type of Ejecta	Mass (M⊙)	Vej/c	Ye	Direction	Duration
Dynamical ejecta	O(10−3)	∼0.2	0.05–0.5	θ ≳ 45°	t – tmerge ≲  10 ms
Early viscosity-driven ejecta	$\sim {10}^{-2}({\alpha }_{\mathrm{vis}}/0.02)$	∼0.15 − 0.2	0.2–0.5	$\theta \gtrsim 30^\circ$	$t-{t}_{\mathrm{merge}}\lesssim 0.1\,{\rm{s}}$
Late-time viscosity-driven ejecta (polar)	∼10−3 (tν/s)	∼0.15	0.4–0.5a	θ ≲ 30°	t − tmerge ∼ tν ∼ 10 s
Late-time viscosity-driven ejecta (equatorial)	≳10−2	∼0.05	0.3–0.4a	θ ≳ 30°	t − tmerge ∼ 1–10 s

Note. The table shows the components of ejecta, ejecta mass, average velocity, electron fraction, direction of the mass ejection, and ejection duration. Here t_ν and t − t_merge denote the duration of the neutrino emission and the time after the onset of the merger, respectively.

^aWe note that for the EOS in which the value of M_max is not as high as that for the DD2 EOS, the MNS could collapse to a black hole within ∼1 s. Even for such an EOS, the late-time viscosity-driven mass ejection should continue after the black hole formation, but because of shorter neutrino irradiation time, the value of Y_e is likely to be smaller than 0.3.

Download table as:  ASCIITypeset image

Here we discuss the properties of the ejecta for individual mass ejection processes based on the results of our numerical-relativity simulations.

We summarize possible mass ejection processes in the merger and post-merger phases in Figure 11. First, during the merger, the dynamical mass ejection occurs (e.g., Hotokezaka et al. 2013b; Sekiguchi et al. 2015, 2016). This mass ejection proceeds for ∼10 ms primarily toward the equatorial direction. The ejecta mass depends on the stiffness of the neutron star EOS, total mass, and mass ratio of the binary, and it is in the range 0.001–$0.02\,\,{M}_{\odot }$. For the DD2 EOS, which is used in this work, it is 0.002–0.005 $\,{M}_{\odot }$ for the total mass of $2.7\ \,{M}_{\odot }$ with the mass ratio between 0.85 and 1 (Sekiguchi et al. 2016). In this work, we consider the equal-mass merger of $1.35\ \,{M}_{\odot }$ neutron stars, in which the dynamically ejected mass is $\sim 0.002\ \,{M}_{\odot }$. In contrast to the ejecta mass, the typical velocity and the profile of the electron fraction of the dynamical ejecta depend only weakly on the neutron star EOS, and they are in the ranges 0.15–$0.25c$ and 0.05–0.5, respectively (Sekiguchi et al. 2015).

Zoom In Zoom Out Reset image size

In the post-merger phase, an MNS surrounded by a torus is typically formed. However, the long-term evolution process of this system depends on the neutron star EOS and total mass of the binary. If the EOS is stiff enough (i.e., the maximum mass for cold spherical neutron stars, M_max, is high enough) or the total mass of the binary is small enough, the MNS survives for seconds or longer after the onset of the merger as shown in this paper (Case I of Figure 11). Otherwise, the MNS collapses to a black hole after it loses its thermal energy and/or angular momentum (Case II of Figure 11). The time at which the collapse takes place would be determined by the value of M_max.

The magnetic field strength inside the remnant MNS is likely to be significantly enhanced due to the Kelvin–Helmholtz instability in the shear layer at the onset of the merger (Kiuchi et al. 2015a, 2017). Then, it is natural to suppose that MHD turbulence would be induced, and, consequently, the MHD-driven viscosity is likely to arise. If the lifetime of the MNS is sufficiently long, i.e., ≳ several tens of ms, the early viscosity-driven mass ejection could occur due to an induced sound wave and resulting shock wave generated associated with the variation of the density profile of the MNS caused by the angular momentum transport driven by the MHD turbulence. Our present study indicates that the early viscosity-driven mass ejection proceeds for $\lesssim 0.1\,{\rm{s}}$ in a moderately isotropic manner with $\theta \gtrsim 30^\circ$. The mass of this ejecta depends on the magnitude of the turbulent viscosity, but the mass of $\sim 0.01\ \,{M}_{\odot }$ could be ejected for a reasonable value of the viscosity parameter of α_vis = 0.01–0.04 (Kiuchi et al. 2017). For the fiducial model (α_vis = 0.02), the mass of the early viscosity-driven ejecta is $\sim 0.01\ \,{M}_{\odot }$. We note that this mass can be larger if the initial torus mass is larger. The typical velocity is in the range 0.15–0.2 c.

The left two panels of Figure 12 show the density and electron fraction profiles of the ejecta at t = 0.1 s for α_vis = 0.02 and 0.04. For both models, the low-electron fraction ejecta exists near the equatorial plane (region 5 in the left panels of Figure 12). This is composed partly of the dynamical ejecta and mainly of the early viscosity-driven ejecta. Since the power for the early viscosity-driven mass ejection is higher for α_vis = 0.04, the low-electron fraction ejecta is located at a more distant zone for this model. We find that the electron fraction of the early viscosity-driven ejecta is in the range 0.2–0.5 with the peak at Y_e ∼ 0.3–0.4 (see the right panels of Figure 12).

Zoom In Zoom Out Reset image size

In the presence of a long-lived MNS, the late-time viscosity-driven mass ejection occurs toward the polar direction for t ≳ 0.2 s. The top left and middle panels of Figure 13 show the density and electron fraction profiles of this polar ejecta for the fiducial model at t = 1 s. The ejecta for θ ≲ 30° is launched from the region near the central MNS in assistance with the neutrino irradiation toward the polar direction (regions 1 and 2 in the left panels of Figure 13); hence, the velocity of the ejecta is ∼0.15 c. Due to the strong neutrino irradiation, the electron fraction of the ejecta is increased to be quite high as 0.4–0.5. The properties of the polar ejecta for the model with α_vis = 0.04 are similar to those of the ejecta for the fiducial model.

Zoom In Zoom Out Reset image size

The late-time viscosity-driven mass ejection from the expanded torus occurs as we found for the model with α_vis = 0.04 (see the bottom left and middle panels of Figure 13). Even for the lower viscosity parameter models, this type of ejecta is likely to be induced for a later phase, as discussed in Section 3.4. This mass ejection is likely to occur irrespective of the presence of the long-lived MNS, and the velocity of this ejecta is ∼0.03−0.05 c (Metzger & Fernández 2014). This velocity is appreciably smaller than that for the polar ejecta component because it is launched from the outer region of the torus for which the typical velocity scale is <0.1 c. In the presence of the long-lived MNS, the effect of the neutrino irradiation is significant enough to increase the electron fraction above ∼0.3 for this ejecta. The bottom middle panel of Figure 13 shows that the late-time viscosity-driven ejecta is moderately neutron-rich as Y_e = 0.3–0.4 for $\theta \gtrsim 30^\circ$. In the absence of the long-lived MNS, the neutrino irradiation is quite minor; hence, the electron faction of the ejecta would be relatively low, as shown by Metzger & Fernández (2014) and Just et al. (2015). Note that for a shorter lifetime of the MNS with ≲1 s, the amount of the neutron-rich matter with Y_e < 0.3 would be larger than what we found.

Because the mass accretion onto the MNS is significantly suppressed for the late time $t\gtrsim 1\,{\rm{s}}$, an appreciable fraction of the torus material is likely to be ejected as the late-time viscosity-driven ejecta. The torus mass for the late time ($t\gtrsim 1$ s) is $\sim 0.05\,\,{M}_{\odot }$ (see the top panel of Figure 4); hence, the mass of this late-time ejecta could be $\sim 0.05\,\,{M}_{\odot }$. We note that this mass can be larger if the initial torus mass is larger.

In Table 3, we summarize the mass, typical velocity, electron fraction, direction of the ejection, and duration for each ejecta component. This is likely to show a universal picture of the mass ejection process from the merger remnant for the case that a long-lived MNS is formed (Case I in Figure 11). On the other hand, if a long-lived MNS is not formed after the merger (Case II in Figure 11), early viscosity-driven ejecta and late-time polar viscosity-driven ejecta would be minor, and the late-time equatorial viscosity-driven ejecta would have a much lower value of Y_e. We also note that if the lifetime of the MNS is shorter than $\lesssim 1\,{\rm{s}}$, the value of Y_e for the late-time viscosity-driven ejecta could be smaller than 0.3 (Metzger & Fernández 2014; Lippuner et al. 2017).

At the end of this subsection, we emphasize that the viscous hydrodynamics we employ in this work is an effective approach to take into account the angular momentum transport and heating due to the MHD turbulence; hence, the values of the viscosity parameter we assumed in this work should be justified by performing very high-resolution MHD simulations that resolve a variety of MHD instabilities in the future.

Here we discuss the elemental abundance in the post-merger ejecta. We note that we do not consider the dynamical ejecta in this subsection.

Figure 14 shows the mass histogram of the ejecta as a function of the electron fraction and specific entropy at t = 1 s for α_vis = 0.02 and 0.04. As found in this figure, the electron fraction of the ejecta is widely distributed, but the mass of the ejecta component with Y_e ≲ 0.25 is minor. In this ejecta, it is expected that the r-process elements heavier than the second peak, including lanthanide elements, are not significantly synthesized (e.g., Wanajo et al. 2014; Martin et al. 2015; Tanaka et al. 2018).

Zoom In Zoom Out Reset image size

Note that the mass histogram of the ejecta for α_vis = 0.04 exhibits remarkable excess in Y_e ≈ 0.3–0.4 and s/k_B ≲ 10, compared to those for α_vis = 0.02. The reason for this is that for α_vis = 0.02, the ejecta at t = 1 s is composed mainly of the early viscosity-driven and late-time polar viscosity-driven ejecta, while for α_vis = 0.04, it is additionally composed of the late-time equatorial viscous-driven ejecta. Thus, the excess found for α_vis = 0.04 indicates that the late-time equatorial viscosity-driven ejecta has Y_e = 0.3–0.4 and $s/{k}_{B}\lesssim 10$.

We performed a nucleosynthesis calculation as a post-process to the ejecta as done in Wanajo et al. (2014). In our scheme, the temporal evolution of temperature, density, and Y_e are obtained by using a tracer particle method (see Nishimura et al. 2015 for details and methodology). In this work, we employed a nuclear reaction network by Nishimura et al. (2016), of which the base theoretical nuclear mass formula is the Finite-Range Droplet Model (FRDM; Möller et al. 1995).

Figure 15 shows the mass fraction of the nuclei as a function of atomic number Z for five angle bins shown in Figures 12 and 13 for the model with α_vis = 0.04, in which the late-time equatorial viscosity-driven ejecta is found clearly. As expected from the Y_e histogram, the r-process nucleosynthesis does not proceed sufficiently, and, remarkably, the mass fraction of lanthanide elements is very small. Especially, the material ejected to angular regions 1–4 in Figure 12 (i.e., θ < 60°) is approximately lanthanide-free. In Table 4, the mass fraction of lanthanide and actinide elements in the individual angular regions is shown. We found that the mass fraction is ≲10⁻⁷ for angular regions 1–4. If the ejecta is contaminated only minorly by the lanthanide elements, the opacity of the ejecta is ≪10 cm² g⁻¹ and would be ∼0.1–1 cm² g⁻¹ (Tanaka et al. 2018). According to the standard macronova/kilonova model (Li & Paczyński 1998; Metzger et al. 2010), if we observe the post-merger ejecta directly, the time to reach the peak emission, peak luminosity, and its effective temperature are estimated to give

$$\begin{eqnarray}{t}_{\mathrm{peak}} & \approx & \sqrt{\displaystyle \frac{\kappa {M}_{\mathrm{ej}}}{4\pi {{cV}}_{\mathrm{ej}}}\xi }\approx 1.5\,\mathrm{days}\ {\left(\displaystyle \frac{{V}_{\mathrm{ej}}}{0.1c}\right)}^{-1/2}\\ & & \times {\left(\displaystyle \frac{{M}_{\mathrm{ej}}}{0.03{M}_{\odot }}\right)}^{1/2}{\left(\displaystyle \frac{\kappa }{0.3{\mathrm{cm}}^{2}{{\rm{g}}}^{-1}}\right)}^{1/2}{\xi }^{1/2},\end{eqnarray} \tag{ 42 }$$

$$\begin{eqnarray}{L}_{\mathrm{peak}} & \approx & \displaystyle \frac{{{fM}}_{\mathrm{ej}}{c}^{2}}{{t}_{\mathrm{peak}}}\approx 4.3\times {10}^{41}\mathrm{erg}\,{{\rm{s}}}^{-1}\ \left(\displaystyle \frac{f}{{10}^{-6}}\right){\left(\displaystyle \frac{{V}_{\mathrm{ej}}}{0.1c}\right)}^{1/2}\\ & & \times {\left(\displaystyle \frac{{M}_{\mathrm{ej}}}{0.03{M}_{\odot }}\right)}^{1/2}{\left(\displaystyle \frac{\kappa }{0.3{\mathrm{cm}}^{2}{{\rm{g}}}^{-1}}\right)}^{-1/2}{\xi }^{-1/2},\end{eqnarray} \tag{ 43 }$$

$$\begin{eqnarray}{T}_{\mathrm{eff},\mathrm{peak}} & \approx & {\left[\displaystyle \frac{{L}_{\mathrm{peak}}}{4\pi {({V}_{\mathrm{ej}}{t}_{\mathrm{peak}})}^{2}{\sigma }_{\mathrm{SB}}}\right]}^{1/4}\approx 8\times {10}^{3}\,{\rm{K}}\ {\left(\displaystyle \frac{f}{{10}^{-6}}\right)}^{1/4}\\ & & \times {\left(\displaystyle \frac{{V}_{\mathrm{ej}}}{0.1c}\right)}^{1/8}{\left(\displaystyle \frac{{M}_{\mathrm{ej}}}{0.03{M}_{\odot }}\right)}^{-1/8}{\left(\displaystyle \frac{\kappa }{0.3{\mathrm{cm}}^{2}{{\rm{g}}}^{-1}}\right)}^{-3/8}{\xi }^{-3/8},\end{eqnarray} \tag{ 44 }$$

where f is the radioactive energy deposition factor,⁷ κ is the opacity of the material, ξ ( ≤ 1) is a geometric factor, and σ_SB is the Stefan–Boltzmann constant. Here we suppose that the average velocity is 0.1 c. These estimates show that if we observe this post-merger ejecta directly (i.e., from a low opening angle θ ≲ 45°), the electromagnetic signal would be of a short-timescale, high-luminosity, and blue transient.

Zoom In Zoom Out Reset image size

We note that in the ejecta component for which Y_e ≳ 0.4, the nucleosynthesis products are likely to have a smaller heating rate (or smaller value of f). Wanajo et al. (2014) showed that the specific heating rate for the ejecta with Y_e ≳ 0.35 is much smaller than that of more neutron-rich ejecta. Thus, the high-electron-fraction ejecta may play a minor role as the energy source of the electromagnetic signal (i.e., f could be much smaller than 10⁻⁶) even if the ejecta mass is much larger than that of the dynamical ejecta. We should also note that Wanajo et al. (2014) considered only dynamical ejecta. The late-time ejecta irradiated by neutrinos has higher entropy than the dynamical ejecta. Even in the high-electron-fraction material, heavy elements can be synthesized if the material has sufficiently high entropy and expansion velocity (see, e.g., Hoffman et al. 1997).

As we found in this paper, the early viscosity-driven ejecta and late-time equatorial viscosity-driven ejecta could have large mass $\gtrsim 0.01\,\,{M}_{\odot }$ and moderately small values of Y_e (0.2–0.5 and 0.3–0.4, respectively). Such ejecta is likely to be the major heating source and contribute to electromagnetic counterparts as the energy source. Note that the mass of the late-time polar viscosity-driven ejecta is likely to be much smaller, and Y_e is large as ≳0.4; hence, their contribution would be minor.

We note that the maximum mass for cold spherical neutron stars for the DD2 EOS is M_max ≈ 2.4 M_⊙. If the value of M_max is not as high as this value for the neutron star EOS in nature, the MNS could collapse into a black hole in a few seconds after the merger. In this case, the electron fraction of the late-time equatorial viscosity-driven ejecta becomes lower than that in the presence of the MNS (Metzger & Fernández 2014; Just et al. 2015; Lippuner et al. 2017); hence, the viscosity-driven ejecta would be more lanthanide-rich. We plan to perform simulations for such an EOS in future work.

In our simulation, we make several approximations and assumptions that could affect the results. In Just et al. (2015), the neutrino energy density and flux obtained by the M1 scheme are compared to those calculated by a ray-tracing method for a black hole–torus system. Their result suggests that the neutrino energy density and flux are overestimated by a factor of ≈2 in the polar direction (with the polar angle θ ≲ 30°) and underestimated in the more equatorial direction. In the phase where the neutrino emission from the torus is stronger than that of the MNS, the profile of the neutrino emissivity is nonspherical, so that the neutrino energy density would be overestimated in the polar region, which leads to the overestimation of the neutrino heating rate. Thus, in the earlier phase of the evolution with t ≲ 0.5 s, during which the torus dominates over the whole neutrino emission (see Figure 7), the polar ejecta may be affected by the behavior of the M1 scheme.

However, the outflow launched toward the polar direction is initially slow (see Figure 7 in Fujibayashi et al. 2017), so that the electron fraction of the ejecta in the polar region is settled into an equilibrium value by the neutrino absorption (Qian & Woosley 1996),

$$\begin{eqnarray}&&{Y}_{e,\mathrm{eq}}\approx {\left(1+\displaystyle \frac{{L}_{{\bar{\nu }}_{e}}}{{L}_{{\nu }_{e}}}\displaystyle \frac{{\epsilon }_{{\bar{\nu }}_{e}}-2{\rm{\Delta }}}{{\epsilon }_{{\nu }_{e}}+2{\rm{\Delta }}}\right)}^{-1},\end{eqnarray} \tag{ 45 }$$

where L_νe and _νe are the luminosity and average energy of electron neutrinos, ${L}_{{\bar{\nu }}_{e}}$ and ${\epsilon }_{{\bar{\nu }}_{e}}$ are those of electron antineutrinos, and Δ ≈ 1.293 MeV is the mass difference between neutron and proton. The equilibrium value is unchanged if the ratio of the fluxes of electron neutrinos and antineutrinos and average energy of neutrinos are the same. Thus, the effects of the M1 scheme on the electron fraction of the polar ejecta would be minor if the overestimation of the energy density arises for electron neutrinos and antineutrinos at the same level. On the other hand, in the later phase (t ≳ 0.5 s), the MNS dominates over the whole neutrino emission of the system; hence, the neutrino emissivity profile becomes more spherical. Thus, the artificial effects due to the M1 scheme would be minor in this phase.

We should note that for the neutrino pair annihilation, the heating rate also depends on the angular distribution of the neutrinos; hence, its heating rate would be affected artificially by the M1 scheme. Figure 16 compares the heating rates due to the neutrino absorption and pair-annihilation processes. While the pair-annihilation heating dominates over the whole heating rate for t = 0.25 s, it becomes comparable to the neutrino absorption heating at z ≲ 20 km for later times due to the decrease of the neutrino emission rate. Thus, at least for the early phase (t ≲ 0.3 s), the outflow is primarily powered by the pair-annihilation heating. Pair-annihilation heating has not been considered in most recent works (e.g., Martin et al. 2015; Lippuner et al. 2017). However, the effect could be large, as found here; hence, we need to consider the pair-annihilation heating appropriately to obtain the properties of the ejecta quantitatively.

Zoom In Zoom Out Reset image size

Since we approximately determine the average energy of neutrinos for calculating the neutrino absorption rates, the uncertainty in the estimated energy would affect the equilibrium value of the electron fraction of the ejecta. If the average energy of neutrinos is higher than that estimated in our simulation, the equilibrium value of the electron fraction becomes lower because the mass difference between neutron and proton becomes unimportant for the higher average energy of neutrinos. Figure 17 shows the average energy of electron-type neutrinos along the rotational axis. This figure shows that the difference between the average energy estimated by Equations (40) and (41) in Fujibayashi et al. (2017) is ≲50%. Using Equation (45), the 50% difference in the average energy around 15 MeV leads to the difference of ≈0.06 in the equilibrium value of the electron fraction if the flux of the electron-type neutrinos is the same.

Zoom In Zoom Out Reset image size

We approximate viscous effects due to MHD turbulence by solving viscous hydrodynamics equations. However, if sufficiently strong magnetic fields are globally formed, the Lorentz force would accelerate the ejecta. This would happen in the polar region because of the low density. The left panel of Figure 10 shows that the density at z = 20 km on the rotational axis is ≈10⁷ g cm⁻³, which suggests that the Lorentz force would be important in the region in which the magnetic field strength exceeds 5 × 10¹⁴ G. Global magnetic fields could be formed when the outflow is driven because the field line is stretched in the outflow. Such an outflow is driven possibly by the viscous heating and/or neutrino heating. Thus, another ejecta component could be generated as the magnetically accelerated wind from the strongly magnetized MNS. Metzger et al. (2018) suggested that a significant mass ejection ($\gtrsim 0.01\,\,{M}_{\odot }$) is possible if the MNS is sufficiently magnetized and rapidly rotating. In addition, the density profiles shown in Figure 10 would be modified in the existence of the global magnetic fields (e.g., Metzger et al. 2007 for winds launched from magnetized neutron stars). We should be checking whether global magnetic fields are formed in the relevant timescale after the merger and whether they affect the properties of the ejecta by performing neutrino radiation MHD simulations.

The remnants of binary neutron star mergers have been considered to power relativistic jets that would drive short-duration gamma-ray bursts (Eichler et al. 1989; Nakar 2007). If a relativistic jet is launched from the merger remnant, the jet may inject a part of its kinetic energy into the surrounding ejecta and modify its expansion velocity (cocoon; Murguia-Berthier et al. 2014; Nagakura et al. 2014). Thus, the timescale of the macronova/kilonova emission would be modified to be shorter. The thermal energy in the cocoon injected from the jet would power another component of electromagnetic transient (cocoon emission; Nakar & Piran 2017). Since the post-merger ejecta is lanthanide-poor in our simulation, the timescale of the cocoon emission is likely to be short; hence, this could contribute to the electromagnetic signal in the early phase (≲day).