AFTERGLOW OF A BINARY NEUTRON STAR MERGER

Coalescence of binary neutron stars (BNSs) is one of the most promising sources for next-generation kilometer size gravitational-wave detectors such as the advanced Laser Interferometer Gravitational-wave Observatory, advanced VIRGO, and the Large-scale Cryogenic Gravitational Wave Telescope. Statistical studies have suggested that the detection rate of gravitational waves emitted by BNSs will be ∼1–100 year⁻¹ (Kalogera et al. 2007). The typical signal-to-noise ratio for such detection will be ∼10 or less. Thus, it will be crucial for the detection of gravitational waves to find electromagnetic counterparts to the gravitational-wave signals. Short–hard gamma-ray bursts (SGRBs) have been inferred to accompany the BNS merger (Narayan et al. 1992; Piran 2004). However, this hypothesis relies on many uncertain assumptions, e.g., high magnetic-field strength or efficient pair annihilation of neutrino–antineutrino. In this Letter, we give a conservative estimate for the strength of electromagnetic signals based on a numerical-relativity simulation and show that a strong electromagnetic signal will indeed accompany with the BNS merger.

A BNS evolves due to gravitational radiation reaction and eventually merges. After the merger sets in, there are two possible fates (e.g., Shibata et al. 2005; Kiuchi et al. 2009): if the total mass M is larger than a critical mass M_c, a black hole will be formed, while a hypermassive neutron star (HMNS) will be formed for M < M_c. The value of M_c depends strongly on the equation of state (EOS) of neutron stars, but the latest discovery of a high-mass neutron star with mass 1.97  ±  0.04 M_☉ (Demorest et al. 2010) indicates that the EOS is stiff and M_c may be larger than the typical total mass of the BNS ∼2.7 M_☉ (Stairs 2004). This indicates that HMNS is the likely outcome for many BNS mergers, at least temporarily (Hotokezaka et al. 2011).

Neutron stars in nature have a strong magnetic field with a typical field strength at the stellar surface of 10¹¹–10¹³ G. One of the neutron stars in a BNS often has field strength smaller than this typical value, ∼10¹⁰–10¹¹ G (Lorimer 2008), probably because of the accretion history of the first-formed neutron star during the formation of the second one. However, at least the second one is likely to have the typical magnetic-field strength.

It is reasonable to believe that each neutron star in the inspiral phase (before the merger sets in) has an approximately dipole magnetic field as in the isolated one. During the late inspiral phase and formation of an HMNS in the merger phase, its magnetic-field profile will be modified due to magnetohydrodynamic (MHD) processes. However, in zeroth approximation, it would be safe to suppose that the dipole field is dominant. For this reason, we consider the evolution of an HMNS with dipole magnetic fields in the following.

One of the most important properties of HMNSs is that they are rapidly and differentially rotating (Shibata et al. 2005). The numerical simulations have shown that the typical angular velocity at its center is Ω ≳ 10⁴ rad s⁻¹, much larger than that of ordinary pulsars, while at equatorial surface it is ∼10³ rad s⁻¹.³ Because of the presence of the differential rotation, the winding of magnetic fields is enhanced: toroidal magnetic-field strength B_T in HMNSs increases linearly with time (t) in the presence of seed poloidal (cylindrically radial) magnetic fields B_P (B_T increases as ∼B_PΩt). The increase of the magnetic-field strength results in the increase of magnetic pressure. Because only dilute matter is present in the surface of HMNSs, Alfvén waves are likely to propagate near the rotational axis with ϖ ≲ 10 km, where $\varpi =\sqrt{x^2+y^2}$, transporting electromagnetic energy generated in the HMNS along the rotational axis; tower-type outflow is driven. As far as the HMNS is active and the rapid rotation is present, the amplification of the toroidal magnetic field continues via the winding effect. Then, the electromagnetic energy should increase approximately as $\dot{E}_B \sim B_P^2 V \Omega$ as described in Meier (1999), where V is an effective volume in which the amplification occurs.⁴ If the amplification efficiently occurs near the rotation axis with ϖ ≲ 10 km, V is approximately ∼(αR)³, where R is the equatorial stellar radius ∼10 km and α is a constant of O(0.1). For a typical HMNS formed after the merger

$$\begin{eqnarray} \dot{E}_B &\sim &10^{46} B_{13}^2 R_6^3 \Omega _4 \alpha _{0.1}^3 \,{\rm erg \,s^{-1}}, \end{eqnarray} \tag{ 1 }$$

where B₁₃ = B_P/10¹³ G, R₆ = R/10 km, Ω₄ = Ω/10⁴ rad s⁻¹, and α_0.1 = α/0.1. We suppose a relatively high magnetic-field strength because it is likely to be amplified by compression that occurs during the merger (Rezzolla et al. 2011). Thus, the luminosity of the electromagnetic radiation will be as high as that of quasars for a typical field strength of a progenitor neutron star, B_P ∼ 10¹² G, and of a resulting HMNS, B_P ∼ 10¹³ G. If a substantial fraction of this generated electromagnetic energy is converted to electromagnetic radiation (as suggested, e.g., in Nakar & Piran 2011), the merger event may be detected by telescopes as electromagnetic signals.

Motivated by the fact mentioned above, we performed MHD simulation for an HMNS in general relativity. As the initial condition, we prepare a rapidly and differentially rotating HMNS in axisymmetric equilibrium as in Shibata et al. (2006): we constructed an HMNS model with the following piecewise polytropic EOS: $P=P_{\rm cold}=K_1\rho ^{\Gamma _1}$ for ρ ⩽ ρnuc and $P_{\rm cold}=K_2\rho ^{\Gamma _2}$ for ρ ⩾ ρnuc. Here, P and ρ are the pressure and rest-mass density, respectively. We set Γ1 = 1.3, Γ2 = 2.75, K1 = 5.16 × 1014 cgs, K2 = K1ρΓ1 − Γ2nuc, and ρnuc = 1.8 × 1014 g cm−3. With this EOS, the maximum gravitational mass (rest mass) is 2.01 M☉ (2.32 M☉) for spherical neutron stars and 2.27 M☉ (2.60 M☉) for rigidly rotating neutron stars. These are similar values to those in realistic stiff EOS (e.g., Read et al. 2009 for a review). We prepare an HMNS with the following physical parameters: gravitational mass M = 2.65 M☉, baryon rest mass Mb = 2.96 M☉, maximum density ρmax = 9.0 × 1014 g cm−3, angular momentum J = 0.82GM2/c, ratio of polar to equatorial radius 0.3, central rotation period Pc = 0.202 ms, and rotation period at the equatorial surface 5.4Pc. Here, G is the gravitational constant. The rotation law is specified in the same way as in Baumgarte et al. (2000) with the differential rotation parameter $\hat{A}=0.8$. This HMNS is similar to that found in the BNS merger simulation of Shibata et al. (2005) and Kiuchi et al. (2009) performed with a nuclear-theory-based EOS. In the evolution, we employ an EOS of the form P = Pcold + (Γ1 − 1)ρ(ε − εcold), where ε is the specific internal energy and εcold is its cold part determined from Pcold.

Poloidal magnetic fields, for which the toroidal component of the vector potential has the form A_φ = A₀ϖ₀ϖ²/(r² + ϖ²₀)^3/2, are superimposed on this initial condition. The magnetic field in the inertial frame is given by Bⁱ = ^ijk∂_jA_k, where _ijk is the completely antisymmetric tensor. Here, r² = ϖ² + z² and (ϖ₀, A₀) are constants: ϖ₀ is chosen to be 5R_e/3–20R_e/3, where R_e(= 12.1 km) is the coordinate radius on the equatorial plane. We found that the electromagnetic luminosity shown below depends weakly on this parameter. In the following, we show the results for ϖ₀ = 10R_e/3 and 20R_e/3. A₀ determines the field strength for which we give the maximum magnetic-field strength B_max ≈ 10¹³–4 × 10¹⁴ G. Here, the magnetic-field strength is defined by $B=\sqrt{b_{\mu }b^{\mu }}$, where b^μ is the 4-vector of the magnetic field in the frame comoving with the fluid. With such strength, the magnetic pressure in the HMNS is much smaller than the matter pressure at its center, and thus, the density profile of the HMNS (except for its surface) is not significantly modified by the magnetic-field effect. Because of the presence of differential rotation, the magnetic fields may be amplified in an exponential manner due to magnetorotational instability (MRI; Balbus & Hawley 1998). However, the wavelength for the fastest growing mode is very short (∼10²–10⁴ cm) in the present setting and it is not possible to resolve it in the numerical simulation. Here we do not pay attention to MRI but only to the winding effect. In the presence of the MRI effects, the electromagnetic energy will be increased more rapidly and the luminosity of electromagnetic radiation may even be enhanced. Thus, this work would determine the lower bound of the magnetic luminosity, which is, however, quite high.

The MHD simulation is performed assuming that the ideal MHD condition holds. A conservative shock capturing scheme is employed for solving MHD equations as in Shibata & Sekiguchi (2005): in the present work, a numerical scheme with third-order accuracy in space and fourth-order accuracy in time is employed. Einstein's evolution equations are solved in fourth-order accuracy in space and time in the so-called BSSN-puncture formulation (Shibata & Nakamura 1995; Baumgarte & Shapiro 1999; Campanelli et al. 2006).

Any conservative scheme in MHD cannot handle vacuum, and hence, we have to add an atmosphere of small density outside the HMNS. Because the matter outside an HMNS formed in a real merger would be dilute, the density of the atmosphere should be as small as possible to exclude spurious effects by it. We set the density of the atmosphere as

$$\begin{eqnarray} \rho _{\rm at}=\left \lbrace \begin{array}{@{}ll}f_{\rm at} \rho _{\rm max} & r < 2R_e, \\ \ms f_{\rm at} \rho _{\rm max}(r/2R_e)^{-n} & r \ge 2R_e, \\ \end{array}\right. \end{eqnarray} \tag{ 2 }$$

where we choose n = 2 or 2.5. f_at is constant, for which we typically give 10⁻⁹. We changed the values of f_at from 10⁻¹⁰ to 10⁻⁷ for B_max = 4.2 × 10¹³ G and found that as long as B²/(4πρ_atc²) ≲ 1, our code works well. For a large value of f_at, the evolution of magnetic fields is substantially affected by the inertia of the matter. However, with decreasing the value of f_at to B²/(4πρ_atc²) ∼ 1, the dependence of magnetic-field evolution on the atmosphere density becomes weak, and hence, the effect of the artificial atmosphere does not play a role. The velocity of the atmosphere is initially set to be zero. With this treatment, the magnetic field is initially modified for t ≲ P_c. However, such modification plays a minor role after the winding effect becomes dominant for the magnetic-field amplification. We always perform simulations for a time much longer than P_c and focus on the stage for which a quasisteady state is achieved. Thus, the artificial effect associated with the initial setting does not matter.

Axisymmetric numerical simulation is performed in cylindrical coordinates for MHD and in Cartesian coordinates for Einstein's equation part (using the so-called Cartoon method). The details are described in Shibata & Sekiguchi (2005) and the references therein. The nonuniform grid is prepared as in Kiuchi et al. (2008). We here impose axial symmetry to guarantee a sufficiently high grid resolution although we can perform a nonaxisymmetric simulation. In the present setting, the equatorial coordinate radius of the HMNS, R_e, is covered by 150 uniform grids. Smaller grids with 100 and 120 were also adopted to check the convergence of the numerical results. Outer boundaries along the x- and z-axes are located at ≈170 R_e ≈ 2000 km.

Figure 1 plots the evolution of the density profiles and magnetic-field strength outside the HMNS. Numerical simulation shows that the system evolves in the following manner. Because of the presence of differential rotation, the winding of magnetic fields proceeds. Toroidal field strength is increased linearly with time, and the growth rate is high in particular near the rotational axis (ϖ ≲ 10 km). Hence, Alfvén waves propagate primarily toward the z-direction near the rotational axis, and magnetic-field strength there also increases. After the substantial winding, the magnetic pressure B²/8π becomes larger than the gravitational potential energy density near the polar surface, ∼GMρ/H, where H(= 0.3R_e) is the vertical coordinate radius of the HMNS. Then, the matter of the HMNS in the vicinity of its polar surface is stripped, and an outflow is driven. Because the magnetic energy density B²/4π is comparable to or slightly smaller than the rest-mass density, ∼(H/GM)(B²/8π) > B²/4πc², the outflow is mildly relativistic with the velocity of order 0.1c in the vicinity of the HMNS. However, in the region far from the HMNS, the outflow velocity could be fairly relativistic, ∼0.9c, near the rotation axis (see Figure 3).

Zoom In Zoom Out Reset image size

Figure 2 plots evolution of the matter and electromagnetic energy ejection rates $\dot{E}_M$ and $\dot{E}_B$, respectively, defined by

$$\begin{eqnarray} &&\dot{E}=-\oint _{r={\rm const}} d\Omega \sqrt{-g} T^r_{ t}, \end{eqnarray} \tag{ 3 }$$

where we substitute the stress energy tensor for the matter and electromagnetic fields into T^r_t, respectively. g is the determinant of the spacetime metric. The surface integral is performed for r ≈ 480 km. We checked that the luminosity depends only weakly on the radius of the surface integral. The top-left and top-right panels show the evolution of $\dot{E}_M$ and $\dot{E}_B$ for B_max = 4.2 × 10¹³ G, and the bottom-left panel shows the evolution of $\dot{E}_B$ for B_max = 1.7 × 10¹³–4.2 × 10¹⁴ G. The bottom-right panel shows the ejection rates integrated only for the angle Δθ = π/20 for B_max = 4.2 × 10¹³ G. The top-left panel is for ϖ₀ = 10R_e/3, and n = 2 with three grid resolutions, and the top-right panel is for different values of ϖ₀ and n. The top-left and top-right panels show that irrespective of the grid resolution, ϖ₀, and n, the total ejection rates are

$$\begin{eqnarray} \dot{E}_M &\sim & 10^{48} B_{13}^2 R_6^3 \Omega _4 \,{\rm erg \,s^{-1}}, \end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray} \dot{E}_B &\sim & 10^{47} B_{13}^2 R_6^3 \Omega _4 \,{\rm erg \,s^{-1}}. \end{eqnarray} \tag{ 5 }$$

The bottom-left panel indeed shows that the scaling relation with respect to the magnetic-field strength holds (the same scaling also holds for $\dot{E}_M$). Another important point is that these ejection rates do not significantly vary in time. Thus, a quasisteady outflow is driven.

Zoom In Zoom Out Reset image size

The value of $\dot{E}_B$ agrees approximately with the prediction of Equation (1), implying that the scenario described in Section 1 is correct. The ratio of $\dot{E}_M/\dot{E}_B$ is of order 2c²H/GM ∼ 10. This agrees approximately with the value required for the mass stripping.

Comparison among top-left, top-right, and bottom-right panels of Figure 2 shows that the electromagnetic energy is mainly emitted in the direction near the rotation axis. By contrast, the matter energy is emitted in a fairly isotropic manner. Along the rotational axis, $\dot{E}_M/\dot{E}_B$ is of order unity (∼1–10). This is the reason that the outflow along the rotation axis can be mildly relativistic.

The amount of angular momentum loss by the matter ejection and electromagnetic radiation in the time duration Δt = P_c is much smaller than the total angular momentum because of our choice of magnetic-field strength, ≲ 3 × 10¹⁴ G. This implies that the matter and electromagnetic waves are continuously ejected, and a quasisteady outflow is formed. Figure 3 plots the density and velocity profiles along the rotation axis for B_max = 4.2 × 10¹³ G, ϖ₀ = 10R_e/3, and n = 2. The density profile (left panel) indeed shows that the averaged density does not change significantly with time. The density decreases with the radius. In these examples, the power-law index is roughly n_ρ ∼ 1.5–2 ($\rho \propto r^{-n_{\rho }}$), but this number depends on the initial setting for B_max and ϖ₀ and varies with time. The velocity profile (right panel of Figure 3) shows that the outflow is mildly relativistic. The profile varies in a short timescale. The maximum velocity is ∼0.9c as mentioned above: this maximum depends weakly on the setting of the atmosphere; for the lower atmosphere density, the maximum speed is larger. The averaged magnitude of the outflow velocity in time is ∼0.4–0.6c (which also depends weakly on the setting of the atmosphere). Because the corresponding Lorentz factor of the jet is Γ_j ≲ 2, the relativistic beaming effect (i.e., observable viewing angle is ≲ 1/Γ_j) is relatively small, and thus, this jet may be observable from a large solid angle.

Zoom In Zoom Out Reset image size

The centrifugal force due to rapid and differential rotation plays a crucial role for supporting strong self-gravity of HMNSs. This suggests that any HMNS will eventually collapse to a black hole after a substantial loss of angular momentum by gravitational-wave emission (Shibata et al. 2005) and/or after a substantial angular momentum transport inside it due to magnetic viscous effects (Duez et al. 2006) in a realistic situation. The predicted lifetime is of order 10–100 ms. After the collapse to a black hole, a system composed of a black hole and compact accretion disk will be formed (Shibata et al. 2006; Duez et al. 2006; Rezzolla et al. 2011; Hotokezaka et al. 2011). Then, the electromagnetic radiation is likely to be emitted through the Blandford–Znajek mechanism (Blandford & Znajek 1977; McKinney & Gammie 2004; McKinney 2005). However, the lifetime of the accretion disk will not be longer than ∼100 ms because of the viscous evolution. Hence, the magnetic energy generation and mass ejection are likely to continue only for a short timescale of order 100 ms after the BNS merger.

However, the luminosity is quite high. If the kinetic energy of the matter and/or electromagnetic energy are efficiently converted to an electromagnetic signal, it will be observed even if the merger happens at a distance of several hundred Mpc as discussed, e.g., in Nakar & Piran (2011). Strong radio afterglow emission could be expected when the jet propagates in the matter that may be ejected from neutron stars and HMNSs during the merger. The ejected mass that can be estimated as $\dot{E}_M\Delta T/c^2 \sim 10^{-6}\,M_\odot B_{13}^2 R_6^3 \Omega _4 (\Delta T/0.1\, {\rm s})$ could yield radioactive elements and be observed like dim supernovae (Li & Paczyński 1998).

In this Letter, we focus only on the HMNS with conservative magnetic-field strength B ∼ 10¹³ G (achieved by compression of ordinary field strength ∼10¹² G). If the magnetic-field strength were as high as that of magnetar (Woods & Thompson 2006), i.e., B ∼ 10¹⁵ G, the electromagnetic luminosity would reach 10⁵¹ erg s⁻¹ (see the bottom-left panel of Figure 2). This value is as high as the luminosity of GRBs, and because the expected time duration is less than 1 s, this is also a candidate model for SGRB (Nakar 2007). The canonical peak isotropic luminosity of an SGRB is ∼10⁵¹ erg s⁻¹, which is consistent with this estimation, assuming the jet solid angle as Ω_j ∼ 0.1 and the conversion efficiency from the jet kinetic energy into the gamma rays as η ∼ 0.1. Hence, if one of the neutron stars in a BNS has a large magnetic-field strength or the magnetic-field strength is significantly amplified during the merger process or in the formed HMNS, this may be observed as an SGRB.

This work was supported by Grant-in-Aid for Scientific Research (nos. 19047004, 21340051, 21684014, 22244019, 22244030), by Grant-in-Aid for Young Scientists (B) (no. 22740178), and by Grant-in-Aid for Scientific Research on Innovative Area (no. 20105004), and HPCI Strategic Program of Japanese MEXT.