Turbulent magnetic field amplification in binary neutron star mergers

Magnetic fields are expected to play a key role in the dynamics and the ejection mechanisms that accompany the merger of two neutron stars. General relativistic magnetohydrodynamic (MHD) simulations offer a unique opportunity to unravel the details of the ongoing physical processes. Nevertheless, current numerical studies are severely limited by the fact that any affordable resolution remains insufficient to fully capture the small-scale dynamo, initially triggered by the Kelvin-Helmholtz instability, and later sourced by several MHD processes involving differential rotation. Here, we alleviate this limitation by using explicit large-eddy simulations, a technique where the unresolved dynamics occurring at the subgrid scales (SGS) is modeled by extra terms, which are functions of the resolved fields and their derivatives. The combination of high-order numerical schemes, high resolutions, and the gradient SGS model allow us to capture the small-scale dynamos produced during the binary neutron star mergers, as shown in previous works. Here, we follow the first 50 milliseconds after the merger and, for the first time, we find numerical convergence on the magnetic field amplification, in terms of integrated energy and spectral distribution over spatial scales. Among other results, we find that the average intensity of the magnetic field in the remnant saturates at $\sim {10}^{16}\mathrm{G}$ around 5 ms after the merger. After 20–30 ms, both toroidal and poloidal magnetic field components grow continuously, fed by the winding mechanism that provides a slow inverse cascade, i.e., gradually transferring kinetic into magnetic energy. We find no clear hints for magnetorotational instabilities and no significant impact of the magnetic field on the redistribution of angular momentum in the remnant in our simulations, probably due to the very turbulent and dynamical topology of the magnetic field at all stages, with small-scale components largely dominating over the large-scale ones. Although the magnetic field grows near the rotation axis of the remnant, longer large-eddy simulations are necessary to further investigate the formation of large-scale, helical structures close to the rotational axis, which could be associated to jet formation.

Figure 1

Surfaces of constant density in the remnant. From left to right, we show representative times $t=(4,6,9)\mathrm{ms}$ after the merger. The collision produces a rotating massive remnant that eventually settles down into a spinning neutron star surrounded by an extended torus or outer envelope. The inner surface, at $\rho ={10}^{13}\mathrm{g}/{\mathrm{cm}}^3$, shows the transition between the bulk and envelope. The outer surface, at $\rho ={10}^{12}\mathrm{g}/{\mathrm{cm}}^3$, displays (partially) the shape of the envelope.

Figure 2

Evolution of relevant quantities of the remnant in a meridional plane. The rows show, from top to bottom, the rest-mass density in $\mathrm{g}/{\mathrm{cm}}^3$, magnetic field intensity in Gauss, angular velocity of the fluid $\mathrm{\Omega }$ in $\mathrm{rad}/\mathrm{s}$, and inverse of $\beta$ factor, at $t=(5,10,30,50)\mathrm{ms}$ after the merger, for the medium resolution LES simulation. In the rest-mass density plots (first row), black solid lines represent constant density surfaces at $\rho =(5\times {10}^{10},{10}^{11},{10}^{12},{10}^{13})\mathrm{g}/{\mathrm{cm}}^3$.

Figure 3

Evolution of the magnetic field intensity in the orbital plane. Snapshots at $t=(0.5,1.5,2,2.5,3.5,5,10,15)\mathrm{ms}$ after the merger, displaying the magnetic field strength in Gauss, together with constant density contours at $\rho =({10}^{13},5\times {10}^{14})\mathrm{g}/{\mathrm{cm}}^3$ showing the transition region between the bulk and envelope of the remnant, as well as the location of very high dense regions. The amplification is mostly driven by the KHI, possibly with contributions from the RTI in the initial stages in the envelope (see top left plot). They induce a turbulent small-scale dynamo, generated mostly near the shear layer, that is quickly advected by the fluid motions and diffuses throughout both the bulk and the envelope.

Figure 4

Evolution of the energies. The panels show the rotational, thermal, and magnetic energies, integrated over the full domain. Colors identify the different resolutions, for LES (solid lines) and standard simulations (dashed). Note that these quantities are computed in a postprocessing step from saved 3D data, which is available only every few ms for some portions of the simulations. This leads to an uneven time resolution of the curves, especially for the cases with less output like LR LES.

Figure 5

Evolution of the average magnetic field. Volume average of the magnetic field intensity either in the bulk (top) or in the envelope (bottom) of the remnant. These averages are calculated for all the resolutions (identified by colors), both for the LES (solid lines) and for the standard simulations (dashes). Clearly, the averages in the bulk calculated with LES converge to a well-defined value.

Figure 6

Evolution of the average poloidal and toroidal magnetic field components. Contributions of the poloidal (blue) and toroidal (red) components to the magnetic field intensity, volume averaged in the bulk (top) or in the envelope (bottom), for the medium resolution LES and the high-resolution simulations (line styles).

Figure 7

Evolution of the average magnetic field (and components) in long timescales. Averaged magnetic field intensity and its components, considering the bulk (solid) or the envelope (dashes), for the medium resolution LES. After the amplification growth, both components are comparable, the expected outcome for the saturation state in isotropic turbulence. The poloidal component decreases until $t\sim 30\mathrm{ms}$ while the toroidal one remains constant. At later times, both components grow linearly (i.e., the toroidal component growth 2 times faster than the poloidal one), probably due to the winding mechanism.

Figure 8

Evolution of the (cylindrical) radial distribution of the magnetic field components. The toroidal and poloidal components of the magnetic field, averaged at different cylindrical surfaces $R$ either in the bulk (top) or in the envelope (bottom), for the medium-resolution LES. At later times, both components grow in regions $R\lesssim 6\mathrm{km}$ in the bulk and $R\lesssim 15\mathrm{km}$ in the envelope.

Figure 9

Evolution of the energy spectra. Spectra of the kinetic (solid) and magnetic energy (dashed) at $t=(5,10,20)\mathrm{ms}$, for the medium resolution LES and the high-resolution simulations. The kinetic and magnetic energies almost reach equipartition at high wave numbers (i.e., before the decay induced by numerical dissipation, which is much more prominent in the magnetic spectra). Note that the high-resolution simulations agrees very well with the medium resolution LES. We caution that the power at low wave numbers is based only on few points in discrete 3D Fourier space and corresponds to a small fraction of the total energy. Therefore, it might be subject to stronger fluctuations.

Figure 10

Evolution of the energy spectra at late times. Kinetic and magnetic spectra for the medium resolution cases at $t=(10,20,30,40,50)\mathrm{ms}$. The standard medium resolution simulation is only qualitatively similar to the medium resolution LES. As expected in turbulence with low magnetization, at $t=10\mathrm{ms}$, the kinetic energy spectra show the standard Kolmogorov power law $k^{-5/3}$ (short solid blue line) in the inertial range, while the magnetic energy follows the Kazantsev power law $k^{3/2}$ (short dashed blue line) at low wave numbers. At late times, the magnetic energy spectra follow the Kolmogorov slope at intermediate and high wave numbers. The weighted-average wave numbers $\overline{k}$ are represented by orange and black dots, corresponding to the typical coherent scales of fluid and magnetic structures.

Figure 11

Evolution of the radial distribution of density and angular velocity. Profile of the density and $\mathrm{\Omega }$ in the orbital plane, averaged over the azimuthal direction. Two different times (colors), $t=(10,50)\mathrm{ms}$ are represented for the medium resolution LES (solid) and the standard simulation (dashed). We also include the nonmagnetized case MR B0 (dotted) and show the falloff behavior (green) ${\mathrm{\Omega }}_K\sim R^{-3/2}$ of Keplerian disks for comparison.

Figure 12

Evolution of the inverse of plasma beta. Volume-averaged inverse of plasma beta, $⟨{\beta }^{-1}⟩=⟨B^2/(2P)⟩$, within the bulk and within the envelope, for the medium resolution LES and the high-resolution simulations.

Figure 13

Evolution of the radial distribution of ${\beta }^{-1}$. The profile of the inverse of beta, averaged at different cylindrical surfaces $R$ either in the bulk (top) or in the envelope (bottom), for the medium-resolution LES. At later times, the magnetic pressure becomes more relevant in the envelope for $R\le 15\mathrm{km}$.

Figure 14

Density in a meridional plane. Density at $t=5\mathrm{ms}$ (left panels) and 50 ms (right) after the merger, for the medium resolution simulations: (top) without magnetic field, (middle) standard simulation, and (bottom) LES including the SGS terms. Clearly, stronger magnetic fields increase the matter contamination of the funnel region near the orbital axis.

Figure 15

Distribution of density in a solid angle near the axis. The mass contained in a solid angle of 15° around the $z$ axis, for the medium resolutions at $t=5\mathrm{ms}$ and $t=50\mathrm{ms}$ after the merger. This plot represents how the radial distribution of the mass increases with the magnetic field.

Figure 16

Unbound mass flux and gravitational waveform at different extraction surfaces. Flux of unbound mass (top) and dominant mode (i.e., $l=m=2$) of the Newman-Penrose ${\mathrm{\Psi }}_4$ (bottom), as a function of time. Both quantities are computed at the three spherical surfaces located at $r_{\mathrm{ext}}=(150,300,450)\mathrm{km}$.

Figure 17

Evolution of gravitational wave modes. Top: The module $|C_{l=2,m}|$ of the ${\mathrm{\Psi }}_4$ modes $m=1$ and $m=2$, for the medium resolution simulations. The $m=1$ mode is persistent for all the cases (i.e., LES vs no LES, and magnetized vs unmagnetized). Bottom: Instantaneous frequency of the $l=m=2$ gravitational wave mode for the medium and the high-resolution simulations.