Postmerger Mass Ejection of Low-mass Binary Neutron Stars

We use a fully general-relativistic, neutrino-radiation, viscous-hydrodynamics code developed in our previous work (Fujibayashi et al. 2018). Einstein's equation is solved with a version of the puncture-Baumgarte–Shapiro–Shibata–Nakamura (BSSN) formalism (Shibata & Nakamura 1995; Baumgarte & Shapiro 1999; Marronetti et al. 2008) incorporating the Z4c prescription (Hilditch et al. 2013) for the constraint violation propagation. The quantities solved in this formalism are listed in Table 1. From the three-metric γ_ij and extrinsic curvature K_ij , we define ${\tilde{\gamma }}_{{ij}}={\gamma }^{-1/3}{\gamma }_{{ij}}$, W = γ^−1/6, K = γ^ij K_ij , ${\tilde{A}}_{{ij}}\,={\gamma }^{-1/3}({K}_{{ij}}-{\gamma }_{{ij}}K/3)$, and ${F}_{i}={\delta }^{{jk}}{\partial }_{j}{\tilde{\gamma }}_{{ik}}$. Here γ denotes the determinant of γ_ij and the determinant of ${\tilde{\gamma }}_{{ij}}$ is assumed to be unity. For the gauge conditions, we employ the dynamical lapse and shift gauge conditions described in Equations (1) and (2) in Fujibayashi et al. (2017). We introduce a new variable Θ, which indicates the degree of constraint violation, for constraint propagation in the Z4c formalism (Hilditch et al. 2013). For this change, the equations in the original BSSN formalism are slightly modified. For example, $\hat{K}\equiv K-2{\rm{\Theta }}$ obeys the same evolution equation as K (see, for example, Shibata 2016). We write the evolution equation of Θ as

$$\begin{eqnarray}\begin{array}{rcl}({\partial }_{t}-{\beta }^{k}{\partial }_{k}){\rm{\Theta }} & = & \displaystyle \frac{1}{2}\alpha \left[R-{\tilde{A}}_{{ij}}{\tilde{A}}^{{ij}}\right.\\ & & +\ \left.\displaystyle \frac{2}{3}(\hat{K}+2{\rm{\Theta }})-16\pi {\rho }_{{\rm{h}}}\right]{e}^{-{r}^{2}/{r}_{{\rm{o}}}^{2}},\end{array}\end{eqnarray} \tag{ 1 }$$

where ρ_h ≡ T^αβ n_α n_β , $r=\sqrt{{x}^{2}+{y}^{2}+{z}^{2}}$, and r_o ≈ 1000 km. Here we set κ₁ = κ₂ = 0 in the full description of Z4c in Hilditch et al. (2013). As Equation (1) shows, we incorporate the constraint propagation only for a near zone to avoid the accumulation of the constraint violation only in the strong field zone.

Table 1. List of the Quantities Evolved in Our Formalism

Notation	Definition	Reference
${\tilde{\gamma }}_{{ij}}={\gamma }^{-1/3}{\gamma }_{{ij}}$	Conformal three-metric	Shibata & Nakamura (1995)
W = γ−1/6	Conformal factor	Marronetti et al. (2008)
K = γij Kij	Trace of the extrinsic curvature Kij	Shibata & Nakamura (1995)
${\tilde{A}}_{{ij}}={\gamma }^{-1/3}({K}_{{ij}}-{\gamma }_{{ij}}K/3)$	Trace-free part of Kij	Equation (2.9) in Shibata & Nakamura (1995)
${F}_{i}={\delta }^{{jk}}{\partial }_{j}{\tilde{\gamma }}_{{ki}}$	Auxiliary variable	Equation (2.19) in Shibata & Nakamura (1995)
Θ	Hamiltonian constraint violation	Hilditch et al. (2013)

Download table as:  ASCIITypeset image

We adopt the sixth-order Kreiss–Oliger dissipation for all the geometrical variables except for the lapse function to avoid the numerical instability that could appear when employing a nonuniform grid in the numerical evolution. We solve Einstein's evolution equation in Cartesian coordinates and employ the so-called cartoon method to impose axisymmetric conditions for the geometrical quantities (Shibata 2000; Alcubierre et al. 2001). In this paper, we employ the fourth-order Lagrangian interpolation for the cartoon procedure as in the series of our recent papers (e.g., Fujibayashi et al. 2018, 2020).

We briefly describe the basic equations and our method to solve viscous-hydrodynamics equations with neutrino radiation transfer. We solve the following evolution equations in a two-dimensional (2D) axisymmetric manner. We decompose neutrinos into "streaming" and "trapped" components. Then we also decompose the total energy-momentum tensor of the matter, which includes 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{ 2 }$$

where T^αβ is the sum of the energy-momentum tensors for the fluid ${T}_{\ (\mathrm{fluid})}^{\alpha \beta }$ and trapped neutrinos ${T}_{({\nu }_{i},{\rm{T}})}^{\alpha \beta }$, and ${T}_{({\nu }_{i},{\rm{S}})}^{\alpha \beta }$ is that for streaming neutrinos. The index i specifies the neutrino species. Following our previous work (Fujibayashi et al. 2017, 2018), we consider three species of neutrinos: electron neutrinos ν_e, electron antineutrinos ${\bar{\nu }}_{{\rm{e}}}$, and the other neutrino species ν_x, which represents muon and tau neutrinos and those antineutrinos altogether.

The evolution equations for these energy-momentum tensors are

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

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

where ${Q}_{(\mathrm{leak}){\nu }_{i}}^{\alpha }$ denotes the "leakage" rate of the ith species of neutrinos. Note that this leakage rate includes the heating due to neutrino absorption and pair annihilation, as well as the cooling due to neutrino emission, as described in Fujibayashi et al. (2017). ⁸ To solve the evolution equations for streaming neutrinos, we employ Thorne's truncated moment formalism (Thorne 1981; Shibata et al. 2011a) with a closure relation (Levermore 1984; González et al. 2007). For trapped neutrinos, we employ a leakage-based scheme developed in Sekiguchi (2010). The detailed description of these schemes is found in Sekiguchi (2010) and Fujibayashi et al. (2017).

The energy-momentum tensor of a viscous fluid with trapped neutrinos is written as

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

where ρ is the baryon rest-mass density, h = 1 + ε + P/ρ is the specific enthalpy with the specific internal energy ε and the pressure P, u^α is the four-velocity of the fluid, ν is the viscous coefficient, and ${\tau }^{0}{}_{\alpha \beta }$ is a shear viscous stress tensor, which is a symmetric tensor satisfying the relation ${\tau }^{0}{}_{\alpha \beta }{u}^{\alpha }=0$.

Following Israel & Stewart (1979), we write the evolution equation for the viscous stress tensor as

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

with the shear tensor σ_αβ 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{ 7 }$$

Here, h_αβ  = g_αβ  + u_α u_β . ζ is a nonzero constant of (time)⁻¹ dimension, which is set to 3 × 10⁴ rad s⁻¹ in this work following our previous work (Fujibayashi et al. 2018). The value is approximately four times larger than the largest angular velocity of the system, ≈7000 rad s⁻¹. With this prescription, the tensor ${\tau }^{0}{}_{\alpha \beta }$ approaches the shear tensor σ_αβ in a sufficiently short timescale of ∼ζ⁻¹ ≈ 0.03 ms.

We can rewrite Equation (6) as

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

where ${\tau }_{\alpha \beta }\equiv {\tau }^{0}{}_{\alpha \beta }-\zeta {h}_{\alpha \beta }$. Thus, in addition to the hydrodynamics equations of Equation (3), we solve Equation (8) as a basic equation that describes the evolution of τ_αβ . The details of the formulation of our general-relativistic viscous-hydrodynamics code are found in Shibata et al. (2017b).

For the viscous coefficient, we use the so-called Shakura–Sunyaev description as

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

where c_s is the sound speed and H_tur is a typical scale height of the system (Shakura & Sunyaev 1973). We adopt H_tur = 10 km following Fujibayashi et al. (2018).

To study the property of the ejecta in numerical simulations, we have to identify the unbound material in a finite computational domain. This gives us a nontrivial question about how we should determine the condition for the unbound material. There are several possibilities for this condition. In our previous numerical simulations for binary NS mergers (Sekiguchi et al. 2015, 2016), we employed the condition

$$\begin{eqnarray}&&{u}_{t}+1\lt 0\end{eqnarray} \tag{ 10 }$$

for the entire fluid elements in the computational domain, assuming that the dynamical ejecta has a small internal energy, which is much smaller than the kinetic energy. In our previous axisymmetric simulation for the evolution of the merger remnants (Fujibayashi et al. 2018), we employed the condition

$$\begin{eqnarray}&&{{hu}}_{t}+1\lt 0\end{eqnarray} \tag{ 11 }$$

for the fluid elements that have outgoing velocity, i.e., we took into account the thermal contribution for identifying the ejecta because the kinetic energy per mass of the postmerger ejecta is much smaller than that of the dynamical ejecta.

In the present work, alternatively, we employ the following condition:

$$\begin{eqnarray}&&{{hu}}_{t}+{h}_{\min }\lt 0,\end{eqnarray} \tag{ 12 }$$

where h_min is the minimum value of the specific enthalpy in the chosen tabulated EOS (≈0.9987). The reason for this choice is that the atomic mass unit m_u is used as the mass per baryon in the EOS we employed, and thus, the specific enthalpy can be smaller than unity for the case where the actual mass per baryon is lighter than m _u and where the pressure contribution to the enthalpy is negligible.

We note that the material with the lowest specific enthalpy h_min is composed mainly of heavy nuclei in the EOSs employed in this study, in which nuclear statistical equilibrium (NSE) is assumed for the nuclear composition. Thus, in Equation (12), the energy released by the recombination of free nucleons into heavy nuclei is taken into account. That is, Equation (12) is the condition for the material being unbound only for the case where the temperature of the material is low enough that free nucleons recombine into heavy nuclei. One issue for using Equation (12) is that some fluid components, which stay in the central region and are strongly bounded by the central massive NS, can satisfy this condition. Such components remain composed of free nucleons, and thus, even if the material satisfies Equation (12), it is not ejected eventually. To remove the contribution of this component, we in addition impose the condition of r ≥ 500 km to identify the real ejecta component when we employ Equation (12) as the condition.

Obviously, the total mass of the ejecta becomes largest when we employ Equation (12). Thus, in our previous studies, in particular for the case where we employed Equation (10), the total mass of the ejecta would be underestimated (see Table 2). Also, the property of the ejecta such as the average electron fraction depends weakly on the choice of the condition for the unbound object.

Table 2. List of Models for the Three-dimensional Merger Simulation

Model	EOS	MNS	Rcold	tend−tmerge	Mej	$\langle {V}_{\mathrm{ej}}\rangle$	$\langle {s}_{\mathrm{ej}}\rangle$	$\langle {Y}_{{\rm{e}},\mathrm{ej}}\rangle$	Ntp	Xla
		(M⊙)	(km)	(ms)	(10−2 M⊙)	(c)	(kB)
DD2-125	DD2	1.25	13.14	48.5	0.10	0.19	27	0.20	2203	0.17
SFHo-125	SFHo	1.25	11.96	50.8	0.13	0.18	21	0.27	2511	0.072
DD2-135	DD2	1.35	13.18	61.7	0.15	0.19	25	0.23	2761	0.093
SFHo-135	SFHo	1.35	11.92	66.6	0.83	0.27	16	0.24	2197	0.031

Note. M_NS, R_cold, t_end−t_merge, M_ej, $\langle {V}_{\mathrm{ej}}\rangle$, $\langle {s}_{\mathrm{ej}}\rangle$, $\langle {Y}_{{\rm{e}},\mathrm{ej}}\rangle$, N_tp, and X_la are the mass of each NS, radius of the cold spherical NS with the mass of M_NS, time after merger at which the simulation is terminated, ejecta mass, their average velocity, entropy, and electron fraction, the number of tracer particles we used, and the mass fraction of lanthanides, respectively. The ejecta is defined as the material that satisfies Equation (10).

Download table as:  ASCIITypeset image

As in our previous work (Fujibayashi et al. 2017, 2018), we use hydrodynamical configurations obtained by three-dimensional (3D) numerical relativity simulations for equal-mass binary NS mergers as the initial condition for our 2D simulations. At several tens of milliseconds after the merger, the system, which is composed of a massive NS and an accretion disk, becomes nearly axisymmetric. We averaged the 3D hydrodynamical quantities at ≈50 ms after the merger over the azimuthal angle around the rotational axis to obtain a 2D axisymmetric hydrodynamical configuration used in our simulation. For the viscous tensor, we set ${\tau }^{0}{}_{\alpha \beta }=0$ initially. That is, we assume that the fluid is ideal initially. The initial condition for the geometrical variables is determined in the assumption of conformal flatness for γ_ij and by solving the Hamiltonian and momentum constraints for a given matter configuration.

In this work, we set the origin of time t = 0 to the beginning of our 2D simulation, which is ≈50 ms after the onset of the merger.

In this subsection, we briefly describe models for the 3D merger simulation. The 3D simulations are performed using our numerical relativity code with an approximate neutrino radiation transport (Sekiguchi et al. 2015, 2016). We employ the puncture-BSSN formalism for the evolution of the geometrical quantities and the leakage-based scheme described in Section 2.2 for the neutrino transport. In these 3D simulations, the neutrino absorption is taken into account as the heating source, but not the other processes such as the neutrino pair annihilation.

In this paper, we explore three models with different nuclear EOSs and total masses of the system as shown in Table 2. As mentioned in Section 1, we use two finite-temperature EOSs for nuclear matter, referred to as SFHo (Steiner et al. 2013) and DD2 (Banik et al. 2014), which predict the maximum mass of cold spherical NSs to be 2.06 M_⊙ and 2.42 M_⊙, respectively. We extend the tables of these EOSs to low-density and low-temperature ranges down to ≈0.16 g cm⁻³ and 10⁻³ MeV, respectively, using the EOS of ideal gas with arbitrarily relativistic and degenerate electrons (Timmes & Swesty 2000). We consider the equal-mass binaries with two sets of NS masses, 1.25 and 1.35 M_⊙.

In Table 2, the properties of the dynamical ejecta (defined by Equation (10)) are shown. They depend on the total mass of the binary and adopted EOS (Bauswein et al. 2013; Hotokezaka et al. 2013). The SFHo EOS is "softer" than DD2, in the sense that if we compare the radius of the NS with the same gravitational mass, the radius in the SFHo EOS is smaller than that in the DD2. The softer EOS generally enhances the shock-heating efficiency in the early evolution stage of the massive NS, which plays a role in increasing the fraction of unbound material. This leads to a larger ejecta mass for SFHo models than those for DD2 models for equal-mass binary mergers. In addition, the positron capture is enhanced in the shocked ejecta due to its high temperature, and thus the average values of the electron fraction for SFHo models are higher than those for DD2 models. For DD2 models, the mass ejection is driven predominantly by the tidal effect, which leads to the lower electron fraction of the ejecta. The dynamical mass ejection occurs primarily toward the direction of the equatorial plane, in particular for the low-mass cases, because the tidal effect is the main mass-ejection channel.

We note that the ejecta mass depends on its definition (see Section 2.3). For example, the ejecta mass for model DD2-125 evaluated using the condition given by Equation (12) with r ≥ 500 km is ≈0.002 M_⊙, which is approximately twice as large as that using the condition Equation (10) (≈0.001 M_⊙; see Table 2). On the other hand, the average electron fraction of the ejecta changes only slightly (by less than 5%).

The models for our 2D simulations in this study are listed in Table 3. For each model, the same EOS tables as those employed in the 3D simulation are used. To evolve the radiation viscous-hydrodynamics equations in the cylindrical coordinates (x, z), we employ the same nonuniform grid as that in our previous work (Fujibayashi et al. 2017, 2018), in which the grid structure is determined by the uniform grid spacing for the inner region (Δx₀), the range of the inner region (R_star), the increase rate of the grid spacing for the outer region (1 + δ), and the total grid number in each direction (N_grid). For all models in this study, we adopt R_star = 30 km and δ = 0.0075. We list Δx and N_grid for each model, together with the size of the computational domain, L, in Table 3.

Table 3. List of Models for the Two-dimensional Postmerger Simulation

Model	EOS	MNS	Δx0	Ngrid	L	αvis	tend	Mdisk	Mej,tot	$\langle {V}_{\mathrm{ej}}\rangle$	$\langle {s}_{\mathrm{ej}}\rangle$	$\langle {Y}_{{\rm{e}},\mathrm{ej}}\rangle$	Ntp	Xla
		(M⊙)	(m)		(km)		(s)	(10−2 M⊙)	(10−2 M⊙)	(c)	(kB)
DD2-125M	DD2	1.25	200	921	8411	0.04	6.2	32.8 (26.6)	11.4	0.089	15.0	0.322	39808	0.00441
SFHo-125H	SFHo	1.25	160	993	8745	0.04	3.8	14.8 (12.7)	6.1	0.108	18.7	0.342	24576	0.00261
DD2-135M	DD2	1.35	200	921	8411	0.04	6.1	24.8 (21.2)	8.6	0.098	17.3	0.337	39168	0.00276
DD2-125L a	DD2	1.25	250	857	8148	0.04	7.6	32.6 (26.3)	10.6	0.089	15.5	0.326	44416	0.00441
DD2-135L a	DD2	1.35	250	857	8148	0.04	6.9	24.3 (21.2)	8.5	0.095	18.8	0.346	48768	0.00378
DD2-125M-h b	DD2	1.25	200	921	8411	0.10	1.5	35.4 (31.6)	19.8	0.113	13.8	0.282	9856	0.0298
DD2-135M-v14 c	DD2	1.35	200	921	8411	0.04	2.8	27.9 (23.7)	8.7	0.113	20.4	0.341	⋯	⋯
DD2-135M-irr d	DD2	1.35	200	921	8411	0.04	5.6	24.7 (21.1)	6.9	0.106	16.3	0.287	⋯	⋯
DD2-135M-v0 e	DD2	1.35	200	921	8411	0.00	1.1	40.9 (11.0)	0.3	0.164	21.1	0.262	⋯	⋯

Notes. M_NS, Δx₀, N_grid, L, α_vis, t_end, M_disk, M_ej,tot, $\langle {V}_{\mathrm{ej}}\rangle$, $\langle {s}_{\mathrm{ej}}\rangle$, $\langle {Y}_{{\rm{e}},\mathrm{ej}}\rangle$, N_tp, and X_la are the mass of each NS, the grid spacing of the inner computational region, the grid number in each direction, the size of the computational domain, viscous parameter, time at which the simulation is terminated, initial disk mass, ejecta mass, its average velocity, entropy, electron fraction, the number of tracer particles we used, and the mass fraction of lanthanides, respectively. The initial disk mass is defined as the mass outside the NS radius R_NS (see Section 3.2 for its definition) at t = 10 ms. In the parentheses, the mass of the material with ρ < 10¹² g cm⁻³ is shown. For M_ej,tot, $\langle {V}_{\mathrm{ej}}\rangle$, $\langle {s}_{\mathrm{ej}}\rangle$, and $\langle {Y}_{{\rm{e}},\mathrm{ej}}\rangle$, the values at t = t_end are shown.

^a The model with a lower resolution. ^b The model with a higher viscous parameter. ^c The model without viscosity in the region with ρ > 10¹⁴ g cm⁻³. ^d The model without neutrino irradiation. ^e The model without viscosity.

Download table as:  ASCIITypeset image

Our fiducial models are DD2-125M, SFHo-125H, and DD2-135M. For DD2-125M and DD2-135M, we set Δx₀ = 200 m. On the other hand, for SFHo-125H, we employ a finer grid resolution with Δx₀ = 160 m than those in other models, because the stellar radii for the SFHo EOS are smaller than those for the DD2 EOS, and thus, the conservation of the total energy is less satisfied in a lower-resolution model. We also perform several other simulations: lower-resolution ones (DD2-125L and DD2-135L) with Δx₀ = 250 m in order to confirm the convergence of the numerical results with the adopted grid resolutions, one without the viscosity in the region for which the density is higher than 10¹⁴ g cm⁻³ (DD2-135M-v14) to clarify that the viscosity in the innermost region of the NS does not play an important role for the early viscosity-driven ejecta (Section 3.1), one without the viscosity in the whole region (DD2-135M-v0) to examine the effect of the viscosity on the evolution of the remnant and the ejecta, and one without neutrino heating/irradiation in the optically thin region (DD2-135M-irr) to understand the importance of neutrino absorption in the evolution of the remnant and the ejecta.

A number of magnetohydrodynamics simulations have shown that the effective value of α_vis becomes O(0.01) for disks rotating in the vicinity of central remnants (e.g., Kiuchi et al. 2018; Fernández et al. 2019). In addition, in our previous work (Fujibayashi et al. 2018), we found that (if α_vis ≤ 0.04) features of the postmerger ejecta do not depend strongly on the viscous parameter except for the mass-ejection timescale, which is approximately proportional to ${\alpha }_{\mathrm{vis}}{}^{-1}$. Therefore, we set α_vis = 0.04 for our fiducial models in this work. For the DD2-125 model, we also perform a simulation with a very high viscosity α_vis = 0.10 (DD2-125M-h) to explore the dependence of the result on the viscosity strength (i.e., the timescale of the viscous evolution of the system).

The evolutions of our grid-based hydrodynamical quantities should be converted into Lagrangian particle motions for the nucleosynthesis calculations (Section 4). For this purpose, we adopt a tracer particle method (Wanajo et al. 2014). As initial positions of the tracer particles, we select 128 points with polar angles in the range of θ = 0–π/2 on the arc with the radius of r_esc = 8000 km. We repeat the same procedure every Δt_set = 0.05 s to set tracer particles. The mass of each particle is set based on the radial mass flux of its initial position as

$$\begin{eqnarray}&&{\rm{\Delta }}m={r}_{\mathrm{esc}}{}^{2}{\rm{\Delta }}{\rm{\Omega }}{\rho }_{* }{v}^{r}{\rm{\Delta }}{t}_{\mathrm{set}},\end{eqnarray} \tag{ 13 }$$

where ΔΩ is the solid angle element, ρ_* = W⁻³ αu^t ρ, and v^r  = dr/dt. It depends on the polar angle and the time at which the particle is set initially. In this manner, the total mass of the particles, M_traj, is compatible with the ejecta mass defined in Section 3.3. We then evolve them backward in time. The numerical method of tracing particles is basically the same as that in Nishimura et al. (2015) except for the difference in the coordinate system employed. The time integration is based on the predictor-corrector method with second-order accuracy. We use bilinear interpolation to convert grid-based hydrodynamical quantities into those of the tracer particles. The particle position (x⁽ⁿ⁾, z⁽ⁿ⁾) at time t⁽ⁿ⁾ is obtained from

$$\begin{eqnarray}&&{x}^{(n)}={x}^{(n+1)}+{v}^{x* }{\rm{\Delta }}t,\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&{z}^{(n)}={z}^{(n+1)}+{v}^{z* }{\rm{\Delta }}t,\end{eqnarray} \tag{ 15 }$$

where ${v}^{x* }$ and ${v}^{z* }$ denote the modified velocity components of the particle estimated by the predictor-corrector method and Δt = t⁽ⁿ⁾ − t⁽ⁿ⁺¹⁾ (<0).

Although the qualitative evolution process of the postmerger remnant composed of a massive NS and an accretion disk was already described in our previous paper (Fujibayashi et al. 2018), we briefly summarize it again in this subsection. Figure 1 shows snapshots of the rest-mass density, electron fraction, temperature, and entropy per baryon for model DD2-125M. The top-left panel of Figure 1 shows the initial profile of the simulation. There is a dense disk extending to ∼100 km around the massive NS. The disk has high density (ρ ∼ 10¹² g cm⁻³), high temperature (k_B T ∼ 10 MeV), and low electron fraction (Y_e ≲ 0.1). The material of low electron fraction is located in the region with the density of ∼10⁶ g cm⁻³ around the equatorial plane. This material was expelled during the dynamical merger phase without experiencing significant positron capture by neutrons because of its low temperature and resulting in poor electron–positron pair creation. Part of this material is the dynamical ejecta, and the rest is still mechanically bound.

Zoom In Zoom Out Reset image size

The massive NS formed as a remnant after the binary NS merger is differentially rotating and it is first evolved by the viscous effect (before the disk is evolved). The timescale is estimated as

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{R}_{\mathrm{eq}}{}^{2}}{\nu } & \approx & 5\,\mathrm{ms}\,{\left(\displaystyle \frac{{\alpha }_{\mathrm{vis}}}{0.04}\right)}^{-1}{\left(\displaystyle \frac{{c}_{{\rm{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{array}\end{eqnarray} \tag{ 16 }$$

where R_eq is the equatorial radius of the NS, and we use Equation (9) for the viscous coefficient ν. In this timescale, the rotational profile of the massive NS is modified due to the viscous angular momentum transport. This produces the early viscosity-driven ejecta, powered by the kinetic energy of differential rotation of the NS, in particular in its outer region. ⁹

The range of the electron fractions of this component is approximately the same as that of dynamical ejecta (see Section 3.3.3). The amount of this component depends on the viscous parameter α_vis as found in Section 3.4.

In the polar direction, the high-entropy material is ejected mainly due to the heating by neutrinos emitted from the remnant (panels for t > 0 s of Figure 1). This polar mass ejection is efficient in particular for t ≲ 0.1 s during which time neutrinos are copiously emitted from the disk (see Section 3.2.3). After that, this mass ejection becomes weak with the decrease of the neutrino luminosity. However, the polar mass ejection still continues throughout our simulation time, because the massive NS is present and continuously emits neutrinos (see Figure 3).

The disk gradually expands due to the outward viscous angular momentum transport and viscous heating, while the temperature of the disk decreases because of the expansion and neutrino cooling. When the neutrino cooling rate becomes lower than the viscous heating rate, the disk material begins to be ejected because most of the viscous heating is used for the disk expansion. This late-time viscosity-driven ejecta consists mainly of the material with Y_e = 0.3–0.35 and s/k_B ≈ 10 (panels for t > 0 s of Figure 1; see also Section 3.3.3). The dynamics of the system among the explored models is quite similar to each other, although the properties of the ejecta depend quantitatively on the EOS, mass of the massive NS, and the magnitude of α_vis.

3.2.1. Radius of the NS Surface

In this section, we present the time evolution of quantities related to the massive NS and the disk surrounding it for each model. First, we describe how we distinguish the disk from the outer part of the massive NS.

It is not trivial to define the boundary that divides these two components, because the material of the massive NS is continuously connected to the disk (i.e., the density varies continuously). We use the angular velocity profile to define the NS and disk components, respectively. The top panel of Figure 2 shows the angular velocity profile along the equatorial plane at selected time slices for DD2-125M. At t = 0, the system including the high-density region of the remnant (x ≲ 15 km) is entirely in a differentially rotating state. The angular velocity is small in the central region, peaks at ≈10 km, and then decreases with the radius. This is the typical rotational profile of a merger remnant (e.g., Shibata et al. 2005). Within a few milliseconds after the simulation begins, the rotational profile of this region changes to become rigid due to a viscous angular momentum transport process, the timescale of which is estimated by Equation (16). This is clearly illustrated in the middle panel of Figure 2, which shows the logarithmic derivative of the angular velocity with respect to the cylindrical radius, i.e., $d\mathrm{ln}{\rm{\Omega }}/d\mathrm{ln}x$. For t > 0.01 s, the innermost region (with x ≲ 15 km) settles into a rigidly rotating state, i.e., $d\mathrm{ln}{\rm{\Omega }}/d\mathrm{ln}x\approx 0$. On the other hand, in the outer region with x ≳ 15 km, its value gradually approaches the Keplerian value of −1.5, which is the expected feature of the disk. We define the NS radius, R_NS, as the innermost cylindrical radius that satisfies $d\mathrm{ln}{\rm{\Omega }}/d\mathrm{ln}x=-0.1$ on the equatorial plane. Although the value is defined along the equatorial direction, we employ this in all polar angle directions to define the NS region (i.e., the NS region is assumed to be spherical).

Zoom In Zoom Out Reset image size

We note that defining the NS surface in terms of density is nontrivial particularly in an early stage of the evolution (t ≲ 0.3 s), wherein the density gradient at the NS surface is not sufficiently large. However, for t ≳ 0.3 s, the density can also be used to separate the two regions because the density gap between those regions becomes very large.

The bottom panel of Figure 2 shows the angular velocity evolution for model DD2-135M-v14, in which the viscosity is switched off for the region with the rest-mass density above 10¹⁴ g cm⁻³ to clarify that the inner region of the NS does not have a significant impact on the early viscosity-driven ejecta (for this model, $d\mathrm{ln}{\rm{\Omega }}/{dx}$ is positive or close to zero inside the NS). In this model, $d\mathrm{ln}{\rm{\Omega }}/d\mathrm{ln}x$ also steeply decreases for x ≳ 15 km and thus we can define the NS surface in the method mentioned above. On the other hand, the viscous angular momentum transport does not occur inside a radius of ≈10 km, in which the rest-mass density is higher than 10¹⁴ g cm⁻³.

We note that, as illustrated in Appendix B, the properties of the early viscosity-driven ejecta for this model are not significantly different from those of DD2-135M. This shows that the early ejecta component is driven by releasing the rotational kinetic energy near the NS surface (not the main body of the NS). We find that the slow increase of the angular velocity for t ≳ 0.1 s in this model is due to the convective motion inside the stellar body (see Appendix B).

3.2.2. Diagnostics

With the radius of the NS surface determined above, we define the baryon rest mass, angular momentum, and neutrino luminosity of the NS as

$$\begin{eqnarray}&&{M}_{\mathrm{NS}}={\int }_{r\leqslant {R}_{\mathrm{NS}}}{d}^{3}x\,{\rho }_{* },\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&{J}_{\mathrm{NS}}={\int }_{r\leqslant {R}_{\mathrm{NS}}}{d}^{3}x\,{\rho }_{* }{{hu}}_{\phi },\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&{L}_{\nu ,\mathrm{NS}}={\int }_{r\leqslant {R}_{\mathrm{NS}}}{d}^{3}x\,\alpha {W}^{-3}{Q}_{(\mathrm{leak})t},\end{eqnarray} \tag{ 19 }$$

where Q_(leak)t is the time component of the leakage rate (see Section 2.2). Here we multiply the lapse function α to take the gravitational redshift into account for the neutrino luminosity, and thus, this luminosity corresponds to the value observed at infinity. The quantities for the disk (M_disk, J_disk, and L_ν,disk) are also defined by integrating the same quantities in the region of r > R_NS. We also define the average electron fraction, average entropy, and total viscous heating rate of the disk by

$$\begin{eqnarray}&&{Y}_{{\rm{e}},\mathrm{disk}}=\displaystyle \frac{1}{{M}_{\mathrm{disk}}}{\int }_{r\gt {R}_{\mathrm{NS}}}{d}^{3}x\,{\rho }_{* }{Y}_{{\rm{e}}},\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&{s}_{\mathrm{disk}}=\displaystyle \frac{1}{{M}_{\mathrm{disk}}}{\int }_{r\gt {R}_{\mathrm{NS}}}{d}^{3}x\,{\rho }_{* }s,\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&{L}_{\mathrm{vis},\mathrm{disk}}={\int }_{r\gt {R}_{\mathrm{NS}}}{d}^{3}x\,\alpha {W}^{-3}\displaystyle \frac{1}{2}\rho h\nu {\tau }^{0\alpha \beta }{\tau }^{0}{}_{\alpha \beta }{u}_{t},\end{eqnarray} \tag{ 22 }$$

where $(1/2)\rho h\nu {\tau }^{0\alpha \beta }{\tau }^{0}{}_{\alpha \beta }$ is the viscous heating rate in the fluid rest frame (see Section 2.6 in Fujibayashi et al. 2018).

3.2.3. Properties of the NS

The top and middle panels of Figure 3 show the baryon mass and angular momentum of the remnant massive NSs (M_NS, J_NS) for models DD2-125M, SFHo-125H, and DD2-135M. M_NS and J_NS increase with time in an early stage of its evolution, t ≲ 1 s, in response to mass accretion. This mass accretion is caused primarily by the cooling of the disk material: the material located near the NS surface does not have Keplerian angular momentum and is sustained predominantly by the thermal pressure (see Figure 2). Through the neutrino cooling, the thermal pressure decreases and then the material falls onto the NS. The increase of the baryon mass is smaller for the SFHo case than for the DD2 cases simply due to the smaller disk mass for the SFHo case (see Table 3). For t ≳ 0.5 s, the total angular momentum of the NS begins to decrease, while the increase of the NS mass saturates. This is a consequence of the fact that neutrinos emitted from the rotating NS carry its angular momentum away (see Appendix A). The angular momentum loss rate is ∼10⁴⁸ erg s s⁻¹, which is consistent with an analytic estimation presented in Appendix A.

Zoom In Zoom Out Reset image size

The bottom panel of Figure 3 shows the neutrino luminosity of the massive NS and the disk for each model. The luminosity of the NS is ∼(5–7) × 10⁵² erg s⁻¹ for t ≲ 0.1 s and gradually decreases for t > 0.1 s. The slight increase of the luminosity for t ≳ 0.5 s is partly due to the finite grid resolution (see Appendix B). The neutrino luminosity of the disk is (1–3) × 10⁵³ erg s⁻¹ for t ≲ 0.1 s for each model, which is several times higher than that of the NS. The high luminosity is sustained for t ≲ 0.1 s by the viscous heating in the disk. However, as the temperature of the disk decreases due to the neutrino cooling and viscous expansion, the luminosity decreases for t ≳ 0.1 s. By comparing models DD2-135M and DD2-135M-v0, we find that the viscous effect enhances the neutrino luminosity of the NS and disk by a factor of 3–5 for t ≲ 0.1 s (for our choice of the viscous parameter, α_vis = 0.04). This enhancement of the neutrino luminosity results in the larger amount of the neutrino-driven ejecta (see Section 3.3.2).

3.2.4. Properties of the Disk

The quantities relevant to the disks for models DD2-125M, SFHo-125H, and DD2-135M are shown in Figure 4. Note that this "disk" component includes unbound material (see the definition in Section 3.2.2). The baryon mass and angular momentum of the disk (M_disk and J_disk) rapidly drop during an early time of t ≲ 1 s due to the mass accretion onto the NS caused by the neutrino cooling (see Section 3.2.3). The timescale of this accretion (i.e., the neutrino cooling timescale) is estimated by

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{M}_{{\rm{g}},\mathrm{NS}}{M}_{\mathrm{disk}}}{{R}_{\mathrm{NS}}{L}_{\nu ,\mathrm{disk}}}\approx 0.64\,{\rm{s}}\\ \ \ \times \ \left(\displaystyle \frac{{M}_{{\rm{g}},\mathrm{NS}}}{2.5{M}_{\odot }}\right)\left(\displaystyle \frac{{M}_{\mathrm{disk}}}{0.1{M}_{\odot }}\right){\left(\displaystyle \frac{{R}_{\mathrm{NS}}}{10\mathrm{km}}\right)}^{-1}{\left(\displaystyle \frac{{L}_{\nu ,\mathrm{disk}}}{{10}^{53}\mathrm{erg}{{\rm{s}}}^{-1}}\right)}^{-1},\end{array}\end{eqnarray} \tag{ 23 }$$

where M_g,NS is the gravitational mass of the NS. The numerical result is consistent with this timescale.

Zoom In Zoom Out Reset image size

The baryon mass and angular momentum of the disk for the SFHo case are smaller than those for the DD2 cases because of the smaller tidal effect on the merging NSs for the former. The disk mass is slightly larger for model DD2-125M than for DD2-135M, because the former has a smaller compactness and thus the self-gravity of the central massive NS is weaker. We note that the disk component has a nonnegligible fraction of the angular momentum of the remnant. The fraction of the disk angular momentum is 25%–30% of the total angular momentum of the remnant depending weakly on the EOS and the total mass of the system. Thus, only less than 75% of the total angular momentum of binary NSs at merger goes into the remnant massive NSs.

The disk mass is often defined by a certain density criterion (e.g., Shibata et al. 2017a; Radice et al. 2018). In the top-left panel of Figure 4, we also show the mass of the material with ρ < 10¹² g cm⁻³, which agrees approximately with the mass using our criterion, r > R_NS, at late times t ≳ 0.1 s.

We define the average radius of the disk R_disk by

$$\begin{eqnarray}&&{R}_{\mathrm{disk}}=\displaystyle \frac{1}{{M}_{{\rm{g}},\mathrm{NS}}}\displaystyle \frac{{J}_{\mathrm{disk}}^{2}}{{M}_{\mathrm{disk}}^{2}},\end{eqnarray} \tag{ 24 }$$

where we assumed that J_disk can be approximated by ${M}_{\mathrm{disk}}\sqrt{{M}_{{\rm{g}},\mathrm{NS}}{R}_{\mathrm{disk}}}$, and for simplicity, we use the Arnowitt–Deser–Misner mass (ADM mass; Arnowitt et al. 1960) obtained initially in the simulation as the gravitational mass of the NS. We note that R_disk depends only on the average specific angular momentum, J_disk/M_disk, and does not change by the heating or cooling of the disk. Moreover, this is a good indicator of the disk radius only for the case where the disk has Keplerian motion and does not appropriately indicate the disk radius when the material in the disk starts to be ejected from the system.

The middle-left panel of Figure 4 shows that R_disk gradually increases, indicating that the average specific angular momentum of the disk increases. For the earlier phase t ≲ 1 s, this occurs due to the accretion of the disk material with small specific angular momentum onto the NS due to the neutrino cooling. After the disk mass saturates for t ≳ 1 s, it is due to the viscous angular momentum transport from the outer envelope of the NS to the disk.

In the middle-right panel of Figure 4, we plot the ratios of the total (net) neutrino luminosity of the disk, L_ν,disk, which also includes the heating by neutrino interaction, to the viscous heating rate in the disk, L_vis,disk, for each model. The ratio increases with time for t ≲ 0.1 s because the optical depth of the disk, which is ∼10–100 at the beginning of the simulation, decreases due to the accretion and expansion of the disk. It peaks at t ∼ 0.1 s, when the optical depth of the disk decreases to ∼1 and neutrinos can carry away the thermal energy generated by the viscous heating and stored in the opaque disk. Then, the ratio decreases to ∼1 for 0.1 s ≲ t ≲ 1 s, when the disk is optically thin and the viscously generated thermal energy can be emitted immediately. For t ≲ 1 s, the ratio is close to or larger than unity, indicating that viscous heating does not contribute substantially to the disk expansion. However, after the accretion timescale of the disk (∼1 s; see Equation (23)), the neutrino emission and the disk accretion become weaker. Then, for t ≳ 1 s, this ratio drops steeply. This implies that the neutrino cooling is no longer efficient and most of the viscous heating is used for the disk expansion. We note that t ∼ 1 s also corresponds to the time at which the weak-interaction timescale becomes longer than the timescale of the viscous expansion of the disk (i.e., the weak interaction freezes out; see Section 3.3).

In Figure 4, we also plot the evolution of the average entropy (bottom left) and electron fraction (bottom right) of the disk, which increase gradually with time. The increase of the entropy is due to the viscous and neutrino heating. The electron fraction of the disk is determined by an equilibrium condition of electron/positron capture. Its equilibrium value becomes higher during the disk expansion because the electron degeneracy becomes weaker (the neutrino irradiation plays a minor role for the increase of the electron fraction; see Section 3.3). Its increase stagnates at t ∼ 0.1 s because of the feedback effect of the efficient neutrino emission (Chen & Beloborodov 2007; Kawanaka & Mineshige 2007): as the middle-right panel of Figure 4 shows, the efficiency of the neutrino emission at t ∼ 0.1 s is quite high. In such a situation, the increase of the specific entropy through viscous heating is suppressed (bottom-left panel of Figure 4). As a result, the electron degeneracy is kept strong, suppressing the increase of Y_e (bottom-right panel of Figure 4). Due to the neutrino irradiation, the value of the electron fraction at t ∼ 0.1 s in our models (∼0.2) is higher by ∼0.1 than that found in the context of the disk around a BH (Siegel & Metzger 2018). In the absence of the neutrino irradiation, the electron fraction becomes lower (see Section 3.2.5). After the neutrino cooling becomes inefficient, the entropy increases continuously, while the electron fraction saturates after the gradual increase until t ∼ 1 s, because the timescale of electron/positron capture in the disk becomes longer than the viscous evolution timescale due to the decrease of the temperature (see Section 3.3). The final average value of Y_e is 0.30–0.35.

3.2.5. Effect of the Neutrino Irradiation on the Properties of the Disk

3.3.1. Diagnostics

In this subsection, we describe the properties of the ejecta. First of all, we summarize the method to calculate the ejecta quantities. We define the mass of the ejecta that escaped from a sphere, M_ej,esc, by integrating the unbound material outflowing through the surface of a radius r_esc,

$$\begin{eqnarray}&&{M}_{\mathrm{ej},\mathrm{esc}}={\int }^{t}{dt}\displaystyle \frac{{{dM}}_{\mathrm{ej},\mathrm{esc}}}{{dt}},\end{eqnarray} \tag{ 25 }$$

where the mass-ejection rate is defined by

$$\begin{eqnarray}&&\displaystyle \frac{{{dM}}_{\mathrm{ej},\mathrm{esc}}}{{dt}}=\int {{dS}}_{k}{v}^{k}{\rho }_{* }H(-{{hu}}_{t}-{h}_{\min }),\end{eqnarray} \tag{ 26 }$$

with the area element on the sphere dS_k and the Heaviside function H. Hereafter, as in Section 2.7, r_esc is set to 8000 km, at which the temperature of the material is sufficiently low and the recombination of free nucleons has already occurred. Thus, Equation (12) can be a suitable condition for the material being unbound.

We also define the rate of the total energy (sum of the rest mass, and internal and kinetic energy) of the ejecta passing through the radius r_esc as

$$\begin{eqnarray}&&\displaystyle \frac{{{dE}}_{\mathrm{ej},\mathrm{esc}}}{{dt}}=\int {{dS}}_{k}{v}^{k}{\rho }_{* }\hat{e}H(-{{hu}}_{t}-{h}_{\min }),\end{eqnarray} \tag{ 27 }$$

where $\hat{e}=h\alpha {u}^{t}-P/(\rho \alpha {u}^{t})$ and the total energy of the ejecta is calculated by

$$\begin{eqnarray}&&{E}_{\mathrm{ej},\mathrm{esc}}={\int }^{t}{dt}\displaystyle \frac{{{dE}}_{\mathrm{ej},\mathrm{esc}}}{{dt}}.\end{eqnarray} \tag{ 28 }$$

This quantity is decomposed in the region far from the source of the gravity into

$$\begin{eqnarray}&&{E}_{\mathrm{ej},\mathrm{esc}}\approx {M}_{\mathrm{ej},\mathrm{esc}}+{T}_{\mathrm{kin}}+U+\displaystyle \frac{{M}_{\mathrm{ej},\mathrm{esc}}{M}_{{\rm{g}},\mathrm{NS}}}{{r}_{\mathrm{esc}}},\end{eqnarray} \tag{ 29 }$$

where T_kin and U are the kinetic and internal energies of the ejecta at a sufficiently distant region, and the last term approximately denotes the contribution from the gravitational binding energy. For the gravitational mass of the NS, we again use the ADM mass with the value obtained at the beginning of the simulation. From these expressions, we define the average ejecta velocity by

$$\begin{eqnarray}&&\langle {V}_{\mathrm{ej}}\rangle =\sqrt{2\left(\displaystyle \frac{{E}_{\mathrm{ej},\mathrm{esc}}}{{M}_{\mathrm{ej},\mathrm{esc}}}-1-\displaystyle \frac{{M}_{{\rm{g}},\mathrm{NS}}}{{r}_{\mathrm{esc}}}\right)},\end{eqnarray} \tag{ 30 }$$

where we assumed that the internal energy of the ejecta is eventually transformed into its kinetic energy during its expansion, i.e., ${T}_{\mathrm{kin}}+U={M}_{\mathrm{ej},\mathrm{esc}}\langle {V}_{\mathrm{ej}}{\rangle }^{2}/2$ at infinity. Because there are several ejecta components having different velocities, we define the velocity of the material passing through the sphere at a given time by

$$\begin{eqnarray}&&{V}_{\mathrm{ej}}=\sqrt{2\left({\left(\displaystyle \frac{{{dM}}_{\mathrm{ej},\mathrm{esc}}}{{dt}}\right)}^{-1}\displaystyle \frac{{{dE}}_{\mathrm{ej},\mathrm{esc}}}{{dt}}-1-\displaystyle \frac{{M}_{{\rm{g}},\mathrm{NS}}}{{r}_{\mathrm{esc}}}\right)},\end{eqnarray} \tag{ 31 }$$

where we used a differential form of Equation (29) and assumed

$$\begin{eqnarray}&&\displaystyle \frac{d}{{dt}}({T}_{\mathrm{kin}}+U)=\displaystyle \frac{1}{2}\displaystyle \frac{{{dM}}_{\mathrm{ej},\mathrm{esc}}}{{dt}}{V}_{\mathrm{ej}}{}^{2}.\end{eqnarray} \tag{ 32 }$$

With this definition, we can derive the velocity of the material ejected at different times.

There is still material satisfying Equation (12) inside the sphere with r = r_esc at the end of the simulations. Thus, we define the total mass by

$$\begin{eqnarray}\begin{array}{rcl}{M}_{\mathrm{ej},\mathrm{tot}} & = & {M}_{\mathrm{ej},\mathrm{esc}}\\ & & +\ {\displaystyle \int }_{500\mathrm{km}\leqslant r\lt {r}_{\mathrm{esc}}}{d}^{3}x{\rho }_{* }H(-{{hu}}_{t}-{h}_{\min }).\end{array}\end{eqnarray} \tag{ 33 }$$

This value is interpreted approximately as the maximum mass that the merger remnant can eject. We define the mass-weighted average electron fraction of this material by

$$\begin{eqnarray}\begin{array}{rcl}\langle {Y}_{{\rm{e}},\mathrm{ej}}\rangle & = & \displaystyle \frac{1}{{M}_{\mathrm{ej},\mathrm{tot}}}\left({\displaystyle \int }^{t}{dt}\displaystyle \int {{dS}}_{k}{v}^{k}{\rho }_{* }{Y}_{{\rm{e}}}H(-{{hu}}_{t}-{h}_{\min })\right.\\ & & +\ \left.{\displaystyle \int }_{500\mathrm{km}\leqslant r\lt {r}_{\mathrm{esc}}}{d}^{3}x{\rho }_{* }{Y}_{{\rm{e}}}H(-{{hu}}_{t}-{h}_{\min }\right),\end{array}\end{eqnarray} \tag{ 34 }$$

which includes both the ejecta that already escaped from and still stays inside the sphere of r = r_esc. In order to derive the electron fraction of the different ejecta components, we also define the flux-weighted average electron fraction of the material passing through the sphere at a given time by

$$\begin{eqnarray}&&{Y}_{{\rm{e}},\mathrm{ej}}=\displaystyle \frac{1}{{{dM}}_{\mathrm{ej},\mathrm{esc}}/{dt}}\int {{dS}}_{k}{v}^{k}{\rho }_{* }{Y}_{{\rm{e}}}H(-{{hu}}_{t}-{h}_{\min }).\end{eqnarray} \tag{ 35 }$$

3.3.2. Contribution of Different Ejecta Components

In Figure 5, the ejecta mass (top left), mass-ejection rate (bottom left), ejecta velocity (top right), and average electron fraction (bottom left) for models DD2-135M, DD2-135M-irr, and DD2-135M-v0 are compared to explore the contribution of the effects of neutrino irradiation and viscosity. Here, we focus on the models of the remnant of the merger with a large total gravitational mass of 2.7 M_⊙. The following discussion is also provided for models DD2-125M and SFHo-125M. At the very beginning of the simulation, unbound material of M_ej,tot ≈ 0.004 M_⊙ was present, reflecting the merger process for this model. This amount is somewhat larger than that shown in Table 2 because of the difference in the definition of ejecta (see Sections 2.3 and 2.5). For the inviscid model DD2-135M-v0, M_ej,tot peaks at t ∼ 0.1 s and then decreases, indicating that part of the ejecta was not actually the ejecta component, and it falls back to the central region. Its saturated value for t ≳ 1 s, ≈0.003M_⊙, is the actual mass of the ejecta. With the viscosity, on the other hand, a larger amount of material is ejected as described below.

Zoom In Zoom Out Reset image size

For model DD2-135M, it is found that M_ej,tot increases at three different times, t ∼ 0.01, ∼0.1, and ∼1 s (see the top right panel of Figure 5). The first increase of M_ej,tot at t ∼ 0.01 s is also found for model DD2-135M-irr but absent for model DD2-135M-v0, and therefore, this increase is due to the viscous effect. This is the contribution of the early viscosity-driven ejecta. By comparing models DD2-135M and DD2-135M-v0, the contribution of this component for our choice of the viscous parameter (α_vis = 0.04) is found to be ∼0.01 M_⊙. The second increase at t ∼ 0.1 s is absent for model DD2-135M-irr, and hence, we confirm that this ejecta is driven by neutrino irradiation. By comparing models DD2-135M and DD2-135M-irr, the contribution of this neutrino-driven component is found to be ∼0.01 M_⊙. The contribution of neutrino irradiation is fairly large compared to DD2-135M-v0 because of the large neutrino luminosity of the remnant due to the viscous heating (see Figure 3). The third one at t ∼ 1 s, which is found in both models DD2-135M and DD2-135M-irr, is the contribution of the late-time viscosity-driven ejecta.

In the bottom-left panel of Figure 5, only two (not three) distinct mass-ejection phases are found, one of which is found in an early phase of t ≲ 0.5 s and the other is in a late phase of t ≳ 1 s. The former is composed of the dynamical, early viscosity-driven, and neutrino-driven components. The value of dM_ej,esc/dt for DD2-135M starts to deviate from that for model DD2-135M-irr at t ≈ 0.08 s, indicating that the neutrino-driven component starts to contribute to M_ej,esc at this time. At t ∼ 0.1 s, the properties of the ejecta are determined by both the dynamical and neutrino-driven ejecta, the former of which has a lower electron fraction of ≈0.2 and the latter a higher value of ≈0.55 (see Figure 8). Y_e,ej (the bottom-right panel of Figure 5) then has the average value of the two components, which is ≈0.40–0.45.

After that, the early viscosity-driven ejecta begins to contribute predominantly to the properties of the ejecta at t ≈ 0.2 s, at which a peak in dM_ej,esc/dt is remarkable in particular for model DD2-135M-irr. This component has a low electron fraction of ≈0.25 as found for model DD2-135M-irr. Y_e,ej also decreases to ≈0.3 at this time for model DD2-135M as a result of the large contribution of the early viscosity-driven ejecta to dM_ej,esc/dt.

In the top-right panel of Figure 5, a peak of V_ej at t ∼ 0.1 s is found for model DD2-135M but absent for DD2-135M-irr, indicating that the neutrino-driven component has a velocity of >0.3 c, which is larger than that of the dynamical and early viscosity-driven components at t ≈ 0.2 s. This results in the simultaneous contribution of these ejecta components to dM_ej,esc/dt at t ∼ 0.1 s, although they are ejected at different times as found in the top-left panel of Figure 5.

Finally, the late-time viscosity-driven ejecta associated with the viscous effect in the disk sets in. dM_ej,esc/dt has a large value of ≳0.01 M_⊙ s⁻¹ for t = 1–3 s, and after that, it decreases with time. This component has a smaller velocity of V_ej ∼ 0.05–0.10c and a range of the average electron fraction of Y_e,ej ≈ 0.30–0.35, reflecting the electron fraction of the disk (Y_e,disk; see Figure 4).

For model DD2-135M, the mass-weighted average electron fraction over all ejecta components is $\langle {Y}_{{\rm{e}},\mathrm{ej}}\rangle \approx 0.34$ at the end of the simulation due to the domination of the late-time ejecta (see Table 3). By comparing models DD2-135M and DD2-135M-irr in Table 3, $\langle {Y}_{{\rm{e}},\mathrm{ej}}\rangle$ is lower by 0.05 in the absence of neutrino irradiation (reflected by the lower electron fraction of the disk; see Section 3.2.5).

It is worth mentioning that even in the absence of neutrino irradiation, the electron fraction of the ejecta increases to ≈0.3 due to the positron capture after the decrease of the electron degeneracy in the disk (see Section 3.3.4). This indicates that, even in the absence of a long-lived massive NS, the electron fraction of the ejecta would not be very low as pointed out in recent work for BH–disk systems (Fernández et al. 2020; Fujibayashi et al. 2020).

3.3.3. Ejecta Properties for Fiducial Models

Figure 6 compares the ejecta properties for models DD2-125M, SFHo-125H, and DD2-135M. The curves are qualitatively similar to each other, showing that there are three mass-ejection mechanisms as described in Section 3.3.2 irrespective of the adopted EOS and the total mass of the model. The baryon mass of the unbound material for t ≲ 0.03 s, composed of dynamical and early viscosity-driven ejecta, depends on models. The baryon mass of these early ejecta components correlates with the baryon mass of the dynamical ejecta (see Table 2). The amount of the subsequent neutrino-driven and late-time viscosity-driven ejecta, on the other hand, correlates with the disk mass (see the top-left panel of Figure 4). The value of M_ej,tot saturates at t = 2–3 s, and by comparing with the top-left panel of the Figure 4, the saturated values indicate that approximately 30% of the initial disk mass becomes the ejecta.

Zoom In Zoom Out Reset image size

For all fiducial models, the average ejecta velocity is $\langle {V}_{\mathrm{ej}}\rangle \sim 0.1\,c$ (see Table 3), which is somewhat larger than that of the late-time ejecta (see the top-left panel of Figure 6) because of the large kinetic energy of the dynamical, early viscosity-driven, and neutrino-driven ejecta components.

Figure 7 shows the mass histograms of the electron fraction (top), entropy (middle), and expansion timescale (bottom; t_exp = r/v^r ) of the tracer particles at T = 5 GK for the fiducial models. The distributions for all models are quite similar and peak at Y_e ≈ 0.3, s/k_B ≈ 10, and t_exp ≈ 0. 1 s. These peaks for each model are produced by the late-time viscosity-driven ejecta. In this figure, the histogram in the corresponding 3D simulations is also shown. This illustrates that the low electron fraction and short expansion timescale sides of the distributions are determined by the dynamical ejecta. ¹⁰

Zoom In Zoom Out Reset image size

The temporal Y_e distributions for the tracer particles that reach r_esc by a given time are shown in Figure 8 for models DD2-125M, SFHo-125H, and DD2-135M. At t = 0.2 s, there are two peaks at Y_e ∼ 0.2 and 0.55, which are determined by the dynamical plus early viscosity-driven ejecta and neutrino-driven ejecta, respectively. The subsequent change of the histograms clearly illustrates that the peak of Y_e in the final distribution is developed by the late-time viscosity-driven ejecta for t > 1 s.

Zoom In Zoom Out Reset image size

The distributions of the ejecta in the s–Y_e plane are shown in Figure 9 for models DD2-125M, SFHo-125H, and DD2-135M. It is found that the distribution profiles depend very weakly on the EOSs and binary total mass. We also find a positive correlation between s and Y_e. The ejecta with Y_e ≈ 0.3 has an entropy of s/k_B ≲ 10, while the ejecta with higher values of Y_e has a higher entropy. The electron fraction for most of the material ejected toward the nonpolar direction is determined predominantly by the equilibrium between electron/positron capture (see Section 3.3.4), and the equilibrium value becomes higher for weaker electron degeneracy. Because the entropy negatively correlates with the electron degeneracy, the material with the higher entropy has the higher equilibrium value of the electron fraction. This results in the positive correlation between s and Y_e. In other words, the material experiencing the longer-time viscous heating has the higher electron fraction reflecting its higher entropy. By comparing the bottom two panels of Figure 9, it is found that the low-Y_e component is more abundant in the absence of neutrino irradiation.

Zoom In Zoom Out Reset image size

In Figure 9, there is a component that has a low electron fraction of Y_e = 0.1–0.2 and a moderate entropy of s/k_B ∼ 50. This component is a part of the dynamical ejecta, which is initially expelled by the tidal effect and then heated up by a shock formed in spiral arms. Because this component experiences shock heating in a distant region from the center, the increase in the temperature due to the shock is not enough for the weak interaction (positron capture) to work in a sufficiently short timescale. Thus, the electron fraction of this component remains low.

3.3.4. Processes that Determine the Electron Fraction of the Ejecta

The electron fraction of the material is approximately determined by the absorption of electron-type (anti)neutrinos and the electron/positron capture on free nucleons. The time evolution of the electron fraction is described by

$$\begin{eqnarray}&&{\dot{Y}}_{{\rm{e}}}={\dot{Y}}_{{\rm{e}},\nu }+{\dot{Y}}_{{\rm{e}},\mathrm{pc}}-{\dot{Y}}_{{\rm{e}},\bar{\nu }}-{\dot{Y}}_{{\rm{e}},\mathrm{ec}},\end{eqnarray} \tag{ 36 }$$

where ${\dot{Y}}_{{\rm{e}},\nu }$, ${\dot{Y}}_{{\rm{e}},\bar{\nu }}$, ${\dot{Y}}_{{\rm{e}},\mathrm{ec}}$, and ${\dot{Y}}_{{\rm{e}},\mathrm{pc}}$ are the change rates of the electron fraction due to the absorption (or capture) of electron neutrinos, electron antineutrinos, electrons, and positrons, respectively, on free nucleons. These reaction rates depend on the local density, temperature, electron fraction, and fluxes (and energies) of electron neutrinos and antineutrinos. In the top panel of Figure 10, the variations of the timescales of these reactions with decreasing temperature in tracer particles are shown together with their expansion timescales defined by t_exp = r/v^r with the coordinate radius r and radial velocity v^r (=dr/dt). Here we define the timescale of the electron/positron capture by ${t}_{\mathrm{cap}}=\min ({Y}_{{\rm{e}}}/{\dot{Y}}_{{\rm{e}},\mathrm{ec}},{Y}_{{\rm{e}}}/{\dot{Y}}_{{\rm{e}},\mathrm{pc}})$ and that of the neutrino absorption by ${t}_{\nu }=\min ({Y}_{{\rm{e}}}/{\dot{Y}}_{{\rm{e}},\nu },{Y}_{{\rm{e}}}/{\dot{Y}}_{{\rm{e}},\bar{\nu }})$. It is found that until the timescales of weak-interaction processes become longer than the expansion timescale, the electron fraction of the tracer particle is determined mainly by the electron/positron capture (red) and subdominantly by the neutrino absorption (blue).

Zoom In Zoom Out Reset image size

In the middle panel of Figure 10, the variation of the electron fraction with decreasing temperature in the tracer particles is shown together with its equilibrium values for neutrino and antineutrino absorption, Y_e,ν , and electron and positron captures, Y_e,cap. These equilibrium values are determined by the equations ${\dot{Y}}_{{\rm{e}},\nu }={\dot{Y}}_{{\rm{e}},\bar{\nu }}$ and ${\dot{Y}}_{{\rm{e}},\mathrm{ec}}={\dot{Y}}_{{\rm{e}},\mathrm{pc}}$. Here, we assumed that the material is not optically thick to neutrinos. That is, in our analysis, we pay attention only to the temperature range in which the beta equilibrium is not established. This assumption is reasonable for low temperatures of k_B T ≲ 3 MeV. It is found that the electron fraction of the tracer particle is determined approximately by the equilibrium value for electron and positron captures for k_B T ∼ 1–3 MeV. We note that the electron fraction deviates from Y_e,cap for k_B T ≳ 3 MeV because the material is opaque to neutrinos and the electron fraction should be determined by the beta equilibrium; the electron fraction is determined not only by the electron/positron captures but also by electron (anti)neutrino absorption. By contrast, at k_B T ∼ 1 MeV, the electron fraction freezes out because the weak-interaction timescale becomes longer than the expansion timescale (see the top panel of Figure 10). This illustrates how the electron fraction is determined in the material ejected toward the nonpolar direction (see Figure 8). The effect of the neutrino absorption on the electron fraction is minor for most of the material ejected toward the nonpolar direction. On the other hand, for the material ejected toward the polar direction, the electron fraction is determined predominantly by the neutrino absorption, and it is higher than that of the ejecta toward the nonpolar direction (see Figure 1).

Because the equilibrium value of the electron fraction changes steeply in the temperature region of k_B T ≈ 1–2 MeV, at which the freeze out of weak interaction occurs, the electron fraction of the ejecta should be affected significantly by the expansion timescale. In the presence of strong viscous effect with α_vis > 0.04 or the other short-timescale mass-ejection processes (for example, the Lorentz force by aligned magnetic fields; see Section 3.5), the electron fraction of the ejecta would become lower (see Section 3.4).

In the bottom panel of Figure 10, the specific heating rates due to viscous dissipation and neutrino absorption are shown together with the specific cooling rate due to the neutrino emission. We note that the cooling rate shown in this figure is equal to the leakage rate in Equation (4), in which neutrino trapping is taken into account. The heating rate is suppressed in the optically thick region in our leakage-based neutrino radiation transfer method because part of the neutrino absorption is assumed to balance out the neutrino emission (see Fujibayashi et al. 2017 for a detailed description of our method). As a result, the heating rate for k_B T ≳ 2 MeV decreases with increasing temperature (i.e., high-density and high-optical depth region). This figure illustrates that viscous heating approximately balances out neutrino cooling for k_B T ≳ 1 MeV (i.e., for the early phase of the evolution), while viscous heating dominates over neutrino cooling for k_B T ≲ 1 MeV (i.e., for the later phase after the neutrino cooling efficiency drops). The late-time viscosity-driven mass ejection occurs in such conditions.

To explore the dependence of the results on the viscous coefficient in particular for the case of extremely high viscosity, we perform a simulation of model DD2-125M-h (α_vis = 0.1). In the top-left panel of Figure 11, the ratio of neutrino luminosity to viscous heating rate for model DD2-125M-h is compared to that for model DD2-125M. Due to the faster viscous expansion of the disk, the neutrino cooling efficiency drops and the mass ejection sets in earlier than for DD2-125M. In addition, due to the larger viscous effect, the ejecta mass becomes larger: the total ejecta mass for this model is M_ej,tot ≈ 0.2 M_⊙ at the end of the simulation, which is approximately twice as large as that for DD2-125M. The reason for this enhancement is that the viscous timescale is shorter, and hence, a substantial fraction of material is accelerated outward by the neutrino and viscous heating before the material accretes onto the NS.

Zoom In Zoom Out Reset image size

Due to the earlier mass ejection from the disk, the equilibrium value of the electron fraction at its freeze out (t ≈ 1 s for DD2-125M-h and t ≈ 2 s for DD2-125M) becomes lower (Fujibayashi et al. 2020), and thus, the electron fraction of the ejecta becomes lower as well (see the bottom panels of Figure 11). The bottom-right panel of Figure 11 shows that the mass fraction with low electron fraction (Y_e < 0.25) for model DD2-125M-h is remarkably larger than for model DD2-125M.

Figure 12 shows the mass distribution of the ejecta in the s–Y_e plane for model DD2-125M-h. It is clearly found that the ejecta with such low electron fraction is present in the low-entropy (s/k_B ∼ 10) region. This illustrates that the electron fraction of the ejecta is sensitive to the onset time of the mass ejection. As we discuss in Section 3.5, the mass ejection may proceed in a shorter timescale in the presence of strong magnetic fields. The ejecta mass may be enhanced and the electron fraction of the ejecta may be decreased in the presence of such an efficient mass-ejection process.

Zoom In Zoom Out Reset image size

The remnant NS and disk formed after the NS mergers would be strongly magnetized due to magnetohydrodynamical instabilities like the Kelvin–Helmholtz instability (Price & Rosswog 2006) and the MRI (Balbus & Hawley 1991), and thus the magnetic field effect could play an important role for determining the ejecta dynamics and the activity of the remnant. In this subsection, we discuss the possible effects of the magnetic field in its presence.

Figure 1 shows that the density in the polar region decreases with time. This indicates that the magnetic pressure by the strong magnetic field would modify the fluid dynamics in such a region. Figure 13 shows the pressure profile for model DD2-125M along the polar direction for t = 0.01–6 s together with the hypothetical magnetic pressure P_B  = B²/8π assuming that the magnetic field has a dipole structure with the magnetic strength at the polar surface of the NS as B_p = 10¹⁵ G. This figure suggests that the pressure along the polar direction substantially decreases after ≈0.3 s and that the magnetar-level magnetic field could determine the fluid dynamics along the pole: nearly force-free magnetosphere is likely to be produced. Using the result of a force-free electromagnetic simulation by Spitkovsky (2006), the electromagnetic luminosity is

$$\begin{eqnarray}\begin{array}{rcl}L & \approx & \displaystyle \frac{1}{4}{B}_{{\rm{p}}}^{2}{R}^{6}{{\rm{\Omega }}}^{4}=6.6\,\times \,{10}^{49}\,\mathrm{erg}\,{{\rm{s}}}^{-1}\\ & & \times \ {\left(\displaystyle \frac{{B}_{{\rm{p}}}}{{10}^{15}{\rm{G}}}\right)}^{2}{\left(\displaystyle \frac{{R}_{\mathrm{NS}}}{15\mathrm{km}}\right)}^{6}{\left(\displaystyle \frac{{\rm{\Omega }}}{5000\mathrm{rad}{{\rm{s}}}^{-1}}\right)}^{4}.\end{array}\end{eqnarray} \tag{ 37 }$$

Because the force-free condition would be satisfied only in the polar region, part of this value would be released toward the polar direction. If ≳10% of L is released, the long-lived and highly magnetized remnant NS, which should be rapidly rotating, could drive a high-Lorentz factor outflow and be an engine of gamma-ray bursts because the baryon-loading problem is likely to be naturally avoided, and a funnel structure would be helpful for confining the relativistic jet.

Zoom In Zoom Out Reset image size

Because of the differential rotation of the remnant system, the toroidal magnetic field would be amplified by the winding effect. It would result in the increase in magnetic pressure and the development of a tower-like outflow along the rotational axis (Meier 1999). The electromagnetic luminosity due to this mechanism would be

$$\begin{eqnarray}&&\begin{array}{l}L\sim {10}^{50}\,\mathrm{erg}\,{{\rm{s}}}^{-1}\\ \qquad \ \times \ {\left(\displaystyle \frac{B}{{10}^{15}{\rm{G}}}\right)}^{2}{\left(\displaystyle \frac{{R}_{\mathrm{NS}}}{10\mathrm{km}}\right)}^{3}\left(\displaystyle \frac{{\rm{\Omega }}}{{10}^{4}\,\mathrm{rad}\,{{\rm{s}}}^{-1}}\right){\left(\displaystyle \frac{\delta }{0.1}\right)}^{3},\end{array}\end{eqnarray} \tag{ 38 }$$

where ${(\delta {R}_{\mathrm{NS}})}^{3}$ is an effective volume for which the magnetic field amplification occurs. At the same time, the material at the polar surface of the NS would be stripped by the strong magnetic field, developing outflow toward the polar direction (Shibata et al. 2011b).

The effective viscosity may not be the only magnetohydrodynamical effect on the disk evolution. If there is a sufficiently strong aligned poloidal magnetic field in the disk, the Lorentz force would expel the disk material, in particular the neutron-rich material located in the region far from the central region, on a shorter timescale than the viscosity. The material ejected by such an effect could be neutron rich if its ejection timescale is short enough for the electron fraction of the material to freeze out earlier (Siegel & Metzger 2018; Fernández et al. 2019; Fujibayashi et al. 2020).

Nucleosynthetic yields for each tracer particle (Section 2.7) are calculated in a postprocessing step by using a nuclear reaction network code, rNET, as described in Wanajo et al. (2018). The network consists of 6300 isotopes (Z = 1–110) that are connected by experimentally evaluated rates when available (JINA REACLIB V2.0, ¹¹ Cyburt et al. 2010; Nuclear Wallet Cards ¹² ) and those from theoretical models otherwise. Theoretical rates for neutron, proton, and alpha capture (TALYS; Goriely et al. 2008) and β-decay (GT2; Tachibana et al. 1990) are based on a microscopic nuclear-mass prediction (HFB-21; Goriely et al. 2010). Neutrino-induced reactions are not included in the nucleosynthesis calculations, because they are expected to play only minor roles in our present models (except for setting the values of Y_e at ≳10 GK; see the middle panel of Figure 10).

Each nucleosynthesis calculation starts when the temperature decreases to 10 GK with the initial composition of 1 − Y_e and Y_e for free neutrons and protons, respectively. At such high temperature, NSE is immediately established. Although the temperatures in the viscosity-driven ejecta exceed 10 GK, those in the early dynamical ejecta (∼1% of the total amount) already have decreased below 10 GK at ≈50 ms (for the 2D cases). For a tracer particle in which the temperature does not reach 10 GK, therefore, the nucleosynthesis calculation starts from the beginning of each postmerger evolution, corresponding to ≈50 ms after the merger. The initial composition for such cases is adopted from the nucleosynthesis abundances at 50 ms in the dynamical ejecta of Wanajo et al. (2014) with (almost) the same initial value of Y_e, because an NSE condition is not assured (in particular for those with <5 GK).

The calculated ejecta masses of nucleosynthesis yields for models DD2-125M, DD2-125M-h, SFHo-125H, and DD2-135M are compared in the left panel of Figure 14 (solid lines). Here, the result for DD2-125M-h is added as a limiting case with a large value of α_vis = 0.10. The r-process residuals to the mass spectrum of the solar system abundances (hereafter the solar r-residuals; Prantzos et al. 2020) are also plotted, which are scaled to fit the ejecta mass of ⁸²Se (one of the first peak abundances) in DD2-125M. Note that our 2D postmerger simulations contain the dynamical ejecta. In fact, the independent nucleosynthesis results (displayed by the dotted lines in the left panel of Figure 14) by using the tracers from corresponding 3D simulations (Table 2) are in reasonable agreement with those from 2D simulations (solid lines, except for DD2-125M-h) for the heavy r-process nuclei (A > 130). For DD2-125M, the masses of these nuclei in 2D (blue solid line) are about a factor of a few larger than those in 3D (blue dotted line), indicating that the marginally bound material in the dynamical phase is pushed out by the subsequent early viscosity-driven outflows (Section 3.1). It is interesting to note that the resultant abundance patterns (except for DD2-125M-h) are quite similar to those in the ejecta from BH accretion disks (expected to be formed after massive binary NS mergers) explored in Fujibayashi et al. (2020).

Zoom In Zoom Out Reset image size

Overall, the abundance patterns are similar among the models with the same viscous parameter (i.e., except for DD2-125M-h), being independent of the NS masses and EOSs adopted. While the abundances lighter than A = 110 are in good agreement with the solar r-residual pattern, the heavier nuclei are sizably underabundant (for α_vis = 0.04). This is a consequence of the fact that the relatively low entropies (∼10–20 k_B; Figure 7) and mild neutron richness (Y_e ∼ 0.25–0.5, Figure 8) in the bulk of postmerger ejecta give the nucleosynthesis-relevant conditions for either of NSE, nuclear quasi-statistical equilibrium, or only a weak r-process (e.g., Wanajo et al. 2011, 2018). Conversely, model DD2-125M-h results in the solar r-process-like abundance pattern because of its relatively neutron-rich ejecta (Y_e ∼ 0.1–0.35; bottom-right panel in Figure 11).

The outcomes for the models with our fiducial choice of α_vis = 0.04 (DD2-125M, SFHo-125H, and DD2-135M) conflict with the robustness of elemental abundance patterns over a wide range of atomic numbers (in particular for Z ≥ 56 and to a lesser extent for Z > 38) among all the r-process-enhanced (or "main" r-process) stars in the Galaxy (e.g., Cowan et al. 2019). In the right panel of Figure 14, the measured abundances of such a star, CS 31081–001 (crosses; Siqueira Mello et al. 2013), are compared with those of our models, along with the solar r-residual pattern (gray line; Prantzos et al. 2020), which are normalized with respect to Zr (Z = 40). We find that the calculated abundance patterns except for DD2-125M-h are rather consistent with those of HD 122563 (circles; Cowan et al. 2005; Honda et al. 2006; Roederer et al. 2012b), one of the r-process-deficient metal-poor stars showing a descending abundance trend and referred to as a "weak" r-process star (Wanajo & Ishimaru 2006). Thus, low-mass NS binaries may be the sources of such weak r-process-like signatures found in metal-poor stars. ¹³ Alternatively, the viscosity in the disk should be effectively as large as α_vis = 0.10, which is adopted in DD2-125M-h for reproducing the solar r-process abundance pattern.

As displayed in Figure 15, each model with α_vis = 0.04 exhibits a solar r-process-like abundance pattern over a wide range of A (∼80–200) only when the time elapsed for the postmerger phase is below 0.2 s. Even if a central remnant collapsed at this epoch (e.g., by the replacement of the binary mass with a larger one or the EOS with a softer one), the subsequent mass ejection from a BH accretion disk would add a substantial amount of material with A < 130 as found in Fujibayashi et al. (2020). We conclude, therefore, that the binary NS systems explored in this study, that is, those ejecting small (∼0.001 M_⊙) and large (>0.05 M_⊙) masses in the dynamical and postmerger phases, respectively, cannot be the predominant source of the Galactic r-process material, given that our choice of α_vis = 0.04 is representative. This implies that the frequency of the low-mass binary NS mergers, leading to long-lived massive NSs as remnants, would be rather low, given that the binary NS mergers were the predominant site of r-process nucleosynthesis. However, if the viscous effect is effectively large with α_vis ∼ 0.1, or in other words, in the presence of an efficient mass-ejection process in the early postmerger phase, this conclusion is significantly modified, as we already mentioned in the previous paragraph.

Zoom In Zoom Out Reset image size

We inspect the radioactive decays of the nuclei synthesized in the postmerger ejecta, which give rise to kilonova emission. For this purpose, the nuclear (specific) heating rates ($\dot{q}$; in units of erg g⁻¹ s⁻¹) are computed from the temporal evolution of nucleosynthesis for DD2-125M (blue), DD2-125M-h (cyan), SFHo-125H (green), and DD2-135M (red) as displayed in the left panel of Figure 16 (solid lines). The heating is chiefly due to β-decays; the contributions from fission and α-decays are unimportant except for DD2-125M-h because of the small amount of trans-lead species synthesized in the ejecta (Figure 14). The heating rates for all models resemble each other, in particular for DD2-125M, SFHo-125H, and DD2-135M, because of their similar nucleosynthesis outcomes. However, these deviate from the empirical power-law-like heating (with the index −1.3; dotted line) that has been found in previous studies as a result of the decaying second-peak nuclei (A ∼ 130) with various half-lives (Metzger et al. 2010; Wanajo et al. 2014).

Zoom In Zoom Out Reset image size

The reason for this deviation is evident from the right panel of Figure 16 that shows the heating rates from the β-decays of individual isotopes, ${\dot{q}}_{\beta }$, which have dominant contributions to the total heating rate, $\dot{q}$ (for DD2-125M as representative). The main contributors are those between the iron group and the first peak nuclei (A ∼ 50–90) that dominate in the nucleosynthesis yields (Figure 14). Only two second-peak nuclei (A = 131 and 132) marginally contribute to the heating rates. Most of these isotopes (listed in Table 4) have half-lives (second column) either of a few hours, a few days, or several tens of days. As a result, these nuclei contribute to heating at different durations, that is, for t < 1 day, t = 1–10 days, and t > 10 days, making three bump-like structures in the curve of $\dot{q}$ (black solid line). For instance, the contributions of the two decay chains ⁶⁶Ni ${\to }^{66}$Cu ${\to }^{66}$Zn and ⁷²Zn ${\to }^{72}$Ga ${\to }^{72}$Ge with similar half-lives (∼2 days; Table 4) stand out in the heating rate between 1 and 10 days (as also found in Wanajo 2018; see the discussion on the contribution of various isotopes in different astrophysical conditions in Metzger et al. 2010; Grossman et al. 2014; Lippuner & Roberts 2015; Martin et al. 2015). DD2-125M-h results in a larger heating rate than those in the other models because of the abundant second-peak elements in its ejecta coming into play.

Table 4. Dominant Decay Chains and Energy Partitions (MeV) ^a

Isotope		Half-life	Gamma-Ray	Electron	Neutrino
66Ni $\to$ 66Cu		2.27 d	0	0.0733	0.179
	$\to$ 66Zn	5.12 m	0.0978	1.07	1.48
72Zn $\to$ 72Ga		1.94 d	0.145	0.0805	0.194
	$\to$ 72Ge	14.1 h	2.77	0.471	0.759
77Ge $\to$ 77As		11.3 h	1.08	0.642	0.982
	$\to$ 77Se	1.62 d	0.00808	0.226	0.448
78Ge $\to$ 78As		1.47 h	0.278	0.227	0.450
	$\to$ 78Se	1.51 h	1.31	1.26	1.68
88Kr $\to$ 88Rb		2.84 h	1.95	0.364	0.606
	$\to$ 88Sr	17.8 m	0.677	2.05	2.58
89Sr $\to$ 89Y		50.5 d	0	0.585	0.912
91Sr $\to$ 91Y		9.63 h	0.707	0.641	0.992
	$\to$ 91Zr	58.5 d	0.00313	0.603	0.938
131I $\to$ 131Xe		8.03 d	0.380	0.182	0.396
132Te $\to$ 132I		3.20 d	0.212	0.0670	0.173
	$\to$ 132Xe	2.30 h	2.26	0.486	0.837

Note.

^a Data are taken from the ENDF/B-VII.1 library (Chadwick et al. 2011).

Download table as:  ASCIITypeset image

The total heating rates, $\dot{Q}$, in the ejecta for each model are calculated as the product of the specific heating rate, $\dot{q}$, the thermalization factor, f_th, and the mass of the ejecta, M_ej,tot. The thermalization factors for the γ-rays, electrons, α-particles, and fission fragments (the latter two are not shown) are obtained by using the analytic formula in Barnes et al. (2016) with the ejecta mass and the average velocity estimated at the end of the simulation for each model (10th and 11th columns in Table 3), as shown in the left panel of Figure 17. The factors (f_th) for γ-rays rapidly decay after several days, while those for electrons slowly decrease at later times (see also Hotokezaka & Nakar 2020). In general, f_th is larger for more massive or slowly expanding ejecta. In our models, the average velocities (∼0.1c) as well as the heating rates are similar to each other, and thus, the ejecta mass is the main source of differences in $\dot{Q}$.

Zoom In Zoom Out Reset image size

The computed curves of $\dot{Q}$ are compared with the bolometric luminosity of the kilonova (Kasliwal et al. 2017, 2019; Smartt et al. 2017; Waxman et al. 2018) associated with the NS merger GW170817. The right panel of Figure 17 displays the results for models DD2-125M (blue), DD2-125M-h (cyan), SFHo-125H (green), and DD2-135M (red). According to the inferred lower bound of the binary mass for GW170817 (2.73 M_⊙; Abbott et al. 2017b), the former three should be taken as predictions for future low-mass binary events, while the latter may be taken as representative for this event. In fact, we find that the curve of $\dot{Q}$ for DD2-135M is in good agreement with the observed light curve (different symbols with error bars) between 1 and 10 days, including the bump-like signature at several days after the merger. During this epoch, the radioactive heating is predominantly due to the decay chains from ⁶⁶Ni and ⁷²Zn, which also is suggested in Wanajo et al. (2018). ¹⁴ Despite the heating rate (Figure 16) being about a factor of 2 smaller than the empirical rate (dotted line) expected from a solar r-process-like abundance distribution (e.g., Metzger et al. 2010; Wanajo et al. 2014), the large M_ej,tot (=0.086 M_⊙) and the subsequent high f_th keep $\dot{Q}$ sufficiently high. Conversely, the ejecta masses for the other models are larger or smaller to match the light curve of the merger GW170817. The behaviors of $\dot{Q}$ are, however, similar to each other.

For early times (t < 1 days), the radioactive energy has been lost by the adiabatic expansion of the ejecta and thus the heating rate is larger than the bolometric luminosity of the kilonova. During this epoch, the two β-decay chains near the first peak (A ∼ 80) at the isobars A = 78 and 88 (Table 4 and the right panel of Figure 16) play dominant roles in brightening the kilonova emission. It is interesting to note that the identification of the Sr absorption line in the ejecta of the merger GW170817 (Watson et al. 2019) may also be indicative of the contribution of the β-decay chain at A = 88 to the luminosity of the kilonova near its peak.

For late times (t > 10 days), the main contributors to the kilonova emission are also the isotopes near the first peak (except for DD2-125M-h), as a result of the successive β-decay at the isobars A = 89 and 91 (Table 4 and the right panel of Figure 16). The relatively long half-lives of ⁸⁹Sr (50.5 days) and ⁹¹Y (58.5 days) can give rise to the long-lasting kilonova emission. At this epoch, the Spitzer 4.5 μm detection (the right panel of Figure 17; triangles) places a tight lower bound, although the bolometric luminosity itself is highly uncertain. The half-lives of ⁸⁹Sr and ⁹¹Y are similar to that of the spontaneous fission of ²⁵⁴Cf (60.5 days), which is suggested to be a dominant heating source of kilonovae at late times (Wanajo et al. 2014; Hotokezaka et al. 2016; Zhu et al. 2018; Wu et al. 2019). This would make it difficult to confirm the presence of the heaviest r-process nuclei, ²⁵⁴Cf, from the observed kilonova light curves at late times as suggested by Zhu et al. (2018). Note that the decay chain (successive electron capture) of ⁵⁶Ni ${\to }^{56}$Co ${\to }^{56}$Fe (with the half-lives of 6.08 days and 77.2 days) plays a subdominant role to the late-time heating because of their energy deposition exclusively in the form of γ-rays with the diminishing f_th for t > 10 days (dotted lines in the left panel of Figure 17).

In our explored models, the dynamical ejecta have negligible contributions to the radioactive heating that powers the kilonova emission. Therefore, the comparison of $\dot{Q}$ with the bolometric luminosity of the kilonova solely does not provide any constraint on the production of heavy r-process elements. However, the mass fraction of lanthanides in each model, X_la, serves as another constraint. The values for the models with α_vis = 0.04, X_la ∼ 0.002–0.004 (last column in Table 3), fall within the observational value, ∼0.001–0.01, in the kilonova ejecta of the merger GW170817 (e.g., Chornock et al. 2017; Waxman et al. 2018; Even et al. 2020). This indicates that the heavy r-process elements synthesized in the dynamical ejecta become the dominant opacity source of the kilonova, although the amount is insufficient to reproduce the solar r-process-like pattern.

The similar behaviors of the $\dot{Q}$ evolution among our explored models (Figure 17) imply that the low-mass binary events (to be measured in the future) likely exhibit GW170817–like kilonova light curves. The brightness depends, however, on the binary mass as well as the true EOS and viscous parameter.