MAGNETIC ENERGY PRODUCTION BY TURBULENCE IN BINARY NEUTRON STAR MERGERS

The inspiral and coalescence of binary neutron star (BNS) systems is a topic of increasingly intensive research in observational and theoretical astrophysics. It is anticipated that the first direct detections of gravitational wave (GW) will be from compact binary mergers. BNS mergers are also thought to produce short-hard gamma-ray bursts (SGRBs; Eichler et al. 1989; Narayan et al. 1992; Belczynski et al. 2006; Metzger et al. 2008). Simultaneous detections of a prompt GW signal with a spatially coincident electromagnetic (EM) counterpart dramatically increase the potential science return of the discovery. For this reason, there has been considerable interest as to which, if any, detectable EM signature may result from the merger (Metzger & Berger 2012; Piran et al. 2013). Other than SGRBs and their afterglows, including those viewed off-axis (Rhoads 1997; van Eerten & MacFadyen 2011), suggestions include optical afterglows associated with the radioactive decay of tidally expelled r-process material (Li & Paczyński 1998; Metzger et al. 2010; though detailed calculations indicate that they are faint (Barnes & Kasen 2013)), radio afterglows following the interaction of a mildly relativistic shell with the interstellar medium (Nakar & Piran 2011), and high-energy pre-merger emission from resistive magnetosphere interactions (Hansen & Lyutikov 2001; Piro 2012).

Merging neutron stars possess abundant orbital kinetic energy (∼10⁵³ erg). A fraction of this energy is certain to be channeled through a turbulent cascade triggered by hydrodynamical instabilities during the merger. Turbulence is known to amplify magnetic fields by stretching and folding embedded field lines in a process known as the small-scale turbulent dynamo (Vaĭnshteĭn & Zel'dovich 1972; Chertkov et al. 1999; Tobias et al. 2011). Amplification stops when the magnetic energy grows to equipartition with the energy containing turbulent eddies. An order of magnitude estimate of the magnetic energy available at saturation of the dynamo can be informed by global merger simulations. These studies indicate the presence of turbulence following the nonlinear saturation of the Kelvin–Helmholtz instability activated by shearing at the NS surface layers (Price & Rosswog 2006; Liu et al. 2008; Anderson et al. 2008; Rezzolla et al. 2011; Giacomazzo et al. 2011). The largest eddies produced are on the ∼1 km scale and rotate at ∼0.1c, setting the cascade time t_eddy and kinetic energy injection rate ε_K ≡ E_K/t_eddy at 0.1 ms and 5 × 10⁵⁰ erg s⁻¹, respectively. When kinetic equipartition is reached, each turbulent eddy contains 5 × 10⁴⁶ erg of magnetic energy, and a mean magnetic field strength

$$\begin{equation} B_\mathrm{RMS} \gtrsim 10^{16}\, \mathrm{G} \left(\frac{\rho }{10^{13}\,\mathrm{{\rm g\,cm^{-3}}}} \right)^{1/2} \left(\frac{v_\mathrm{eddy}}{0.1c} \right). \end{equation} \tag{ 1 }$$

Whether such conditions are realized in merging neutron star systems depends upon the dynamo saturation time t_sat and equipartition level E_sat/E_K. In particular, if t_sat ≲ t_merge then turbulent volumes of neutron star material will contain magnetar-level fields throughout the early merger phase. Once saturation is reached, a substantial fraction of the injected kinetic energy, 0.7ε_K, is resistively dissipated (Haugen et al. 2004) at small scales. Magnetic energy dissipated by reconnection in optically thin surface layers will accelerate relativistic electrons (Blandford & Eichler 1987; Lyutikov & Uzdensky 2003), potentially yielding an observable EM counterpart, independently of whether the merger eventually forms a relativistic outflow capable of powering a short gamma-ray burst (GRB).

In this Letter, we demonstrate that the small-scale turbulent dynamo saturates quickly, on a time t_sat ≲ t_merge, and that ≳ 10¹⁶ G magnetic fields are present throughout the early merger phase. This implies that the magnetic energy budget of merging BNSs is controlled by the rate with which hydrodynamical instabilities randomize the orbital kinetic energy. Our results are derived from simulations of the small-scale turbulent dynamo operating in the high-density, trans-relativistic, and highly conductive material present in merging neutron stars. We have carefully examined the approach to numerical convergence and report grid resolution criteria sufficient to resolve aspects of the small-scale dynamo. Our Letter is organized as follows. The numerical setup is briefly described in Section 2. Section 3 reports the resolution criterion for numerical convergence of the dynamo completion time and the saturated field strength. In Section 4, we assess the possibility that magnetic reconnection events may convert a sufficiently large fraction of the magnetic energy into high-energy photons to yield a prompt EM counterpart detectable by high-energy observatories including Swift and Fermi.

The equations of ideal relativistic magnetohydrodynamics (MHD) have been solved on the periodic unit cube with resolutions between 16³ and 1024³:

$$\begin{eqnarray} \nabla _\nu N^\nu &= 0 \end{eqnarray} \tag{ 2a }$$

$$\begin{eqnarray} \nabla _\nu T^{\mu \nu } &= S^\mu \end{eqnarray} \tag{ 2b }$$

$$\begin{eqnarray} \frac{\partial \mathbf {B} }{\partial t} &= \mathbf {\nabla } \times (\mathbf {v} \times \mathbf {B}). \end{eqnarray} \tag{ 2c }$$

Here, $T^{\mu \nu} = ph^{*} u^{\mu } u^{\nu } + p^{*} g^{\mu \nu } - b^{\mu } b^{\nu }, b^{\mu } = F^{\mu }_{\nu } u^{\nu }$ is the magnetic field four-vector, and h* = 1 + e* + p*/ρ is the total specific enthalpy, where p* = p_g + b²/2 is the total pressure, p_g is the gas pressure, and e* = e + b²/2ρ is the specific internal energy. The source term S^μ = ρa^μ − ρ(T/T₀)⁴u^μ includes injection of energy and momentum at the large scales and the subtraction of internal energy (with parameter T₀ = 40 MeV) to permit stationary evolution. Vortical modes at k/2π ⩽ 3 are forced by the four-acceleration field a^μ = du^μ/dτ which smoothly decorrelates over a large-eddy turnover time, as described in Zrake & MacFadyen (2012, 2013).

We have employed a realistic micro-physical equation of state (EOS) appropriate for the conditions of merging neutron stars. It includes contributions from high-density nucleons according to a relativistic mean field model (Shen et al. 1998, 2010), a relativistic and degenerate electron–positron component, neutrino and anti-neutrino pairs in beta equilibrium with the nucleons, and radiation pressure. For our conditions, all the components make comparable contributions to the pressure. We have also employed a simpler gamma-law EOS and found close agreement for the conditions explored in this Letter, indicating that the MHD effects are insensitive to the EOS for transonic conditions. The models presented in our resolution study use the far less expensive gamma-law EOS.

All of the simulations presented in this study use the HLLD approximate Riemann solver (Mignone et al. 2009), which has been demonstrated as crucial in providing the correct spatial distribution of magnetic energy in MHD turbulence (Beckwith & Stone 2011). The solution is advanced with an unsplit, second-order MUSCL-Hancock scheme. Spatial reconstruction is accomplished with the piecewise linear method configured to yield the smallest possible degree of numerical dissipation. The divergence constraint on the magnetic field is maintained to machine precision at cell corners using the finite volume CT method of Tóth (2002). Full details of the numerical scheme may be found in Zrake & MacFadyen (2012).

3.1. Growth and Saturation of the Magnetic Energy

3.2. Numerical Convergence

The same driven turbulence model was run through magnetic saturation at resolutions 16³, 32³, 64³, 128³, 256³, and 512³. Another model at 1024³ was run through the end of the kinematic phase, but further evolution was computationally prohibitive. Figure 1 shows the development of B_RMS, E_M, and E_K as a function of time at each resolution.

Zoom In Zoom Out Reset image size

We find that sufficiently resolved runs (⩾512³) attain mean magnetic field strengths of 10¹⁶ G within two large-eddy rotations. All models with resolutions ⩾32³ eventually attain mean fields of ≳ 3 × 10¹⁶ G. The saturated field strength increases until resolution 256³. We find that the kinematic growth rate is higher at each higher resolution, while the timescale for the nonlinear saturation converges at 256³ to roughly five large-eddy rotation times, or about 0.5 ms for the physical parameters of BNS mergers.

In order to quantitatively describe the time development of magnetic energy E_M(t), we describe it with an empirical model,

$$\begin{equation} E_M(t) =\left\lbrace \begin{array}{@{}l@{\quad }l@{}}E_M(0) e^{t^2/\tau ^2_0} & 0 < t < t_1 \\ E_M(t_1) e^{(t-t_1)/\tau _1} & t_1 < t < t_\mathrm{nl} \\ E_\mathrm{nl} + (E_\mathrm{sat} - E_\mathrm{nl}) (1 - e^{-(t-t_\mathrm{nl})/\tau _2}) & t_\mathrm{nl} < t < \infty \end{array}\right., \end{equation} \tag{ 3 }$$

where the six parameters (summarized in Table 1) are obtained by a least-squares optimization. Figure 2 shows the empirical model given in Equation (3) applied to a representative run at 128³. The first phase, t < t₁, is a startup transient and lasts until the hydrodynamic cascade is fully developed at t₁ ≈ t_eddy. The kinematic dynamo phase lasts between t₁ and t_nl, during which the magnetic energy exponentiates on the timescale τ₁. At t_nl, the smallest scales reach kinetic equipartition and the growth rate slows. In the final stage, E_M asymptotically approaches E_sat on the timescale τ₂. We define the dynamo completion time t_sat as t_nl + τ₂. Figure 3 shows the best-fit E_sat, t_sat, and τ₁ as a function of the resolution.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1. Empirical Model for Magnetic Energy Growth

Fit Parameter	Description	Numerically Converged Value
τ0	Startup timescale	None, artifact
t1	Hydrodynamic cascade fully developed	teddy
τ1	Kinematic growth timescale	None, ∝N−2/3
tnl	End of kinematic phase	2teddy
τ2	Nonlinear saturation timescale	tnl + τ2 ≈ 5teddy, 2563
Esat	Saturated magnetic energy	0.6 EK, 5123

Notes. Model parameters characterizing the growth of magnetic energy before and during the dynamo saturation. Numerical convergence is attained for the fully saturated magnetic energy E_sat and the dynamo completion time t_sat defined as t_nl + τ₂.

Download table as:  ASCIITypeset image

The magnetic energy E_sat at dynamo completion shows signs of converging to a value of 0.6 E_K by resolution 512³. The timescale τ₂ on which the magnetic energy asymptotically approaches E_sat is consistently ≈3t_eddy at different resolutions. The dynamo completion time t_sat is numerically converged at ≈5t_eddy by 256³. The best-fit kinematic growth time follows a power law in the resolution, τ₁∝N^−0.6. This is consistent with the value of −2/3 expected if the dynamo time is controlled by shearing at the smallest scale, the cascade is Kolmogorov (i.e., u_ℓ∝ℓ^1/3), and the viscous cutoff ℓ_ν occurs at a fixed number of grid zones. In that case, $\tau _1 \sim t_\nu \propto \ell _\nu ^{2/3} \sim N^{-2/3}$.

3.3. Power Spectrum of Kinetic and Magnetic Energy

The time development of kinetic and magnetic energy power spectra has been studied for a single run with resolution 512³. We present three-dimensional, spherically integrated power spectra with the dimensions of erg cm⁻³ m defined as

$$\begin{eqnarray} P_K(k_i) &= \frac{1}{\Delta k_i}\sum _{ \mathbf {k} \in \Delta k_i} {\left|\right. \mathcal {F}_{ \mathbf {k} }\left[ \mathbf {v} \sqrt{\rho /2}\right]\left|\right. ^2} \end{eqnarray} \tag{ 4a }$$

$$\begin{eqnarray} P_M(k_i) &= \frac{1}{\Delta k_i}\sum _{ \mathbf {k} \in \Delta k_i} {\left|\right. \mathcal {F}_{ \mathbf {k} }\left[ \mathbf {B} /\sqrt{8\pi }\right]\left|\right. ^2}, \end{eqnarray} \tag{ 4b }$$

where the Newtonian versions of kinetic and magnetic energy are appropriate since the conditions are only mildly relativistic. The definitions in Equations (4) satisfy ∫P(k)dk = 〈E〉 for P_K and P_M. Figure 4 shows the power spectrum of kinetic and magnetic energy at various times throughout the growth and saturation of magnetic field. During the kinematic phase, the kinetic energy has a power spectrum P_K(k)∝k^−5/3 consistent with the Kolmogorov theory for incompressible hydrodynamical turbulence, while $P_M(k,t) \propto e^{t/\tau _1} k^{3/2}$ consistent with Kazanstev's model. P_M(k) maintains the same shape, but exponentiates in amplitude at the timescale τ₁ which is controlled by shearing at the resistive scale. According to Kazantsev's model, P_M(k) should peak at the resistive scale. This is consistent with the observed peak in the magnetic energy at roughly 10 grid zones, the same scale at which we observe the viscous cutoff. This is also consistent with Pm = 1 expected from the numerical scheme employed.

Zoom In Zoom Out Reset image size

When the magnetic energy at the resistive scale surpasses the level of the kinetic energy at that scale, P_M(k) changes shape. The equipartition scale ℓ_{K, M} moves into the inertial range and migrates to larger scale until full saturation occurs with ℓ_{K, M} ≈ L/4. The movement of ℓ_{K, M} to larger scale is associated with the formation of coherent and dynamically substantial magnetic structures of increasing size. For Kolmogorov turbulence, ℓ_{K, M}∝t^3/2 (Beresnyak et al. 2009; Schleicher et al. 2013). In the fully saturated state, the magnetic field is in scale-by-scale super-kinetic equipartition throughout the inertial range, with P_M(k) ≈ 4P_K(k). The largest scales remain kinetically dominated so that the numerically converged saturation level is E_M ≈ 0.6 E_K.

In this Letter we have determined the timescale and saturation level of the small-scale turbulent dynamo operating in the conditions of BNS mergers. We have presented numerically converged simulations showing that magnetic fields are amplified to the ∼10¹⁶ G level within a small fraction of the merger dynamical time (t_sat ≲ t_merge), independently of the seed field strength. If hydrodynamical instabilities create fluctuating velocities on the order of 0.1c as indicated by global simulations, then each 1 km³ turbulent volume dissipates ≳ 10⁴⁶ erg of magnetic energy per 0.1 ms. If f_turb represents the fraction of the merger remnant that contains such turbulence, the magnetic energy dissipated during the merger is at least

$$\begin{equation} E_\mathrm{diss} \sim 10^{49}\, \mathrm{erg} \times \left(\frac{f_\mathrm{turb}}{10^{-2}} \right) \left(\frac{t_\mathrm{merge}}{1\, \mathrm{ms}} \right). \end{equation} \tag{ 5 }$$

A fraction of that dissipation will occur through magnetic reconnection in optically thin surface layers, supplying relativistic electrons which synchrotron radiate in the merger remnant magnetosphere. If 5% of that magnetic energy dissipation creates radiation in the 15–150 keV band, then the fluence at 200 Mpc would be 10⁻⁷ erg cm⁻², potentially rendering most merging neutron stars in the advanced LIGO and Virgo detection volumes detectable by Swift BAT. If so, then merging neutron stars are accompanied by a prompt EM counterpart, independently of whether a later merger phase yields a collimated outflow capable of powering a short GRB.

We suggest that merger flares may be present in the current sample of short GRBs and may be roughly isotropic on the sky since they are seen to distances where the cosmological matter distribution becomes homogeneous (e.g., Tanvir et al. 2005). Searches for merger flares should seek to identify short flares, not unlike soft-gamma repeaters, among the short burst population. If mergers also produce short GRBs, then the flaring activity could be detectable as a separate component of the emission, which may appear as a precursor (e.g., Troja et al. 2010) if there is a considerable delay between merger onset and the formation of short GRB engine (e.g., Baiotti et al. 2008). The possibility of prompt flares from BNS mergers has also been suggested to follow from a shock breakout scenario (Kyutoku et al. 2012) or tidal interaction in the late inspiral (Dall'Osso & Rossi 2013).

The presence of strong magnetic fields may also aid in the ejection of neutron-rich material from surface layers of the merger remnant, possibly enhancing the enrichment of interstellar medium by r-process nuclei (Freiburghaus et al. 1999; Rosswog et al. 1999; Hotokezaka et al. 2013). Enhanced production of r-process nuclei also increases the likelihood of EM detection by radioactive decay powered afterglows, or "kilonovae" (Li & Paczyński 1998; Metzger et al. 2010).

Finally, it has been shown that magnetic fields will directly influence the GW signature and remnant disk mass if they exist at the 10¹⁷ G level (Etienne et al. 2012). Such strong fields are unlikely on evolutionary grounds in merging BNSs, but our results suggest they may be revived up to the scale of turbulent cells during the merger. To exert influence, such turbulent cells would have to fill a considerable fraction of the merger volume. As we have shown in this Letter, the overall magnetic energy budget is controlled by the prevalence (f_turb) and vigor (ε_K) of the turbulent volumes. This fact motivates the use of higher resolution global simulations aimed at measuring f_turb and ε_K.

This research was supported in part by the NSF through grant AST-1009863 and by NASA through grant NNX10AF62G issued through the Astrophysics Theory Program. Resources supporting this work were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center.