Electromagnetic Counterparts of Binary-neutron-star Mergers Leading to a Strongly Magnetized Long-lived Remnant Neutron Star

Kyohei Kawaguchi; Sho Fujibayashi; Kenta Hotokezaka; Masaru Shibata; Shinya Wanajo

doi:10.3847/1538-4357/ac6ef7

1. Introduction

The merger of neutron stars (NSs) is currently one of the most interesting multimessenger phenomena in high-energy astrophysics, in which physical processes in extreme (strongly self-gravitating, high-density, and high-temperature) environments are realized. During the inspiral motion just before the merger, gravitational waves reflecting the mass spin, and the internal structures of the NSs are emitted. At the onset of the merger, the NS matter is ejected by tidal disruption and collisional shock heating (e.g., Rosswog et al. 1999; Ruffert et al. 2001; Hotokezaka et al. 2013). After the binary merger, a massive neutron star (MNS) or a black hole (BH) surrounded by a strongly magnetized hot and dense accretion torus is formed (Price & Rosswog 2006; Kiuchi et al. 2018). The accretion torus is considered to launch a relativistic jet and outflows by magnetic pressure and tension, the viscous heating due to magnetohydrodynamical turbulence, and the neutrino irradiation. In addition, the r-process nucleosynthesis of heavy elements is expected to proceed in the neutron-rich ejecta (Lattimer & Schramm 1974; Eichler et al. 1989; Freiburghaus et al. 1999; Cowan et al. 2021). In such a situation, the weak interaction processes play an important role in determining the thermodynamic properties of the merger remnants, the post-merger environment, and the abundance of the elements synthesized in the ejecta (e.g., Metzger et al. 2010; Goriely et al. 2010; Wanajo et al. 2014; Just et al. 2015; Sekiguchi et al. 2015, 2016; Radice et al. 2016; Miller et al. 2019; Fujibayashi et al. 2018, 2020a, 2020b, 2020c; Foucart et al. 2020; Just et al. 2021).

The observations of gravitational waves from an NS–NS merger (GW170817; Abbott et al. 2017a) and its multiwavelength electromagnetic (EM) counterparts (Abbott et al. 2017b) have demonstrated that such observations provide an invaluable opportunity to study the fundamental physics under extreme conditions. For example, the gravitational wave detection enables us to infer the degree of the NS tidal deformation in close orbits, leading to constraining the NS equation of state (EOS; Abbott et al. 2017a; De et al. 2018). Observations in the optical and near-infrared (NIR) wavelengths (Andreoni et al. 2017; Arcavi et al. 2017; Coulter et al. 2017; Cowperthwaite et al. 2017; Díaz et al. 2017; Drout et al. 2017; Evans et al. 2017; Hu et al. 2017; Valenti et al. 2017; Kasliwal et al. 2017; Lipunov et al. 2017; Pian et al. 2017; Pozanenko et al. 2018; Smartt et al. 2017; Tanvir et al. 2017; Troja et al. 2017; Utsumi et al. 2017) are consistent with the theoretical predictions for a kilonova; the emission powered by the radioactive decay heating of r-process elements (Li & Paczynski 1998; Kulkarni 2005; Metzger et al. 2010; Kasen et al. 2013; Tanaka & Hotokezaka 2013). The rapid color evolution and spectra strongly suggest that indeed the r-process elements have been synthesized in the ejecta (e.g., Cowperthwaite et al. 2017; Kasen et al. 2017; Kasliwal et al. 2017; Perego et al. 2017; Tanaka et al. 2017; Villar et al. 2017; Rosswog et al. 2018; Kawaguchi et al. 2018).

A synchrotron afterglow was discovered in the radio, optical, and X-ray bands (Hallinan et al. 2017; Haggard et al. 2017; Margutti et al. 2017; Troja et al. 2017; Lyman et al. 2018). A superluminal motion of the radio source reveals that a narrowly collimated relativistic jet is produced in GW170817 (Mooley et al. 2018). In addition to the jet afterglow, a long-term synchrotron flare is expected to arise from the merger ejecta on timescales of ∼1–10⁴ yr (Nakar & Piran 2011; Hotokezaka & Piran 2015; Hotokezaka et al. 2018; Margalit & Piran 2020). Interestingly, the X-ray flux observed around 3.5 yr after GW170817 exceeds the expected flux of the jet afterglow, suggesting that the ejecta afterglow starts to dominate over the jet component (Hajela et al. 2022; Troja et al. 2022). However, the origin of the X-ray excess is sill under debate because the radio flux is still consistent with the prediction of the jet afterglow (Balasubramanian et al. 2021).

There are a variety of open questions associated with GW170817. One of the most important questions is whether and, if so, when the remnant NS has gravitationally collapsed into a BH. This question is connected to the determination of the maximum mass of the NS, which places a constraint on the NS EOS (e.g., Margalit & Metzger 2017; Rezzolla et al. 2018; Shibata et al. 2019). The amount of the ejecta mass and its ejection mechanism are also suggestive for understanding the underlying physics (e.g., Bauswein et al. 2017; Coughlin et al. 2018; Radice & Dai 2019; Kiuchi et al. 2019), while they are also still in debate. To address these questions for GW170817 and also for the observations of the future events, a quantitative understanding of the relation between the mass ejection mechanism and the EM counterparts is important. In this article, among the various types of EM counterparts for the NS mergers, we focus particularly on the EM counterparts that are launched by the consequence of subrelativistic and mildly relativistic mass ejection associated with the presence of a long-lived MNS.

Various studies in terms of the numerical simulations have been performed for NS mergers, revealing the detailed property of the ejecta and the resulting element abundances of the nucleosynthesis together with the dependence on the mass ejection mechanism and the binary parameters, such as the NS mass and NS EOS (Hotokezaka et al. 2013; Bauswein et al. 2013; Wanajo et al. 2014; Sekiguchi et al. 2015; Foucart et al. 2016; Sekiguchi et al. 2016; Radice et al. 2016; Dietrich et al. 2017; Bovard et al. 2017; Kiuchi et al. 2018; Dessart et al. 2009; Metzger & Fernández 2014; Perego et al. 2014; Just et al. 2015; Wu et al. 2016; Siegel & Metzger 2017; Shibata et al. 2017; Lippuner et al. 2017; Fujibayashi et al. 2018; Siegel & Metzger 2018; Ruiz et al. 2018; Fernández et al. 2019; Christie et al. 2019; Perego et al. 2019; Miller et al. 2019; Fujibayashi et al. 2020a, 2020b, 2020c; Bernuzzi et al. 2020; Ciolfi & Kalinani 2020; Nedora et al. 2021; Foucart et al. 2020; Fernández et al. 2020; Mösta et al. 2020; Shibata et al. 2021a, 2021b; see Shibata & Hotokezaka 2019 for a review). Based on or motivated by the knowledge of the ejecta profile and the element abundances obtained by those simulations, radiative transfer simulations with the realistic heating rate and/or the detailed opacity calculations have been performed to predict the kilonova light curves in the last decade (e.g., Kasen et al. 2013, 2015; Barnes et al. 2016; Wollaeger et al. 2018; Tanaka et al. 2018; Wu et al. 2019; Kawaguchi et al. 2018; Hotokezaka & Nakar 2020; Kawaguchi et al. 2020; Korobkin et al. 2021; Bulla et al. 2021; Zhu et al. 2021; Barnes et al. 2021; Nativi et al. 2020; Kawaguchi et al. 2021; Wu et al. 2022; Just et al. 2022; Curtis et al. 2021).

Here, it is to be emphasized that, for accurately deriving the light curve and spectrum of kilonovae, a long-term hydrodynamics evolution with the timescale of ≫10 s is crucial. At the time of the ejecta formation (≲10 s), the ejected matter still has a nonnegligible amount of internal energy compared to its kinetic energy, and thus, the subsequent ejecta trajectory can be modified by the thermal pressure (Kastaun & Galeazzi 2015). This hydrodynamical effect in the subsequent evolution could be reflected in the light curve of the EM counterparts. Currently, there are only a limited number of studies that consider the long-term hydrodynamics evolution of ejecta to predict the property of the EM counterparts (Rosswog et al. 2014; Grossman et al. 2014; Fernández et al. 2015, 2017; Foucart et al. 2021; Kawaguchi et al. 2021; Wu et al. 2022).

Recently, Shibata et al. (2021b) performed a long-term simulation for binary-neutron-star (BNS) merger remnants using a general-relativistic radiation magnetohydrodynamics (GRRMHD) simulation code with a mean-field dynamo term (Shibata et al. 2021a). They paid attention to the models in which the remnant MNS survives for a long period after the merger with ≳3 s. For these systems, it is indicated that the intrinsic magnetohydrodynamics (MHD) effects such as the magnetocentrifugal effect (Blandford & Payne 1982) and magnetic-tower effect may play an important role, for the mass ejection if the magnetic-field strength in the merger remnant is amplified, and the high field strength is preserved for a timescale longer than the mass ejection timescale. For example, the magnetic field lines anchored at the remnant MNS induce the magnetocentrifugal effect (Blandford & Payne 1982) to the surrounding matter, which receives angular momentum from the MNS and a part of which becomes a fast outflow. In this paper, we employ the numerical results of Shibata et al. (2021b) to predict a variety of the hypothetical EM counterparts. We find that, in the presence of the significant MHD effects, a synchrotron afterglow with high luminosity can be associated with the merger event. We also perform a radiative transfer simulation for kilonovae and find that they can also be bright in the early epoch, while it shows the optical emission rapidly declining in a few days and the long-lasting (∼10 days) emission which is very bright in the NIR wavelength.

This paper is organized as follows: In Section 2, we describe the method employed in this study. In Section 3, we describe the BNS models we study in this work. The results of the nucleosynthesis calculation are also presented in this section. In Section 4, we present the property of the ejecta obtained by the long-term hydrodynamics evolution and the numerical results for the light curves of the EM counterparts. Finally, we discuss the implication of this paper in Section 5. Throughout this paper, c denotes the speed of light.

2. Method

To obtain the ejecta profile in the homologously expanding phase, we evolve the ejecta for ≈0.1 d by solving the axisymmetric hydrodynamics equations employing the outflow data obtained by numerical-relativity (NR) simulations as the boundary condition (Shibata et al. 2021b). In the following, to distinguish between the present simulation and NR simulation, we refer to the present hydrodynamics simulations as the HD simulations.

We basically use the same method as that described in Kawaguchi et al. (2021) for the HD simulations except for an update on the treatment of the radioactive heating effect. In this method, relativistic hydrodynamics equations in the spherical coordinates are solved taking into account the effect of fixed-background gravity of a nonrotating BH metric in the isotropic coordinates. We employ the ideal-gas EOS with the adiabatic index of Γがんま = 4/3. In the previous work, the radioactive heating effect is incorporated by adding a source term in the energy and momentum equations. However, as noted in Appendix A of Kawaguchi et al. (2021), this naive treatment violates the energy and momentum conservation of the system: Because the radioactive heating is the consequence of the release of the binding energy of nuclei, the total energy and momentum of the system do not change with this process. Although the error due to the naive treatment is small for our setup, we have corrected our HD simulation code to self-consistently solve the system. The detail of the modification is summarized in Appendix A.

For the HD simulations, the uniform grid spacing with N_θしーた grid points is prepared for the polar angle θしーた, while for the radial direction, the following nonuniform grid structure is employed; the jth radial grid point is given by

$\begin{eqnarray}&&\mathrm{ln}\,{r}_{j}=\mathrm{ln}\left(\displaystyle \frac{{r}_{\mathrm{out}}}{{r}_{\mathrm{in}}}\right)\displaystyle \frac{j-1}{{N}_{r}}+\mathrm{ln}\,{r}_{\mathrm{in}},j=1\cdots {N}_{r}+1,\end{eqnarray} \tag{ 1 }$

where r_in and r_out denote the inner and outer radii of the computational domain, respectively, and N_r denotes the total number of the radial grid points. We set the same time origin for the HD simulations as in the NR simulations for the post-merger evolution.

The time-sequential hydrodynamics property of the outflow is extracted from the NR simulations of Shibata et al. (2021b) at a certain radius, and is imposed at r = r_in of the HD simulations as the inner boundary condition. Accordingly, r_in is initially set to agree with the radius at which the outflow profile of the NR simulations is extracted. After the NR simulation data are run out at t ∼ 5 s, the HD simulation is continued by setting a very small floor-value to the density of the inner boundary. When the high-velocity edge of the outflow reaches the outer boundary of our HD simulation, the radial grid points are added to the outside of the original outer boundary. At the same time, the innermost radial grid points are removed to keep the total number of the radial grid points. By this prescription, the value of r_in is increased in the late phase of the HD simulations, and r_out is always set to be 10³ r_in. We note that the total mass in the removed grids lost by this prescription is always negligible (≲4 × 10⁻⁴ M_⊙) compared to the post-merger ejecta present around the inner boundary (∼0.1 M_⊙).

For calculating the kilonova light curves and the synchrotron emission that arise from the ejecta fast tail, the ejecta profile at 0.1 d after the onset of the merger is employed. We perform the HD simulations with a larger grid resolution than the previous study, that is, (N_r, N_θしーた) = (3192, 196) to resolve the fast tail of the ejecta with relatively low density and with mildly relativistic velocity. To predict the velocity distribution, which is employed to calculate the nonthermal emission, we switch off the radioactive heating effect in this HD simulation to reduce the computational time. This is justified by the fact that the energy deposition due to radioactive heating is negligible compared to the kinetic energy for the fast tail. Indeed, we confirmed that the light curves of the nonthermal emission shown below are not affected by the presence of the heating terms in the HD simulations (see Appendix B). For calculating the kilonova light curves, we employ the same NR ejecta profile data as that for calculating the fast-tail, but the HD simulations are performed taking the radioactive heating into account, because the internal energy deposited by the radioactive heating during the ejecta expansion partially contributes to the kilonova light curves (Kawaguchi et al. 2021). We note that a relatively small grid resolution of (N_r, N_θしーた) = (1024, 128) is employed for the viscous model (see the next subsection), because we find that it is sufficient to resolve the ejecta with a small amount of the fast tail components (v > 0.6 c with v being the radial velocity).

3. Model

3.1. Numerical Relativity Simulation

As the input for the HD simulations, we employ the outflow profiles obtained by NR simulations for a post-merger evolution of a BNS in Fujibayashi et al. (2020c); Shibata et al. (2021b). The key quantities of each model are summarized in Table 1. The first five models listed in Table 1 are obtained by long-term resistive MHD simulations with a mean-field dynamo term (Shibata et al. 2021b), employing the axisymmetrized merger remnant profile of an equal-mass BNS merger simulation with the DD2 EOS (Banik et al. 2014) and with each NS mass of 1.35 M_⊙; see Table 1 for the employed conductivity σしぐま_c and mean-field dynamo parameters αあるふぁ_d. We note that, while the magnetic field energy reaches ∼10⁵¹ erg (∼1% of the kinetic energy of the system) at the peak for all the MHD models, for the larger values of σしぐま_c the dissipation timescale of the magnetic fields is longer, and for the larger values of αあるふぁ_d, the amplification timescale of the magnetic-field strength associated with the dynamo action is shorter. For the αあるふぁ_d and σしぐま_c values employed in the models listed in Table 1, the magnetic-field amplification timescale and dissipation timescale are found to be ∼10 ms and ∼1 s (σしぐま_c/10⁸ s⁻¹), respectively (see Section IIB and Figure 9 in Shibata et al. 2021b). The last model listed in Table 1 with the label of "αあるふぁ = 0.04" denotes the result for the same BNS model but of a viscous-hydrodynamics simulation with the dimensionless alpha viscous parameter of αあるふぁ = 0.04 (Fujibayashi et al. 2020c).

Table 1. Key Model Parameters

Model	(σしぐま_c [s⁻¹], αあるふぁ_d)	${M}_{\mathrm{eje}}^{0}({M}_{\mathrm{eje}}^{0.1\,})\,[{10}^{-2}{M}_{\odot }]$	${E}_{\mathrm{eje}}({E}_{{\rm{K}},\mathrm{eje}}^{0.1\,{\rm{d}}})\,[{10}^{51}{M}_{\odot }{c}^{2}]$	X_lan+act	X_Y+Zr
MNS70a	(1 × 10⁷, 1 × 10⁻⁴)	4.6 (3.7)	1.0 (0.84)	6.9 × 10⁻³	7.6 × 10⁻²
MNS70b	(1 × 10⁷, 2 × 10⁻⁴)	7.8 (6.1)	2.1 (1.7)	3.3 × 10⁻³	9.3 × 10⁻²
MNS75a	(3 × 10⁷, 1 × 10⁻⁴)	9.4 (8.4)	14 (13)	4.2 × 10⁻³	1.2 × 10⁻¹
MNS75b	(3 × 10⁷, 2 × 10⁻⁴)	12 (12)	16 (15)	1.2 × 10⁻²	1.2 × 10⁻¹
MNS80	(1 × 10⁸, 1 × 10⁻⁴)	13 (13)	41 (42)	9.0 × 10⁻³	1.3 × 10⁻¹
viscous (αあるふぁ = 0.04)	⋯	7.6(6.5)	0.57(0.52)	3.1 × 10⁻³	5.0 × 10⁻²

Note. The columns describe the model name, the phenomenological dynamo parameters of the NR simulations (see Shibata et al. 2021b), ejecta mass, ejecta energy, lanthanide+actinide mass fraction, and Y+Zr mass fraction, respectively. The values for the ejecta mass and energy shown outside and inside the parenthesis denote the values calculated from the input outflow data of the NR simulations and the results of the HD simulations at t = 0.1 day, respectively. The values listed here are obtained by the results for the HD simulations without taking the radioactive heating into account.

Download table as: ASCII Typeset image

Table 1 summarizes the ejecta mass, ${M}_{\mathrm{eje}}^{0}$ , and ejecta energy, ${E}_{\mathrm{eje}}^{0}$ , calculated from the time sequence of the input outflow data of the NR simulations by the following equations:

$\begin{eqnarray}&&{M}_{\mathrm{eje}}^{0}:={\int }_{r={r}_{\mathrm{in}}}{\rho }_{* }{v}^{r}{dtdS},\end{eqnarray} \tag{ 2 }$

$\begin{eqnarray}&&{E}_{\mathrm{eje}}^{0}:={\int }_{r={r}_{\mathrm{in}}}{\rho }_{* }{c}^{2}({hw}-1){v}^{r}{dtdS}.\end{eqnarray} \tag{ 3 }$

Here, ${\rho }_{* }:=\rho {u}^{t}\sqrt{-g}$ with ρろー the rest-mass density, u^t the time component of the four velocity u^μみゅー, and g the determinant of the spacetime metric. v^r=u^r/u^t, w=αあるふぁ_g u^t with αあるふぁ_g the lapse function, and h is the specific enthalpy.⁶ We note that ${E}_{\mathrm{eje}}^{0}$ contains the contribution from both the kinetic- and internal-energy flows. We note that all the matter crossing the inner radius is included in the integration. The contributions from the dynamical ejecta to the total ejecta mass and ejecta energy, which are determined only by the merger phase, are ≈1.5 × 10⁻³ M_⊙ and ≈5 ×10⁴⁹ erg, respectively (Fujibayashi et al. 2020c). Thus, the post-merger ejecta contributes primarily to these quantities.

As found from Table 1, the ejecta become massive and more energetic as the value of αあるふぁ_d or σしぐま_c increases (i.e., as the post-merger MHD activity is enhanced). In particular, the increase in the ejecta energy is significant as the value of αあるふぁ_d increases, while the difference in the ejecta mass is within a factor of 3 among the models. Since the rotational kinetic energy of the MNS (∼10^52–53 erg) is the main source to accelerate the ejecta through the magnetocentrifugal effect (Blandford & Payne 1982) in the MHD simulations, these results indicate that the efficiency of converting the rotational kinetic energy of the MNS to the kinetic energy of the ejecta depends strongly on the magnetic-field evolution. Specifically, the dissipation timescale of the magnetic fields controlled by the conductivity determines the efficiency. On the other hand, the ejecta mass converges within 20%–50% of the disk mass irrespective of the dynamo parameters (Shibata et al. 2021b). The ejecta mass for the viscous model (with αあるふぁ = 0.04) is as large as that for model MNS70b and larger than that for model MNS70a. However the ejecta energy for the viscous model has only a half of that for MNS70a, for which the MHD activity is weakest among the present models. This also shows that the kinetic energy of the ejecta is enhanced by an intrinsic MHD effect.

Figure 1 shows the profiles of the rest-mass density and electron fraction at t = 0 of the HD simulations (note that the HD and NR simulations for the post-merger evolution share the same time origin). We note that, for the electron fraction profile, the value at the temperature of 5 × 10⁹ K is shown. In this initial time slice, only the dynamical ejecta component is present in the computational domain. The rest-mass density profile of the dynamical ejecta exhibits weakly spheroidal morphology, although the electron fraction is distributed in an anisotropic manner. The ejecta in the polar region of θしーた ≲ 45°–60° are dominated by the matter with Y_e ≥ 0.3, while those around the equatorial region (θしーた ≳ 60°) are dominated with Y_e ≤ 0.2. The difference in the ejecta electron fraction reflects the difference in the mass ejection mechanisms: the polar component is driven by the collisional shock heating during the merger of two NSs, while the other component is driven by the tidal torque induced by the nonaxisymmetric matter distribution at the onset of the merger (Wanajo et al. 2014; Sekiguchi et al. 2015, 2016; Radice et al. 2016).

3.2. Nucleosynthesis

We performed a nucleosynthesis calculation in each ejecta particle for a given model. The method of the nucleosynthesis calculation as well as the post-process particle tracing is the same as those described in our latest paper (Fujibayashi et al. 2020c). In this work, we set the radius for the extraction of the ejecta information at 1500 km, which is the initial value of r_in for the HD simulations. The tracer particles are distributed at this radius (r = 1500 km) uniformly in the polar angle, which are traced back in time. The bulk of the dynamical ejecta has already passed the extraction radius at the beginning of the NR (MHD) simulations (see Figure 1). We adopt, therefore, the result of the nucleosynthesis calculation in the NR (viscous) model of Fujibayashi et al. (2020c) for the matter with initial radius 1500 km < r < 8000 km.

Figure 2 shows the results of the nucleosynthesis calculation for all the models listed in Table 1. The left and right panels present the mass fractions as functions of atomic mass number A (after decay) and of atomic number Z (at 1 day), respectively, the latter being used for the radiative transfer simulations in this study. The heavy r-process elements are underproduced than the solar-abundance pattern (r-process residuals to the solar abundance taken from Prantzos et al. 2020) not only for the viscous-hydrodynamics model but also for the MHD models irrespective of the dynamo parameters. This is because the mass of the post-merger ejecta is more than one order of magnitude larger than that of the dynamical ejecta, and the post-merger ejecta predominantly synthesizes relatively light r-process elements with the atomic mass number smaller than 130. On the other hand, significant amounts of Y (Z = 39) and Zr (Z = 40), which contribute a lot to the opacity in the optical wavelength (Tanaka et al. 2020; Kawaguchi et al. 2021; Ristic et al. 2022), are synthesized (see also Table 1 for the lanthanide+actinide and Y+Zr mass fraction of the ejecta). Note that the MHD models synthesize a slightly larger amount of r-process elements (including lanthanides+actinides and Y+Zr) than that in the viscous model (αあるふぁ = 0.04) owing to the higher entropies in the ejecta (and in part higher outgoing velocities) for the formers (Shibata et al. 2021b).

As we discussed in Fujibayashi et al. (2020c), the abundance pattern shown in Figure 2 is universally found in the presence of long-lived MNSs as a merger remnant. This abundance pattern does not agree with the solar-abundance pattern because of the overproduction of the relatively light elements with A ≲ 130 (underproduction of the heavy elements). Therefore, suppose that the solar-abundance pattern is universal in the universe and the main site for the r-process nucleosynthesis is BNS mergers; the channel in which a long-lived MNS is formed should be the minority of the BNS mergers (but see, Wu & Tamborra 2017; Wu et al. 2017; Li & Siegel 2021, for the case that the neutrino oscillation plays an important role in determining the ejecta Y_e). In other words, if the future observation suggests that the formation of the long-lived remnant MNSs is the majority of the BNS mergers, it will indicate that BNS mergers might not be the main production sites of the r-process elements. Besides the event rates, however, this merger channel synthesizes the rich r-process elements and, thus, is the promising source of bright kilonovae, which may be observed even if they appear at a large distance. We show the light-curve models of kilonovae in Section 4.3.

4. Results

4.1. Ejecta Profiles

This subsection focuses on the results of the ejecta profile in the HD simulation. As we found in the previous study, a fraction of matter that passed through the extraction radius in the NR simulations fails to become gravitationally unbound during the subsequent hydrodynamics evolution even if the often-used Bernoulli criteria ( $\equiv {{hu}}_{t}\leqslant -{h}_{\min }$ , see Kastaun & Galeazzi 2015; Vincent et al. 2020; Fujibayashi et al. 2020c; Foucart et al. 2021) was satisfied at the time of passing through. The reason for this is that the pressure of the precedingly outgoing matter can decelerate the matter that catches up. Hence, the total rest mass and energy of the ejecta at the homologously expanding phase can be slightly different from those measured in the NR simulations.⁷

To quantify the rest mass and kinetic energy of the ejecta at the homologously expanding phase, we evaluate the following quantities at t = 0.1 d:

$\begin{eqnarray}&&{M}_{\mathrm{eje}}^{0.1\,{\rm{d}}}:={\left.\int {\rho }_{* }{dV}\right|}_{t=0.1\,{\rm{d}}}\end{eqnarray} \tag{ 4 }$

$\begin{eqnarray}&&{E}_{{\rm{K}},\mathrm{eje}}^{0.1\,{\rm{d}}}:={\left.\int {\rho }_{* }{c}^{2}(w-1){dV}\right|}_{t=0.1\,{\rm{d}}}.\end{eqnarray} \tag{ 5 }$

Here, we simply measure the mass and the energy in the computational domain of the HD simulation, because, at t = 0.1 d, approximately all the fluid elements satisfy the geodesic criteria of unbound matter ( ≡ u_t < −1). The numerical results are listed in Table 1, which are obtained without taking the radioactive heating into account in the HD simulations. We note that the ejecta mass increases by 10%–15% if we take the radioactive heating into account, although it has only a negligible contribution to the EM light curves.

The differences between ${M}_{\mathrm{eje}}^{0}$ and ${M}_{\mathrm{eje}}^{0.1\,{\rm{d}}}$ and between ${E}_{\mathrm{eje}}^{0}$ and ${E}_{{\rm{K}},\mathrm{eje}}^{0.1\,{\rm{d}}}$ denote the rest mass and energy of the fall-back matter. Note that all the matter passing through the inner boundary (including that finally falls back) is counted to ${M}_{\mathrm{eje}}^{0}$ and ${E}_{\mathrm{eje}}^{0}$ and satisfies the Bernoulli criteria ( $\equiv {{hu}}_{t}\leqslant -{h}_{\min }$ , see Fujibayashi et al. 2020c ) at the time of passing through. These quantities decrease as the values of αあるふぁ_d and σしぐま_c increase, i.e., the kinetic energy (and velocity) of the ejecta becomes higher. The fact that ${E}_{{\rm{K}},\mathrm{eje}}^{0.1\,{\rm{d}}}/{E}_{\mathrm{eje}}^{0}$ is larger than ${M}_{\mathrm{eje}}^{0.1\,{\rm{d}}}/{M}_{\mathrm{eje}}^{0}$ implies that the specific energy of the fall-back matter is smaller than the bulk specific energy of the ejecta.

We note that, in the HD simulations, after the outflow data of the NR simulations run out, the rest-mass density on the inner boundary is switched to that for an atmosphere value (= a sufficiently small floor-value). We should note that this treatment could artificially increase the mass of the fall-back matter due to the sudden vanishing of the pressure support on the inner boundary. Nevertheless, as found in our previous paper Kawaguchi et al. (2021), the contribution of such marginally unbound matter to either the synchrotron afterglow or the kilonova emission is minor because it has only low velocity.

4.1.1. 2D Ejecta Profile

Figures 3 and 4 show the profiles of the rest-mass density and electron fraction of the ejecta at t = 0.1 d for the viscous model (αあるふぁ = 0.04) and for models MNS70a, MNS75a, and MNS80, respectively. We note that, as is the case for Figure 1, the value at the temperature of 5 × 10⁹ K is shown for the electron fraction profile. We note that the ejecta profiles shown in Figures 3 and 4 are obtained in the HD simulations that take the radioactive heating into account.

**Figure 3.** Profiles of ejecta rest-mass density (the left panel) and electron fraction (the right panel) at t = 0.1 d for the viscous model (*αあるふぁ* = 0.04). The value at the temperature of 5 × 10⁹ K is shown for the electron fraction profile.
Download figure:
Standard image High-resolution image

**Figure 4.** The same as Figure 4 but for the MHD models. The top, middle, and bottom panels denote the results for models MNS70a, MNS75a, and MNS80, respectively.
Download figure:
Standard image High-resolution image

As Shibata et al. (2021b) pointed out, the outflow property for the MHD models with a small value of the conductivity, σしぐま_c = 1 × 10⁷ s⁻¹, is qualitatively similar to that for the viscous model. Thus, we first discuss the results for the viscous model as the reference; and then, we compare the results for the MHD models with those for the viscous model to clarify the specific features of the MHD effects.

The rest-mass density and electron fraction profiles for the viscous model (αあるふぁ = 0.04) show approximately the same features as those found for the viscous model in the previous study (a low-mass model, DD2-125; Fujibayashi et al. 2020c; Kawaguchi et al. 2021). The density profile of the ejecta with r/ct ≲ 0.1 exhibits a weakly spheroidal shape. This component is driven by the viscous effect from the torus surrounding the central MNS. On the other hand, the density profile of the high-velocity ejecta with r/ct ≳ 0.1 exhibits mildly prolate morphology: for example, the rest-mass density larger than 10⁻¹³ g cm⁻³ at t = 0.1 d is distributed for the velocity space up to r/ct ≈ 0.5 and ≈0.4 toward the polar and equatorial directions, respectively. For the Y_e distribution, primarily two distinct components are seen: one has a torus-like shape located around x/ct ≈ 0.2–0.4 and θしーた ≳ 45° with Y_e ≲ 0.25, and the other has a prolate shape located around x/ct ≲ 0.2 and θしーた ≲ 45° with Y_e ≳ 0.3. The former corresponds to the dynamical ejecta component driven by tidal torque induced by the nonaxisymmetric matter distribution at the onset of the merger. The latter is composed of a part of the dynamical ejecta driven by the shock heating at the merger and neutrino-irradiation-enhanced ejecta (Fujibayashi et al. 2020c; Kawaguchi et al. 2021).

In contrast to the viscous model, an appreciable amount of the ejecta with high velocity (≳0.5 c) is present in the MHD models, in particular for σしぐま_c ≥ 3 × 10⁷ s⁻¹, for which the dissipation time of the magnetic field in the MNS is longer than or as long as the timescale of the post-merger mass ejection of several hundreds milliseconds. Figure 4 shows that the high-velocity post-merger ejecta push the dynamical ejecta. The pushed ejecta component is highly appreciable for σしぐま_c ≥ 3 × 10⁷ s⁻¹, and the highest velocity of the ejecta with the density larger than 10⁻¹³ g cm⁻³ at t = 0.1 d reaches ≳0.9 c in the polar direction. The velocity of the post-merger ejecta can reach ∼0.7 c even in the equatorial direction. As the consequence, the low-Y_e dynamical ejecta is accelerated and confined toward the equatorial plane, and also, the ejecta rest-mass density profile settles into a more isotropic shape than that for the viscous model. By contrast, the MHD model with a relatively small value of σしぐま_c = 1 × 10⁷ s⁻¹ (MNS70a) shows an ejecta profile similar to those for the viscous model, besides the presence of the high-velocity (≳0.5 c) matter in the polar region. This similarity is due to the fact that for such a small value of σしぐま_c the dissipation timescale of the magnetic fields in the MNS is shorter than the post-merger mass-ejection timescale, and the acceleration of the post-merger ejecta by the magnetocentrifugal effect plays only a minor role.

4.1.2. Kinetic Energy Distribution of Ejecta

Figure 5 shows the kinetic energy distributions of ejecta at t = 0.1 d as a function of the radial four-velocity defined by

$\begin{eqnarray}&&{E}_{{\rm{K}}}(\gt {u}^{r})={\int }_{\gt {u}^{r}}{\rho }_{* }{c}^{2}(w-1){dV}.\end{eqnarray} \tag{ 6 }$

As we already found in Table 1, the ejecta become more energetic as the values of αあるふぁ_d and/or σしぐま_c increase. The kinetic energy distributions for the MHD models are composed broadly of two components. One is a trans-relativistic component with the Lorentz factor of ≲2, and the other is a mildly relativistic component with the Lorentz factor of ≈2–5. The former component dominates the total kinetic energy of the ejecta with ∼10⁵¹ erg for σしぐま_c = 1 ×10⁷ s⁻¹ (MNS70a and MNS70b) and ≳10⁵² erg for σしぐま_c ≥ 3 ×10⁷ s⁻¹ (MNS75a, MNS75b, and MNS80), respectively. The latter component has ∼10⁴⁸–10⁴⁹ erg, which is by more than ∼2 orders of magnitudes smaller than the former components.

The bottom panel of Figure 5 presents the angular dependence of the kinetic energy. This shows that the kinetic energy of the ejecta is dominated by the matter in a polar region (θしーた ≤ 60°). In particular, the mildly relativistic component with u^r/c ≳ 3(v^r/c ≳ 0.95) are only present in the polar region with θしーた ≤ 30°. This profile is obtained by the following mechanism: the magnetocentrifugal effect initially induces the acceleration of the post-merger ejecta primarily toward the equatorial direction. However, the straightforward acceleration is blocked by the torus surrounding the MNS. By the interaction with the torus, shocks are formed around the torus surface, and thus, the motion of the post-merger ejecta is changed to the polar direction.

In the viscous model, the kinetic energy distribution of ejecta for u^r/c > 0.05 is similar to that only with the dynamical ejecta (Hotokezaka et al. 2018). This is due to the fact that the post-merger ejecta have only a minor contribution to the fast tail of the ejecta for the viscous model. This is consistent with the fact that the average velocity of the post-merger ejecta is lower than that for the MHD models.

4.2. Synchrotron Afterglow

The enhancement of the rest mass and total kinetic energy of the fast velocity component in the MHD models has a significant impact on the flux of the ejecta afterglow, because the light curve depends strongly on the Lorentz factor and kinetic energy of the ejecta (see, e.g., Nakar & Piran 2018). It has been suggested that the afterglow of a BNS merger is expected to be extremely bright if a fraction of the remnant NS's rotational energy ∼10⁵² erg is converted to the ejecta kinetic energy (Yu et al. 2013; Metzger & Bower 2014; Horesh et al. 2016; Yu et al. 2018; Li et al. 2018; Beniamini & Lu 2021). Here we present the synchrotron afterglow light curve arising from the merger ejecta using the kinetic energy distribution shown in Figure 5.

The light curves for the synchrotron afterglow are calculated using the same analysis as that in Hotokezaka et al. (2018) with the equipartition parameters at the shocks, _e (the fraction of the energy that accelerated electrons have) and _B (the fraction of the energy that the magnetic field has), set to be 0.1 and 0.01, respectively. The energy distribution of the accelerated electrons is assumed to be a power-law distribution of the electron Lorentz factor with the power of p = 2.2. Note that Margalit & Quataert (2021) show that thermal electrons are expected to contribute significantly to the radio flux for mildly relativistic outflows.

Figure 6 shows the light curves of the kilonova radio afterglow at 3 GHz. The top, middle, and bottom panels denote the results for the cases of the different interstellar medium (ISM) densities, n = 10⁻¹, 10⁻³, and 10⁻⁵ cm⁻³, respectively. Broadly speaking, there are two components in the radio light curves: one is the main component that peaks at t ≥ 0.1 yr, and the other is the faint early rising component present for t ≤ 0.1 yr. The emission around the peak luminosity arises from the trans-relativistic ejecta component with the Lorentz factor of ≲2. The mildly relativistic component with the Lorentz factor of ≈3–5 contributes to the early radio light curves for t ≤ 0.1 yr. The contribution of the latter component to the total radio luminosity is minor compared to the former component, and hence, we focus mainly on the former one in the following, although the latter could be important for the follow-up observation in the early phase.

epsilon — **Figure 6.** The light curves of the kilonova radio afterglow at 3 GHz for the MHD models at the distance to the source of 200 Mpc. The top, middle, and bottom panels denote the results for the cases that the ISM density, n, is 10⁻¹, 10⁻³, and 10⁻⁵ cm⁻³, respectively. Here the parameters of _e = 0.1, _B = 0.01, and p = 2.2 are used for modeling the shock. We note that the luminosity for the viscous model ("vis × 10²") is scaled-up by a factor of 100. The data points denote the observation of GW170817 taken from Makhathini et al. (2021).
Download figure:
Standard image High-resolution image

For either of the cases, the radio emission around the peak time in the MHD models is much brighter than that predicted by the model only with the dynamical ejecta (Hotokezaka et al. 2018) and by the viscous model. The peak luminosity varies by more than two orders of magnitude for a fixed ISM density, and the larger luminosities are achieved for the larger values of σしぐま_c and/or αあるふぁ_d. Notably, the peak flux of the models with σしぐま_c = 1 × 10⁸ s⁻¹ exceeds ∼0.1 mJy even for the distance to the source of 200 Mpc and a very low density n ≲ 10⁻⁵ cm⁻³. This peak flux is as bright as the radio afterglow peak observed in GW170817. On the other hand, n ≥ 10⁻³ cm⁻³ is required for the radio emission to be brighter than ∼0.1 mJy for the models with σしぐま_c = 1 × 10⁷ s⁻¹. This suggests that the radio afterglow will carry important information for the magnetic-field enhancement in the merger remnants with a long-lived MNS.

The time of the peak luminosity also depends on the ISM density and the dynamo parameters. In particular, the peak time is delayed for the models with large values of σしぐま_c. This reflects the fact that larger total kinetic energy of the ejecta is achieved for the larger values of σしぐま_c because of the longer dissipation timescale of the magnetic fields (see Table 1). Interestingly, the rising parts of the radio light curves for the models with σしぐま_c ≥ 3 × 10⁷ s⁻¹ have a shape similar to each other because of the similarity in the kinetic energy distribution for u^r/c ≳ 2. Measuring the slope of the early radio light curves may have an important implication in the magnetic-field amplification.

Figure 7 shows the same as Figure 6 but for that in an X-ray band (1 keV). Here, the case of n = 10⁻³ cm⁻³ is shown. In Figure 7, we also plot the data points for GW170817 obtained by extrapolation from the data at 3 GHz assuming a single power-law spectrum with the best-fit index of −0.584 (Makhathini et al. 2021). As is the same in the radio band, the X-ray emission for the models with σしぐま_c ≥ 3 × 10⁷ s⁻¹ is significantly brighter than that observed in GW170817 for n = 10⁻³ cm⁻³. The X-ray light curves become fainter than the flux extrapolated from the radio light curve with νにゅー^−(p−1)/2 once the cooling frequency crosses 1 keV. This occurs around the peak time for σしぐま_c ≥ 3 × 10⁷ s⁻¹ and around t = 10 yrs for σしぐま_c = 1 × 10⁷ s⁻¹. This timescale is shorter for higher ISM densities, e.g., a week for σしぐま_c ≥ 3 ×10⁷ s⁻¹ at n = 10⁻¹ cm⁻³. By contrast, the X-ray luminosity from the ejecta is not as high as that of GW170817 for σしぐま_c = 1 × 10⁷ s⁻¹ (i.e., for the short dissipation timescale) and the viscous model.

4.3. Kilonovae

The light curves of kilonovae are calculated for models MNS80, MNS75a, MNS70a, and the viscous model using the same method as in our previous paper (Kawaguchi et al. 2021), i.e., using a wavelength-dependent radiative transfer simulation code (Tanaka & Hotokezaka 2013; Tanaka et al. 2017, 2018; Kawaguchi et al. 2020, 2021). In this code, the photon transfer is simulated by a Monte Carlo method for the given ejecta profiles composed of the density, velocity, and element abundance. The time-dependent thermalization efficiency is taken into account following an analytic formula derived by Barnes et al. (2016). The ionization and excitation states are determined under the assumption of the local thermodynamic equilibrium (LTE) by using the Saha ionization and Boltzmann excitation equations. We note, however, that the assumption of the LTE is not always valid in the low-density region of the ejecta for the later epoch. We discuss the possible non-LTE effects to the kilonova light curves in Section 4.4.

For the photon–matter interaction, the bound–bound, bound–free, and free–free transitions and electron scattering are taken into account for the transfer of optical and infrared photons (Tanaka & Hotokezaka 2013; Tanaka et al. 2017, 2018). The formalism of the expansion opacity (Friend & Castor 1983; Eastman & Pinto 1993; Kasen et al. 2006) and the updated line list derived in Tanaka et al. (2020) are employed for the bound–bound transitions. The new line list is constructed by an atomic structure calculation for the elements from Z = 26 to Z = 92 (see Tanaka et al. 2020 for details), and supplemented by the Kurucz's line list for Z < 26 (Kurucz & Bell 1995).

The radiative transfer simulations are performed from t = 0.1 d to 30 d. As the background hydrodynamical state, the density profiles of the HD simulations at t = 0.1 d are used with the assumption of the homologous expansion. The initial internal energy and temperature for the radiative transfer simulations are also determined from those obtained by the HD simulations. The spatial distributions of the heating rate and element abundances are determined by the table obtained by the nucleosynthesis calculations referring to the injected time and angle of the fluid elements. Note that the element abundances at t = 1 d (the right panel in Figure 2) are used during the entire time evolution in the radiative transfer simulations to reduce the computational cost, because this simplified prescription (i.e., neglecting the effects that come from the late-time nucleosynthesis) gives only a minor systematic error on the numerical results.

We first focus on models MNS75a, MNS70a, and the viscous model. Figure 8 shows the bolometric luminosity as a function of time for models MNS75a, MNS70a, and the viscous model (αあるふぁ = 0.04). The bolometric luminosity for MNS75a is always larger than those for MNS70a and the viscous model. This is primarily due to the larger ejecta mass of MNS75a than those of MNS70a and the viscous model, but the higher total specific heating rate and thermalization efficiency are also relevant for the difference. The release of the internal energy stored in the matter until the initial period also contributes to the total bolometric luminosity for t ≤ 0.5 d.

For a given model, the bolometric luminosity becomes larger for smaller viewing angles (except for that of model MNS70a observed from 0° ≤ θしーた ≤ 20° as we discuss later). This is primarily due to the presence of the lanthanide-rich dynamical ejecta in the equatorial region, of which a high opacity plays an important role to block the optical photons coming from the inner dense region (known as the lanthanide curtain effect; Kasen et al. 2015; Bulla 2019; Kawaguchi et al. 2020; Zhu et al. 2020; Darbha & Kasen 2020; Korobkin et al. 2021). The viewing angle dependence is weaker for model MNS75a than those for model MNS70a and the viscous model, because the lanthanide-rich dynamical ejecta are confined more prominently in the equatorial region for model MNS75a than for the other models (cf. Figure 4).

Model MNS70a and the viscous model show similar total bolometric light curves. This is also the case for the bolometric light curves observed from several viewing angles (see the right panel of Figure 8), except for the face-on observation: the bolometric luminosity observed from 0° ≤ θしーた ≤ 20° is as small as that observed from 86° ≤ θしーた ≤ 90° for model MNS70a for t ≲ 1 d. This is due to the dominance of the first r-process peak elements (with a small number of relevant radioactive isotopes) in the polar region. As found in Figure 4, the matter is blown up with high velocity in the MHD models. In addition, appreciable amounts of Y and Zr are present in the matter blown up with high velocity (see Table 1 for the Y+Zr mass fraction in the ejecta). Note that Y and Zr (particularly, those neutral and the first ionized atoms) have significant contributions to the opacity in the optical wavelength (Tanaka et al. 2020; Kawaguchi et al. 2021; Ristic et al. 2022). Photons that diffuse out toward the polar direction have to pass through such a region because Y and Zr are dominantly produced in the high-density part of the ejecta (v/c ≲ 0.3). As a consequence, the emission toward the polar direction is suppressed for model MNS70a in comparison to the viscous model. We note that a small fraction of Sr and lanthanide+actinide elements that are synthesized and blown up in the polar region also contributes to suppressing the optical emission. By contrast, such suppression of the bolometric luminosity is not very significant for model MNS75a. Our interpretation for this is that the release of the internal energy from the ejecta in the high-velocity region, for which the opacity is relatively low, dominates the emission in the early epoch (t ≤ 0.5 d).

Note that model MNS70a and the viscous model show broadly similar light curves regardless of the difference of the total ejecta mass. This is because the difference of the total ejecta mass between the two models is due to the difference in the ejecta density in the region where the velocity is less than 0.05 c, and such a difference has only a minor effect on the light curves up to two weeks, particularly for the peak flux in the optical wavelengths.

Figure 9 shows the gzK-band light curves for models MNS75a, MNS70a, and the viscous model at a hypothetical distance to the source of 40 Mpc. In the early epoch (for t ≤ 0.5 d), the kilonova emission for MNS75a is brighter than or as bright as that observed in GW170817. However, the optical (gri-band) emission fades rapidly, and becomes fainter than that observed in GW170817 for t ≥ 0.5–1 d. On the other hand, the NIR (JHK-band) emission is always brighter than that observed in GW170817. Similar features have been found in our previous study for a model in which the post-merger ejecta is accelerated to a high velocity by a hypothetical remnant activity (see the SMNS case in Kawaguchi et al. 2020).

The rapid fading of the optical emission for model MNS75a stems primarily from the rapid exhaustion of the internal energy deposited at the time of ejecta formation. The ejecta optical thickness decreases due to the rapid outgoing expansion of the ejecta for t ≲ 0.3 d, and the release of the internal energy stored by the radioactive heating (and also that deposited at the time of ejecta formation) dominates the emission. For t ≲ 1 d, such internal energy has already been mostly exhausted by the adiabatic expansion, and as a consequence, the optical emission fades rapidly. The large value of opacity in the optical wavelength due to the presence of the first r-process peak elements (including Y and Zr) in the polar region also plays an important role. In particular, the neutral and the first ionized atoms, of which the fraction increases for t ≥ 1 d, have great contributions to the opacity in the optical wavelengths.

As found in Figure 8, the viewing angle dependence of the gzK-band emission for model MNS75a is weaker than those for the other models, reflecting the approximately spherical ejecta profile and the confinement of the lanthanide-rich dynamical ejecta component in the equatorial region. This feature is advantageous for the observation of the kilonovae from the off-axis directions.

Model MNS70a and the viscous model show approximately identical light curves, except for those observed from the polar direction (0° ≤ θしーた ≤ 20°), reflecting a similar ejecta profile. As we described already, the difference in the polar light curves is primarily due to the different density structure and the different fractions of the first r-process peak elements in the polar region. The similarity in the light curves for the low-σしぐま_c MHD models and the viscous model implies that the results of viscous-hydrodynamics simulations can provide a good phenomenological model for the first-principle MHD model, in the case that the intrinsic MHD effects such as the magnetocentrifugal effect (Blandford & Payne 1982) are not very strong.

Kilonova light curves for model MNS70a and the viscous model are in a fair agreement with those for GW170817. This suggests that, after the merger of the BNS in GW170817, a long-lived MNS with an insufficient magnetic-field amplification or with a short dissipation timescale of the amplified magnetic field (shorter than the mass ejection timescale) might be formed. The weak magnetic-field effect is consistent with the nondetection of bright radio waves associated with the fast ejecta component. However, as we already showed in Figure 2, the abundance patterns of the r-process elements derived for model MNS70a and the viscous model do not agree with the solar-abundance pattern. This suggests that GW170817 might be a rare event of the BNS mergers, under the hypothesis that the BNS mergers are the major site for the r-process nucleosynthesis, and the solar-abundance pattern is universal in the universe (see also the next subsection for the current uncertainty).

Figure 10 compares the kilonova and kilonova afterglow light curves in the r band. The r-band emission from the kilonova afterglow reaches its peak magnitude at 0.1–1 yr. Because the kilonova emission has typically already faded at such an epoch, the afterglow emission can be a characteristic feature for the MHD activity. The afterglow emission could be as bright as the kilonova emission and that of GW170817 intrinsically in the optical wavelength if σしぐま_c ≥ 3 × 10⁷ s⁻¹ and n ≳ 10⁻¹ cm⁻³. For such a case, the emission will be observable even up to the distance of 200 Mpc by 1 m class telescopes (Nissanke et al. 2013). On the other hand, the afterglow emission will be ≈21 and 25 mag in the r band for 40 and 200 Mpc, respectively, if σしぐま_c ≤ 3 × 10⁷ s⁻¹ or n ≲ 10⁻³ cm⁻³, and 4/8 m class telescopes are necessary for the follow-up observation.

Next, we pay attention to model MNS80, in which the MHD effects are most significantly found among the models studied in this work. Figure 11 shows the bolometric luminosity and gzK-band light curves for this model in the same formats as in Figures 8 and 9, respectively, except for the enlarged plot range of the vertical axes. The kilonova emission for model MNS80 is always brighter than that for MNS75a. In particular, for t ≲ 1 d, the total luminosity and isotropic luminosity observed from the polar direction (≤20°) for model MNS80 are larger by factors of 3 and 5, respectively, than those for model MNS75a. We note that these differences are much larger than the difference in the total heating rate (see the top left panel of Figure 11), of which the difference is primarily due to the difference in the total ejecta mass (see Table 1). On the other hand, for t ≳ 1 d, the total bolometric luminosity and isotropic luminosity for the viewing angle larger than 40° for model MNS80 are larger than those for model MNS75a only by a factor of 1.5–2, which approximately corresponds to the difference in the total ejecta mass. This implies that the difference in the luminosity in the late phase is primarily due to the difference in the total ejecta mass.

The bright emission in the early phase of model MNS80 is partially due to the high bulk velocity of the ejecta. Figure 12 shows the ejecta mass distribution at t = 0.1 d as a function of the radial velocity (βべーた = v^r/c). While the ejecta matter distribution is mostly concentrated in the velocity smaller than ≈0.2 c for models MNS75a, MNS70, and the viscous model, the ejecta matter is distributed in a higher-velocity region for model MNS80; the ejecta matter distribution has the peaks at the velocity of ≈0.2 c and ≈0.9 c. The high bulk velocity of ejecta for this model helps photons to diffuse out more rapidly from the high-velocity edge of the ejecta, and hence, enhances the luminosity.

The absence of the first r-process peak elements in the polar high-velocity region is the other reason for the bright polar emission of model MNS80. As we already pointed out, some of the first-peak r-process elements (Y and Zr) have a significant contribution to the opacity in the optical wavelength. We find that the total Y + Zr mass in the polar region of θしーた ≤ 30° with r/ct ≤ 0.3 is by an order of magnitude smaller for model MNS80 than those for models MNS70a and MNS75a. As a consequence, the gz-band light curves for model MNS80 are significantly brighter than those for other models, while the K-band light curve does not show an appreciable difference among the models.

The kilonova emission for model MNS80 is also much brighter than that of GW170817 with the hypothetical distance of 40 Mpc. In particular, the gz-band magnitudes of model MNS80 for t ≤ 1 d are 2 mag brighter for the polar view (≤20°), although they decline more rapidly after t ≈ 1 d. As is the case for model MNS75a, this suggests that the MHD effect might not play an important role in the post-merger phase of GW170817. On the other hand, if a BNS merger that results in a long-lived MNS with the MHD effect as significant as for model MNS80 occurs, our result indicates that the kilonova can be as bright as in GW170817 even at the distance of 120 Mpc.

4.4. Possible Non-LTE Effect

For the later phase of the kilonova emission (for t ≳ 1 d), the condition of LTE, which we assumed in the present radiative transfer simulations, could break down for the low-density region in which the ionization of the atoms by the radioactive rays becomes more significant than the recombination of ions (Hotokezaka et al. 2021; Kawaguchi et al. 2021). In fact, the importance of taking the non-LTE effect in the excitation/ionization population into account is well known for supernovae radiative transfer simulations (e.g., Boyle et al. 2017).

Hotokezaka et al. (2021) solved the ionization and electron temperature evolution in the kilonova nebular phase using a model in which the atomic properties of the ejecta are represented by Nd ions. They found that the population of the neutral and first ionization atoms are significantly suppressed in the low-density region because of the high efficiency of radioactive ionization due to r-process nuclei. Pognan et al. (2022b) also obtained similar results on the ionization evolution in the nebular phase. These results indicate that the non-LTE effects may suppress the neutral and first ionized ions in the outer part of the ejecta even in the earlier phases. In fact, Pognan et al. (2022a) show that the non-LTE effect can significantly modify the ejecta opacity particularly in the low-density region even at a few days after the merger. Because computing the ionization/excitation population is challenging due to its computational complexity and lack of the atomic data for the r-process elements (see Hotokezaka et al. 2021; Pognan et al. 2022b), we here provide the qualitative estimates for the impacts of the non-LTE effects on the kilonova light curves by removing the neutral or both neutral and first ionized atoms. For this purpose, we perform the radiative transfer simulations with a hypothetical setup in which the neutral or both neutral and the first ionized atoms obtained by solving the Saha's equations are artificially forced to be ionized to the first or the second ionization states, respectively. Note that this prescription is applied to whole ejecta including high-density regions for simplicity.

Figure 13 shows the gzK-band light curves observed from 0° ≤ θしーた ≤ 20° for models in which the opacity contributions from the neutral ("w/o I"; dashed curves) or both neutral and first ionized atoms ("w/o I and II"; dotted curves) are switched off. The light curves with the default setting shown in Figure 9 are also plotted with the solid curves in Figure 13 as the references. By suppressing the opacity contribution from the neutral or the first ionized atoms, the brighter and fainter emission is realized at the time of the peak magnitudes in the gz bands and K band, respectively. Furthermore, the time of the peak magnitudes moves to later and earlier epochs in the g bands and K band, respectively. The effect of suppressing the opacity contribution from the neutral or the first ionized atoms is more significant for the emission in the shorter wavelengths, and the effect on the g band is significant even from t ≈ 0.5 d. On the other hand, the zK-band light curves for t ≲ 1 d are less affected by the suppression of the opacity contribution from the neutral or the first ionized atoms.

The extinction of the polar optical emission, which is seen for models MNS75a and MNS70a (see Figure 9), is suppressed by switching off the opacity contribution from the neutral and the first ionized atoms. The g-band light curves for both model MNS70a and the viscous model become in better agreement with the observed data points of GW170817 if we suppose that the light curves could vary within the ranges up to those calculated without the opacity contribution from both the neutral and the first ionized atoms. On the other hand, Figure 13 indicates that the K-band emission for t ≳ 2 d can become fainter than the observation by the non-LTE effect; though it is still in the range of the uncertainty. This might imply that a different type of BNS merger that produces brighter NIR photons in the kilonova, such as a BNS with unequal-mass NSs or a BNS of which remnant MNS collapses to a BH within a short timescale (≲100 ms), can be more consistent with the observed data of GW170817.

The light curves for models MNS80 and MNS75a show different features than from those observed in GW170817 even if the opacity contribution from the neutral or the first ionized atoms is suppressed: though the difference in the g-band emission from the observation can be less pronounced by the non-LTE effect, the z and K-band emission is still always brighter than the data points by more than 1 magnitude for t ≲ 10 d. These results support our hypothesis that a long-lived MNS associated with a strong global magnetic field is unlikely to be formed in GW170817.

Comparing with the other models, the effect of switching off the opacity contribution from the neutral and the first ionized atoms is relatively minor for model MNS80, especially for the gz-band light curves. This reflects the fact that Y and Zr in the polar region are less abundant for model MNS80 than for the other models, and the suppression of the optical light curves is not significant in the first place.

Figure 13 shows that both the brightness and color of the kilonova emission can be significantly modified if the non-LTE effect plays an important role for determining the ionization population of the ejecta. Beside the non-LTE effect, we note that the uncertainty in the abundance distribution (and thus radioactive heating rate) could also be the source of the systematic error for the light curves because not only the emissivity but also the ejecta temperature distribution, which is responsible for the ionization structure (Barnes et al. 2021), are modified by it. Hence, the systematic study employing different nucleosynthesis models will be needed for quantitatively understanding the uncertainty in the light curve.

5. Conclusions and Discussion

Our results suggest that the presence of the bright synchrotron flare would be an indicator for the presence of long-lived MNSs with significant amplification of the global magnetic field, if we can distinguish it from the jet afterglow. Such synchrotron emission can be observed even for a far distance for which gravitational waves are not detectable. For example the radio afterglow of the MHD models with σしぐま_c ≥ 3 × 10⁷ s⁻¹ and with an ISM density of >10⁻³ cm⁻³ at a distance of 200 Mpc can be identified with untargeted surveys (see, e.g., Hotokezaka et al. 2016a; Dobie et al. 2021). In fact, Dobie et al. (2022) demonstrated that an untargeted search with ASKAP can identify radio transients on timescales of a few days to a year with a flux level down to ∼200 μみゅー Jy if the localization area is reasonably small ∼30 deg². They found one radio transient in the localization area of GW190814, which is likely unrelated to the merger event.

The ASKAP untargeted search for GW190814 provides an interesting upper limit on the surface density of the radio transients above 170 μみゅー Jy varying on timescales of days to a year as < 0.013 deg⁻² (Dobie et al. 2022). If we suppose that a fraction f of BNS mergers with a rate of 300 Gpc⁻¹ yr⁻¹ produces a radio transient similar to model MNS80 with n = 10⁻³ cm⁻³ (see the middle panel of Figure 6), the surface density of such transient is expected to be $\sim 0.03f\,{\deg }^{-2}{({S}_{\nu }/170\,\mu {Jy})}^{-3/2}$ . This observed upper limit implies f ≲ 0.3. Thus, BNS mergers of which remnant MNSs survive for a long period could be the minority, if significant MHD effects are always induced in the remnant MNSs. This is consistent with our latest finding (Fujibayashi et al. 2020c), which shows that, in the presence of long-lived MNSs, relatively light r-process elements should be overproduced, and the entire solar-abundance pattern of the r-process elements, which is believed to be universal, cannot be reproduced.

On the other hand, we emphasize that not necessarily all the BNS mergers produce ejecta with the solar one; though it should be not a major type of BNSs if the solar r -process-abundance pattern is universal. For example, a metal-poor star for which the elemental abundance pattern is different from the solar one might also be explained by the r-process nucleosynthesis in some rare type of a BNS merger, and we do not always have to consider that it was proceeded in a source different from the BNS merger. We also emphasize that a bright long-lasting radio transient has not been found after short gamma-ray bursts (GRBs) even though the extensive follow-up observations have been conducted (Metzger & Bower 2014; Horesh et al. 2016; Fong et al. 2016; Klose et al. 2019; Schroeder et al. 2020). This implies that BNS mergers that result in long-lived MNSs with significant amplification of the global magnetic field (i.e., magnetars) might not be the central engine of short GRBs.

In the presence of a long-lived MNS with a strong global magnetic field, the kilonova emission will also show distinct features from the cases that the magnetic-field amplification and its retention are not significant or the MNS collapses to a BH in a short timescale (≲100 ms); the optical emission is bright for ≲1 d but becomes steeply weak in a short timescale of a few days, and the NIR emission is always brighter by more than 1 magnitude than those observed in GW170817 for t ≲ 10 d. While we should note that the presence of such features, particularly in the optical wavelength, depends on how strongly the non-LTE effects play a role (see Section 4.4), the significant amplification of the global magnetic field can be examined by observing such features in the early epoch. For this purpose, the rapid kilonova search with a high cadence is important.

The bolometric luminosity for the viscous model and MNS70a (except that observed from θしーた ≤ 30°) agrees fairly well with the observed data of GW170817. The peak magnitudes in the grizJHK bands for these models are also in good agreement with the observation of GW170817. The disagreements of the light curves after the time of the peak magnitudes will be much less pronounced if we consider the possible uncertainty of the light-curve prediction due to the non-LTE effects. On the other hand, we showed that the kilonova light curves for the MHD models with σしぐま_c ≥ 3 × 10⁷ s⁻¹ will exhibit different features from GW170817. This suggests that the MHD effect might not play an important role in the post-merger phase of GW170817. Since the significance of the MHD effect depends strongly on the parameters of the phenomenological dynamo model, it is not clear whether and, if so, at which time the remnant MNS collapsed to a BH. Thus, the investigation for the realistic value for the phenomenological dynamo parameters as well as the quantitative evaluation of the non-LTE effects will be important tasks to interpret the systems of GW170817 and the future events.

K.K. thanks Smaranika Banerjee, Nanae Domoto, and Masaomi Tanaka for the valuable discussions. Numerical computation was performed on Yukawa21 at Yukawa Institute for Theoretical Physics, Kyoto University and the Sakura, Cobra, Raven clusters at Max Planck Computing and Data Facility. This work was supported by Grant-in-Aid for Scientific Research (JP20H00158, JP21K13912) of JSPS/MEXT and Japan Society for the Promotion of Science (JSPS) Early-Career Scientists grant No. 20K14513.

Appendix A: Treatment of Radioactive Heating

In this paper, we update the treatment of radioactive heating in the HD simulation from the previous study, based on the formulation of Takahashi et al. (2016), Uchida et al. (2017). In the previous study (Kawaguchi et al. 2021), the rest-mass density, ρろー_*, was evolved as a conserved quantity in solving the continuity equation. However, strictly speaking, using ρろー_* as the conserved quantity is incorrect in the presence of the radioactive reaction because the rest-mass of baryons can be changed due to the change of the nuclear binding energy. Thus, in this work, we define $\hat{\rho }={m}_{{\rm{u}}}{n}_{{\rm{B}}}$ where m_u is the atomic mass unit and n_B is the baryon number density, and define ${\hat{\rho }}_{* }=\hat{\rho }{u}^{t}\sqrt{g}$ as a conserved quantity, which follows the continuity equation without a source term; that is,

$\begin{eqnarray}&&{\partial }_{t}{\hat{\rho }}_{* }+{\partial }_{i}\left({\hat{\rho }}_{* }{v}^{i}\right)=0,\end{eqnarray} \tag{ A1 }$

where vⁱ = uⁱ/u^t.

The release of the nuclear binding energy by the radioactive decay is described by the change in the mass per baryon, m_b: we define the mass excess of the baryon by

$\begin{eqnarray}&&{\delta }_{\mathrm{exc}}:=\displaystyle \frac{{m}_{{\rm{b}}}-{m}_{{\rm{u}}}}{{m}_{{\rm{u}}}}.\end{eqnarray} \tag{ A2 }$

Using δでるた_exc, the rest-mass density is given by $\rho =(1+{\delta }_{\mathrm{exc}})\hat{\rho }$ . The enthalpy per baryon particle (weighted by 1/m_u), $\hat{h}$ , is then written as

$\begin{eqnarray}&&\hat{h}=\left(1+{\delta }_{\mathrm{exc}}\right){c}^{2}+\hat{\epsilon }+\displaystyle \frac{P}{\hat{\rho }},\end{eqnarray} \tag{ A3 }$

with $\hat{\epsilon }$ being the internal energy per baryon particle (weighted by 1/m_u). Note that $\hat{h}$ is related to the specific enthalpy (enthalpy per rest mass), h, defined in the previous paper with $\rho h=\hat{\rho }\hat{h}$ . Then, the energy-momentum tensor is written as

$\begin{eqnarray}&&{T}_{\mu \nu }=\hat{\rho }\hat{h}{u}_{\mu }{u}_{\nu }+{{Pg}}_{\mu \nu }.\end{eqnarray} \tag{ A4 }$

Considering that a fraction of the energy release by the radioactive decay is carried and lost from ejecta by the neutrino emission, the energy-momentum conservation is written as

$\begin{eqnarray}&&{{\rm{\nabla }}}_{\mu }{T}^{\mu \alpha }=-\hat{\rho }{\dot{q}}_{(\nu )}{u}^{\alpha }(\tau ),\end{eqnarray} \tag{ A5 }$

where ${\dot{q}}_{(\nu )}$ denotes the neutrino energy deposition rate per baryon particle (weighted by 1/m_u) due to the radioactive decay, and τたう denotes the proper time (affine parameter) in the fluid rest frame. The time evolution of the mass excess in the fluid rest frame is connected to the radioactive-heating rate per baryon particle (weighted by 1/m_u), ${\dot{q}}_{\mathrm{tot}}$ , with

$\begin{eqnarray}&&\displaystyle \frac{d{\delta }_{\mathrm{exc}}}{d\tau }=-\displaystyle \frac{1}{{c}^{2}}{\dot{q}}_{\mathrm{tot}}(\tau ).\end{eqnarray} \tag{ A6 }$

We can rewrite this equation into a form that is similar to the conservation form of the hydrodynamics;

$\begin{eqnarray}&&{{\rm{\nabla }}}_{\mu }\left(\hat{\rho }{\delta }_{\mathrm{exc}}{u}^{\mu }\right)=-\displaystyle \frac{1}{{c}^{2}}\hat{\rho }{\dot{q}}_{\mathrm{tot}}(\tau ).\end{eqnarray} \tag{ A7 }$

To evaluate ${\dot{q}}_{(\nu )}$ and ${\dot{q}}_{\mathrm{tot}}$ at each time step, the proper time in the fluid rest frame is needed. For this purpose, we also determine the proper time of each fluid element by solving

$\begin{eqnarray}&&{{\rm{\nabla }}}_{\mu }\left(\hat{\rho }\tau {u}^{\mu }\right)=\hat{\rho },\end{eqnarray} \tag{ A8 }$

which is equivalent to solving d τたう/d τたう = u^μみゅー∇_μみゅー τたう = 1.

Appendix B: Numerical Convergence of the Ejecta Kinetic Energy Distribution

Figure 14 compares the kinetic energy distribution for MNS75a obtained by the HD simulations with various grid resolutions. We found that the model with N_r = 3072 and N_θしーた = 192 is fairly sufficient to resolve the ejecta fast tail with u^r/c ≲ 5, where N_r and N_θしーた denote the numbers of the grid points in r and θしーた directions, respectively.

We also confirmed that the light curves of the kinetic energy distribution are not affected by the presence of the radioactive heating terms in the HD simulations (see the curves with "w/ heating" in Figure 14). This can be understood by the fact that the energy deposition due to radioactive heating is negligible compared to the kinetic energy for the fast tail.

Electromagnetic Counterparts of Binary-neutron-star Mergers Leading to a Strongly Magnetized Long-lived Remnant Neutron Star

Article metrics

Author e-mails

Author affiliations

ORCID iDs

Dates

Abstract

1. Introduction

2. Method