A large-scale magnetic field produced by a solar-like dynamo in binary neutron star mergers

Abstract

The merger of two neutron stars launches a relativistic jet, which must be driven by a strong large-scale magnetic field. However, the magnetohydrodynamical mechanism required to build up this magnetic field remains uncertain. By performing an ab initio super-high-resolution neutrino-radiation magnetohydrodynamics merger simulation in full general relativity, we show that the αΩ dynamo mechanism, driven by the magnetorotational instability, builds up the large-scale magnetic field inside the long-lived remnant of the binary neutron star merger. As a result, the magnetic field induces a Poynting-flux-dominated relativistic outflow with an isotropic equivalent luminosity of ~10⁵² erg s⁻¹ and a magnetically driven post-merger mass ejection of ~0.1 M_⊙. Therefore, the magnetar hypothesis, in which an ultra-strongly magnetized neutron star drives a relativistic jet in binary neutron star mergers, is possible. Magnetars can be the engines of short, hard gamma-ray bursts, and they should be associated with very bright kilonovae, which current telescopes could observe. Therefore, this scenario is testable in future observations.

Main

The observation of GW170817/GRB 170817A/AT 2017gfo made binary neutron star mergers a leading player of multi-messenger astrophysics^1,2,3,4,5. It revealed that at least a part of the origin of short, hard gamma-ray bursts and the heavy elements synthesized by rapid neutron capture (the r process) is binary neutron star mergers^{1,2,3,4,6,7,8,9}.

All the fundamental interactions play an essential role in binary neutron star mergers. Thus, the method chosen to explore them theoretically is a numerical relativity simulation that implements all the effects of the fundamental interactions. Numerical relativity simulations can qualitatively explain AT 2017gfo, that is, the kilonova emissions associated with the radioactive decay of the r-process elements^{10,11,12,13,14,15,16,17,18}. However, a quantitative understanding of this event is still being developed. Moreover, there is no theoretical consensus about how the binary neutron star merger drove the short, hard gamma-ray burst^19,20,21.

A relativistic jet was launched from this binary neutron star merger, which was observed as a short and hard gamma-ray burst^1,2,3,4. Such relativistic jets are most likely driven by a magnetohydrodynamics process. This indicates that a binary neutron star merger remnant must build up a large-scale magnetic field with the dynamo to launch the jet^22,23. However, the mechanism behind the large-scale dynamo in the merger remnant still needs to be clarified^24,25.

Also, the site for generating the large-scale magnetic field after the merger is a riddle. Is it inside a massive torus around a black hole that formed after a merger remnant collapsed into it? Or is it the long-lived remnant of a massive neutron star²⁶? Recent numerical simulations for a black hole and torus system, as a remnant of a binary neutron star merger or a black hole–neutron star merger, suggest that a large-scale magnetic field is established after a long-term evolution of O(1–10) s (refs. ^20,27,28,29). Consequently, a relativistic jet is launched due to the Blandford–Znajek mechanism³⁰. However, the long-lived remnant of a massive neutron star is more computationally challenging than the black hole plus torus system because the requirement to resolve relevant magnetohydrodynamical instabilities numerically is severe³¹. The physical mechanism for generating the large-scale magnetic field and launching the jet could differ in the two scenarios mentioned above. It could be observationally testable if we managed to build a long-lived binary neutron start model that generates a large-scale B field and launches a jet.

We tackle this problem with a super-high-resolution neutrino-radiation magnetohydrodynamics simulation of a binary neutron star merger in full general relativity.

We employ the latest version of our numerical relativity neutrino-radiation magnetohydrodynamics code³². We employ a static mesh refinement with 2:1 refinements to cover a wide dynamic range. For the simulations in this article, the grid structure consists of 16 Cartesian domains. Each domain has 2N × 2N × N elements in the x, y and z directions, and we assume orbital plane symmetry. We set N = 361 and Δx_finest = 12.5 m. The grid resolution is the highest among the binary neutron star merger simulations³¹. The initial orbital separation is ~44 km (Methods).

We employ DD2 (ref. ³³) as an equation of state for the neutron star matter and symmetric binary, which has a total mass of 2.7 M_⊙. With this set-up, a hypermassive neutron star that transiently formed after the merger will survive for >O(1) s (ref. ³⁴).

The purely poloidal magnetic field loop is embedded inside the neutron stars and has a maximum field strength of ${B}_{0,\max }=1{0}^{15.5}$ G (ref. ³⁵). It is much stronger than the 10⁷–10¹¹ G observed for binary pulsars (ref. ³⁶). The available computational resources limit us to a strong field strength and idealized topology. Nonetheless, it is natural to anticipate that in reality the amplification of the magnetic field leads to a high field strength (>10¹⁶ G) in a short timescale after the merger (see below).

Figure 1 plots the evolution of the electromagnetic energy as a function of the post-merger time. As reported in refs. ^31,35, the electromagnetic energy was exponentially amplified shortly after the merger due to the Kelvin–Helmholtz instability, which emerges at the contact interface when the binary merges. A part of the turbulent kinetic energy due to the instability is converted to electromagnetic energy³⁵. In the inset, we generate the same plot for −1 ≤ t − t_merger ≤ 5 ms but different initial magnetic field strengths of ${B}_{0,\max }=1{0}^{15}$ G (red) and ${B}_{0,\max }=1{0}^{14}$ G (orange) while keeping the grid resolution. Also, we plot the results for the other grid resolutions of Δx_finest = 100 or 200 m (brown and purple curves, respectively), while keeping ${B}_{0,\max }=1{0}^{15.5}$ G. The figure clearly shows that the growth rate of the electromagnetic energy for 0 ≲ t − t_merger ≲ 1 ms does not depend on the initial field strength but on the grid resolution, as expected from the properties of the Kelvin–Helmholtz instability^37,38 (Extended Data Fig. 1). The post-merger amplification of the magnetic field due to this instability terminates when the shear layer is dissipated by the shock waves at t − t_merger ≈ 2 ms.

Fig. 1: Overview of the magnetic field evolution in the remnant of the binary neutron star.

Top, Electromagnetic energy as a function of the post-merger time for the total (solid), the poloidal (dashed) and the toroidal (dotted) components. The purple curves are for the simulation with Δx_finest = 200 m. The inset shows how the amplification of the magnetic field due to the Kelvin–Helmholtz instability depends on the initial magnetic field strength and grid resolution. B14 and B15 mean B_0,max= 10¹⁴ G and 10¹⁵ G, respectively. Bottom, Magnetic field lines for a density of ρ < 10¹³ g cm⁻³ at t − t_merger ≈ 130 ms. The core of the hypermassive neutron star is shown for a density of ρ > 10¹³ g cm⁻³.

At t − t_merger ≈ 5 ms, the electromagnetic energy temporarily settles to ~3 × 10⁵⁰ erg. However, the toroidal magnetic field (dotted curve) is subsequently amplified. In particular, its contribution to the total electromagnetic energy becomes prominent for t − t_merger ≳ 20 ms and the growth rate is proportional to t², which indicates that the magnetic winding is efficient with a coherent poloidal magnetic field. The differential rotational energy of ~1–2 × 10⁵³ erg of the remnant of the massive neutron star is the energy budget for the magnetic winding. Once the magnetic braking starts to work, the electromagnetic energy saturates around ~5 × 10⁵¹ erg. We also plot the result of the simulation with Δx_finest = 200 m in Fig. 1. With this resolution, it is hard to resolve the Kelvin–Helmholtz and magnetorotational instability, particularly in the high-density region (Extended Data Figs. 1–3, and the section on convergence in Methods). Consequently, the amplification of the electromagnetic energy is less efficient in the low-resolution case. In particular, there is a striking difference in the poloidal magnetic field between the simulations with Δx_finest = 12.5 and 200 m. Because the magnetic field that is amplified due to the Kelvin–Helmholtz instability is randomly oriented³¹, there must be a mechanism to generate the coherent poloidal magnetic field.

To make it clear which part of the merger remnant is responsible for the generation of the coherent magnetic field, we evaluate the electromagnetic energy in the region where the magnetorotational instability is active, as defined by ρ ≤ 10^14.5 g cm⁻³ and as shown in Fig. 2 (Extended Data Fig. 2). We also decompose the contributions from the mean poloidal and toroidal magnetic fields by taking an axisymmetric average. The electromagnetic energy in this region is subdominant compared to that contained deep inside the core region where ρ ≥ 10^14.5 g cm⁻³ (see also Fig. 1). However, in the simulation with Δx_finest = 12.5 m, the mean poloidal field exponentially grows in the range 20 ≲ t − t_merger ≲ 50 ms. At t − t_merger ≈ 20 ms, its contribution to the total poloidal field energy (solid cyan curve) is ~1% and becomes ~10% after the exponential growth. Such a clear exponential growth is invisible in the simulation with Δx_finest = 200 m for t − t_merger ≲ 100 ms. Also, the mean poloidal field energy differs by one or two orders of magnitude in the simulations with Δx_finest = 12.5 or 200 m. For the toroidal component, the mean and total field energy are of the same order of magnitude for t − t_merger ≳ 30 ms. We also confirm that the magnetorotational instability is responsible for generating the mean poloidal magnetic flux (Extended Data Figs. 4 and 5). These results indicate that the generation of a coherent magnetic field is triggered by the magnetorotational instability.

Fig. 2: Generation of mean electromagnetic field.

Mean and total electromagnetic energy in the region where the magnetorotational instability is active with ρ ≤ 10^14.5 g cm⁻³ as a function of the post-merger time. The blue, green, cyan and purple curves denote the mean poloidal, mean toroidal, total poloidal and total toroidal components, respectively, in the high-resolution simulation. The dashed curves are for the simulation with Δx_finest = 200 m.

The αΩ dynamo, which is driven by the magnetorotational instability³⁹ in the current context, is a potential mechanism for generating the large-scale magnetic field^22,23. In the mean field dynamo theory, we assume that each physical quantity Q is composed of the mean field $\bar{Q}$ and the fluctuation q, that is, $Q=\bar{Q}+q$. Thus, we cast the induction equation as

$$\begin{array}{r}{\partial }_{t}\bar{{{{\bf{B}}}}}=\nabla \times \left(\bar{{{{\bf{U}}}}}\times \bar{{{{\bf{B}}}}}+\bar{{{{\bf{{{{\mathcal{E}}}}}}}}}\right),\end{array}$$

(1)

where $\bar{{{{\bf{{{{\mathcal{E}}}}}}}}}=\overline{{{{\bf{u}}}}\times {{{\bf{b}}}}}$ is the electromotive force due to the fluctuations, ${{{\bf{B}}}}=\bar{{{{\bf{B}}}}}+{{{\bf{b}}}}$ is the magnetic field and ${{{\bf{U}}}}=\bar{{{{\bf{U}}}}}+{{{\bf{u}}}}$ is the velocity field. Note that the magnetorotational instability-driven turbulence produces the fluctuations u and b. In the αΩ dynamo, we express the electromotive force as a function of the mean magnetic field:

$$\begin{array}{r}{\bar{{{{\mathcal{E}}}}}}_{i}={\alpha }_{ij}{\bar{B}}_{j}+{\beta }_{ij}{\left(\overline{\nabla \times \bar{B}}\right)}_{j},\end{array}$$

(2)

where α_ij and β_ij are tensors that do not depend on ${\bar{B}}_{j}$. We calculate the mean field by taking the average in the azimuthal direction.

In the presence of a cylindrical differential rotation, the simplest mean field dynamo is an αΩ dynamo, in which the toroidal magnetic field is generated by the shear of the poloidal magnetic field by the differential rotation ${\bar{U}}_{\phi }$, which is also called the Ω effect. The decrease in the rotation with radius and equation (1) imply that ${\bar{B}}_{\phi }$ should be anti-correlated with ${\bar{B}}_{R}$. To complete the dynamo cycle, the poloidal magnetic field has to be generated by the toroidal electromotive force ${\bar{{{{\mathcal{E}}}}}}_{\phi }$, whose main contribution is from a diagonal component α_ϕϕ, which is the so-called α effect. ${\bar{{{{\mathcal{E}}}}}}_{\phi }$ is, therefore, correlated or anti-correlated to ${\bar{B}}_{\phi }$ depending on the sign of α_ϕϕ (equation (2)). This complete cycle forms a dynamo wave that oscillates with a period of ${P}_{\rm{theory}}=2\pi\left\vert\frac{1}{2}{\alpha }_{\phi \phi }\frac{\mathrm{d}\varOmega}{\mathrm{d}\ln R}{k}_{z}\right\vert^{-1/2}$ (refs. ^22,23) and propagates in the direction α_ϕϕ∇Ω × e_ϕ according to the Parker–Yoshimura rule^40,41. Here k_z is the wavenumber of the dynamo wave in the vertical direction. In this theoretical description, we have supposed that contributions from the other α_ij components and the turbulent resistivity tensor β_ij are subdominant. We will show that it is the case in the following.

First, we confirm that the simulation set-up employed can capture the fastest-growing mode of the magnetorotational instability in the outer region of the hypermassive neutron star (Extended Data Fig. 2). Therefore, the turbulent state is developed. Figure 3 plots the butterfly diagrams for ${\bar{B}}_{\phi },{\bar{B}}_{R}$ and ${\bar{{{{\mathcal{E}}}}}}_{\phi }$ at R = 30 km. The top left and top right panels show the anti-correlation between ${\bar{B}}_{\phi }$ and ${\bar{B}}_{R}$, which indicates the Ω effect. To quantify the correlation between ${\bar{B}}_{\phi }$ and ${\bar{{{{\mathcal{E}}}}}}_{\phi }$, we compute the Pearson correlation C_P(X, Y) between the two quantities X and Y in the bottom right panel of Fig. 3. This figure shows that ${\bar{B}}_{\phi }$ and ${\bar{{{{\mathcal{E}}}}}}_{\phi }$ are anti-correlated for z ≲ 15 km, where the pressure scale height at R = 30 km is ~14.6 km. The correlation between the electromotive force and mean current ${\bar{J}}_{i}={\left(\overline{\nabla \times \bar{B}}\right)}_{i}$ is small, and thus, the β_ij tensor can be neglected (Extended Data Fig. 6). Therefore, the generation of the mean poloidal magnetic field is determined primarily by the α_ij tensor. Although α_ϕR has a non-negligible contribution, the contribution of α_ϕϕB_ϕ dominates in the turbulent region for z ≲ 10 km (Extended Data Fig. 6). This shows that α_ϕϕ is the main component, and it is plotted in the bottom right panel of Fig. 3. We also confirmed that the α²Ω dynamo is irrelevant compared to the αΩ dynamo (Extended Data Fig. 7).

Fig. 3: αΩ dynamo inside the remnant of the binary neutron star merger.

Butterfly diagrams at R = 30 km. Top left, Mean toroidal magnetic field ${\bar{B}}_{\phi }$. Top right, Mean radial magnetic field ${\bar{B}}_{R}$. Bottom left, Toroidal electromotive force ${\bar{{{{\mathcal{E}}}}}}_{\phi }$. Bottom right, Parameters α_ϕϕ (orange) and correlation between ${\bar{{{{\mathcal{E}}}}}}_{\phi }$ and ${\bar{B}}_{\phi }$ (blue).

With these quantities, we can predict the period of the αΩ dynamo, as listed in Table 1. We measured the shear rate $q=-\mathrm{d}\ln\varOmega/\mathrm{d}\ln R$ at the selected radius and choose the wavenumber k_z corresponding to the pressure scale height, which is the most extended vertical length in the turbulent region. The sixth and seventh columns report the predicted period of the αΩ dynamo and the period measured in the butterfly diagram. The agreement at R = 30 km is remarkable. We show the comparison at different radii from R = 20 to 50 km in the table, and the agreement is also reasonable. In addition, since α_ϕϕ is negative and q is positive, the dynamo wave propagates along the direction of the Parker–Yoshimura rule, that is, the z direction (http://www2.yukawa.kyoto-u.ac.jp/~kenta.kiuchi/anime/SAKURA/Br_core_cut_9km.m4v). With these findings, we conclude that the dynamo in our simulation can be interpreted as an αΩ dynamo, and it builds up the large-scale magnetic field in the remnant of the hypermassive neutron star (Fig. 1).

Table 1 The αΩ dynamo period prediction and simulation data at several radii

After the development of the coherent poloidal magnetic field due to the αΩ dynamo and the resultant efficient magnetic winding, a magnetic-tower outflow is launched towards the polar direction (http://www2.yukawa.kyoto-u.ac.jp/~kenta.kiuchi/anime/SAKURA/DD2_135_135_Dynamo.mp4; Methods). The mean poloidal magnetic flux is generated in the region where the magnetorotational instability is active. The mean magnetic field deep inside the core (ρ ≥ 10^14.5 g cm⁻³), which could be a relic of the initial field, stays buried below r ≲ 10 km in the polar region throughout the simulation (http://www2.yukawa.kyoto-u.ac.jp/~kenta.kiuchi/anime/SAKURA/movie_Mean_Poloidal_Flux.mp4; Methods). This indicates that the large-scale field generated by the αΩ dynamo is responsible for launching the jet. We confirm that the low-resolution simulation cannot capture the launching of the strong Poynting-flux-dominated outflow or the enormous post-merger ejecta²¹ (see also Extended Data Fig. 8).

The luminosity of the Poynting flux, which is defined by ${L}_{{{{\rm{Poy}}}}}=$ $-{\oint }_{r\approx 500{\rm{km}}}\sqrt{-g}\,\left({{T}^{r}}_{t}\right)_{\rm{(EM)}}\,\mathrm{d}\varOmega$, where ${T}_{{{{\rm{(EM)}}}}}^{\mu \nu }$ is the stress energy tensor for the electromagnetic field, is ~10⁵¹ erg s⁻¹ at the end of the simulation when t − t_merger ≈ 150 ms. The high Poynting flux is confined in the region with θ ≲ 12^∘ where θ is a polar angle, as plotted in the top left panel of Fig. 4 (see also Extended Data Fig. 8). The luminosity and jet opening angle are broadly compatible with some of the observed short, hard gamma-ray bursts⁴². We estimate the terminal Lorentz factor Γ_∞ by the magnetization parameter σ_LC = b²/4πρc² at the light cylinder radius ${r}_{{{{\rm{LC}}}}}=c/\varOmega \approx 40\,{\rm{km}}$ ($\left(\varOmega /7,000\,{\rm{rad}}\,{\rm{s}}^{-1}\right)^{-1}$) contained in a region θ < 12° (refs. ^43,44). We found the baryon loading is still severe, for example, ${\dot{M}}_{\rm{outflow}} \approx10^{-2}\,{M}_{\odot }\,{\rm{s}}^{-1}$ in the polar region. Nonetheless, Γ_∞ increases with time but fluctuates. It reaches ~10–20 at the end of the simulation. Therefore, if the magnetic reconnection efficiently dissipated the Poynting flux, a Lorentz factor of ~10 could be possible⁴⁴. Also, the evacuation due to the strong Poynting flux in the polar region continues at the end of the simulation. Therefore, a higher Γ_∞ could be possible after the long-term evolution.

Fig. 4: Electromagnetic signals from the remnants of the binary neutron star merger.

Top left, Angular distribution of the Poynting flux on a sphere with r ≈ 500 km at t − t_merger ≈ 150 ms. Top right, Electron fraction distribution for the magnetically driven post-merger ejecta (dashed) and the dynamical ejecta (dotted). Bottom, Same as the top right panel but for the terminal velocity (left) and entropy per baryon (right).

The dynamical ejecta driven by the tidal force and the shock wave during the merger phase represent ~0.002 M_⊙ (ref. ³⁴). After the development of the coherent magnetic field, we observe a new component in the ejecta that is driven by the Lorentz force, namely, the magnetically driven wind⁴⁵. The mass of this component is ~0.1 M_⊙ at the end of the simulation t − t_merger ≈ 150 ms. The electron fraction of the post-merger ejecta shows a peak around ~0.2 with an extension to ~0.5. The terminal velocity of the post-merger ejecta peaks around ~0.08c, where c is the speed of light. Therefore, it could contribute to a bright kilonova emission through the synthesis of r-process heavy elements^26,46.

We speculate that the high-luminosity state and the post-merger mass ejection would continue for O(1) s, which corresponds to the neutrino cooling timescale^34,47. As reported in refs. ^27,34, the torus expands due to the angular momentum transport facilitated by the magnetorotational instability-driven turbulent viscosity after the neutrino cooling becomes inefficient. As a result, the funnel region above the remnant neutron star expands and the jet could no longer be collimated. A long-term simulation for O(1) s of a remnant of a massive neutron star is future work to be pursued.

Also, a simulation with initially weakly magnetized binary neutron stars is necessary to confirm the picture reported in this article because the saturated magneto-turbulent state and the generation timescale of the coherent poloidal magnetic field due to the magnetorotational instability could depend on the initial magnetic field strength and flux⁴⁸. The recent large-eddy simulations of a magnetized binary neutron star merger may give a hint about this issue^49,50. Even if the initial magnetic field has a much weaker strength or different topology from the one we assumed, the saturation strength and field profile due to the Kelvin–Helmholtz instability in the merger remnant are like those that we found in this article. Also, how the mean poloidal magnetic field is set in reality after the merger is an open problem. If the mean poloidal magnetic field just after the merger is a relic of the pre-merger poloidal field, that is, 10¹⁰–10¹¹ G at maximum, it may take $\ln (1{0}^{4})\times 0.02\,{{{\rm{s}}}} \approx0.2\,{{{\rm{s}}}}$ to reach the saturation strength of 10¹⁴–10¹⁵ G, as we assume the mean poloidal field is exponentially amplified with the period of the αΩ dynamo²² (Fig. 2 and Table 1). Consequently, the jet may be launched at O(0.1) s after the merger in reality. However, the interior structure of the magnetic field in the pre-merger neutron stars is not well understood, and the magnetic reconnection of the fluctuating poloidal field generated by the Kelvin–Helmholtz instability may enhance the mean poloidal field after the merger. This situation should be explored as a future task.

To summarize, we generated a large-scale magnetic field in a long-lived remnant of a binary neutron star merger using a super-high-resolution neutrino-radiation magnetohydrodynamics simulation in general relativity. The magnetorotational instability-driven αΩ dynamo generates the large-scale magnetic field. As a result, the launch of the Poynting-flux-dominated relativistic outflow in the polar direction and an enormous amount of magnetically driven wind are induced. Our simulation suggests that the magnetar engine generates short, hard gamma-ray bursts and bright kilonovae emission, which could be observed in the near-future multi-messenger observations.

Methods

Numerical method

Our code implements the Baumgarte–Shapiro–Shibata–Nakamura-puncture formulation to solve Einstein’s equation^51,52,53,54. The code also uses the Z4c prescription to suppress the constraint violation⁵⁵. The fourth-order accurate finite difference is used as a discretization scheme in space and time. The sixth-order Kreiss–Oliger dissipation is also employed. The HLLD Riemann solver⁵⁶ and the constrained transport scheme⁵⁷ are used to solve the equations of motion of the relativistic magnetohydrodynamic fluid. The neutrino-radiation transfer is solved by the grey M1+GR-leakage scheme⁵⁸ to take into account neutrino cooling and heating.

Initial data

We adopted quasi-equilibrium initial data of irrotational binary neutron stars in the neutrino-free beta equilibrium derived in a previous paper³⁴ using the public spectral library LORENE (ref. ⁵⁹). The initial orbital angular velocity is Gm₀Ω₀/c³ ≈ 0.028 with m₀ = 2.7 M_⊙. The residual orbital eccentricity is ~10⁻³. For our high-resolution study, the data are remapped onto the computational domain described in the ‘Grid set-up’ section.

Grid set-up

We employ the static mesh refinement with 2:1 refinement, that is, a coarser domain with a grid resolution twice that of a finer domain. All the domains are composed of concentric Cartesian domains with a fixed grid size N. The size of the grid is 2N × 2N × N in the x, y and z directions, and we assume orbital plane symmetry. We employ 16 domains with N = 361 and the finest grid resolution of Δx_finest = 12.5 m. The first three finest domains, whose sizes are [−4.5: 4.5 km]² × [0: 4.5 km], [−9: 9 km]² × [0: 9 km] and [−18: 18 km]² × [0: 18 km], respectively, are employed to resolve the Kelvin–Helmholtz instability, which emerges on a contact interface (shear layer) when the two neutron stars merge. The initial binary separation is ~44 km, and the coordinate radius of the neutron star is 10.9 km. Thus, the fourth finest domain with [−36: 36 km]² × [0: 36 km] covers the entire binary neutron star.

We started each simulation with Δx_finest = 12.5 m and 16 static mesh refinement domains. At ~30 ms after the merger, we removed the two finest domains with Δx_finest = 12.5 m and Δx = 25 m and continued the simulation with Δx_finest = 50 m until ~50 ms after the merger. Then, we removed the domain with Δx_finest = 50 m and continued the simulation with Δx_finest = 100 m. The timing for removing the finer domains was determined by monitoring the magnetorotational instability quality factor so as not to go below the critical value of 10 after the removal (see below). This strategy allowed us to capture the efficient amplification of the magnetic field through the Kelvin–Helmholtz instability and resolve the magnetorotational instability inside the remnant while saving computational costs.

To check the validity of removing finer domains, we continued the simulation with Δx_finest = 12.5 m up to t − t_merger ≈ 33 ms (see the blue dashed curve in Fig. 1). We observed an ~10% decrease in the poloidal electromagnetic energy in the simulation with Δx_finest = 50 m compared to that with Δx_finest = 12.5 m around this time. Nonetheless, the poloidal electromagnetic energy increased around t − t_merger ≈ 40 ms in the simulation with Δx_finest = 50 m (the green dashed curve in Fig. 1) due to the α effect we discuss in the main text. Also, there was no substantial decrease in the toroidal electromagnetic energy after removing the finer domains at t − t_merger ≈ 30 ms. Similarly, there was no visible deterioration after the second removal of the finer domain at t − t_merger ≈ 50 ms.

For the convergence test, we also performed the simulations with Δx_finest = 100 (200) m and N = 361 (185). The number of domains was 13.

Equation of state

We extended the original DD2 equation of state to the low-density and temperature region with the Helmholtz equation of state⁶⁰. Because any high-resolution shock-capturing scheme does not allow a vacuum state, we employed the atmospheric prescription outside the neutron stars. Specifically, we set a constant atmospheric density of 10³ g cm⁻³ inside r ≤ L_atm and assumed a power-law profile ${\rho}_{\rm{atm}}=\max [10^{3}{({L}_{{{{\rm{atm}}}}}/r)}^{3}\,{\rm{g}}\,{\rm{cm}}^{-3},{\rho }_{{{{\rm{fl}}}}}]$ for r > L_atm where L_atm = 36 km. The Helmholtz equation of state employed determines the floor value of the rest-mass density, which is ρ_fl = 0.167 g cm⁻³. We also assumed a constant atmospheric temperature of 10⁻³ MeV.

Kelvin–Helmholtz instability

The top panel of Extended Data Fig. 1 is the same as the inset in Fig. 1 but with additional data. The solid blue, solid brown and solid purple curves are the results with Δx_finest = 12.5, 100 and 200 m, respectively, with ${B}_{0,\max }=1{0}^{15.5}$ G. The solid red and solid orange curves show the results with Δx_finest = 12.5 m with ${B}_{0,\max }=1{0}^{15}$ and 10¹⁴ G, respectively. The cyan curve shows the result with Δx_finest = 18.75 m and ${B}_{0,\max }=1{0}^{14}$ G. The dotted blue and dotted red curves plot the result with Δx_finest = 12.5 m and ${B}_{0,\max }=1{0}^{14}$ G magnified by the square of the ratio of the initial magnetic field strength, that is, ${(1{0}^{15.5}/1{0}^{14})}^{2}=1{0}^{3}$ and ${(1{0}^{15}/1{0}^{14})}^{2}=1{0}^{2}$, respectively.

This figure suggests the following two points. First, the solid and dotted blue curves and the solid and dotted red curves overlap until the back reaction starts to activate, which happens when the electromagnetic energy reaches ~3 × 10⁴⁹ erg. This implies that the magnetic field is passive for t − t_merger ≲ 1 ms, which is the linear phase. The saturation of the electromagnetic energy due to the Kelvin–Helmholtz instability is likely to be ~1–3 × 10⁵⁰ erg corresponding to O(0.1)% of the internal energy^31,50,61,62.

The second point is that the growth rate of the electromagnetic energy in the linear phase depends on the grid resolution not on the initial magnetic field strength. To quantify the dependence of the growth rate on the grid resolution and initial magnetic field strength, we estimated the growth rate by fitting the electromagnetic energy as ${E}_{{{{\rm{mag}}}}}(t)=A\exp \left({\sigma }_{{{{\rm{KH}}}}}(t-{t}_{{{{\rm{merger}}}}})\right)$ for 0 ≲ t − t_merger ≲ 1 ms, which corresponds to the linear phase. The bottom panel of Extended Data Fig. 1 plots the estimated growth rate as a function of the initial magnetic field strength. The symbols denote the grid resolution. With Δx_finest = 12.5 m, the growth rate was ~7 ms⁻¹, irrespective of the initial magnetic field strength. The growth rate increased with the grid resolution. With Δx_finest = 100 or 200 m, the efficient amplification cannot be captured. All the properties are consistent with those of the Kelvin–Helmholtz instability, that is, the smaller the scale of the vortices is, the larger the growth rate is⁶³. The reality is located in the leftward and upward direction in this diagram, that is, the weak magnetic field observed in binary pulsars (the leftward direction in the diagram) and infinitesimal grid resolution (the upward direction in the diagram). Therefore, we anticipate that the electromagnetic energy will saturate in a very short timescale, ≪1 ms after the binary neutron star merger, irrespective of the magnetic field strength at the onset of the merger unless a black hole is promptly formed^31,35,64.

We estimated how small Δx_finest should be to simulate the amplification due to the Kelvin–Helmholtz instability starting from the ‘realistic’ initial magnetic field strength of 10¹¹ G (ref. ³⁶). The saturation strength found in the simulation with Δx_finest = 12.5 m is likely to be $| \bar{B}| \approx 1{0}^{16.5}$ G. We fitted the growth rate as a function of the grid resolution from Extended Data Fig. 1. It was σ_KH(ms⁻¹) = 90/Δx_finest(m), where we assumed the inverse proportionality would originate due to the Kelvin–Helmholtz instability^37,38. Assuming that the initial magnetic field strength was 10¹¹ G and that the lifetime of the shear layer was 2 ms, we estimated the required growth rate as

$${\sigma }_{{{{\rm{KH}}}}}=\frac{1}{2\,{\rm{ms}}}\ln{\left(\frac{1{0}^{16.5}\,{\rm{G}}}{1{0}^{11}\,{\rm{G}}}\right)}^{2}=12.7\,{\rm{ms}}^{-1}.$$

(3)

Therefore, the required grid resolution is

$$\begin{array}{r}\Delta {x}_{{{{\rm{finest}}}}}=\frac{90}{{\sigma }_{{{{\rm{KH}}}}}}\approx 7.1\,{{{\rm{m}}}}.\end{array}$$

(4)

Magnetorotational instability, neutrino viscosity and neutrino drag

To quantify how well the magnetorotational instability is resolved in our simulation, we estimated the rest-mass-density-conditioned magnetorotational instability quality factor, as defined by

$${\langle {Q}_{\rm{MRI}}\rangle }_{\rho }\equiv \frac{{\langle {\lambda }_{\rm{MRI}}\rangle}_{\rho}}{\Delta x}=\frac{1}{\Delta x}\frac{{\int}_{\rho}{\lambda }_{\rm{MRI}}\,\mathrm{d}^{3}x}{{\int}_{\rho}\,\mathrm{d}^{3}x},$$

(5)

where ${\lambda }_{{{{\rm{MRI}}}}}=2\uppi {B}^{z}/\left(\sqrt{4\pi \rho }\varOmega \right)$ is the fastest-growing mode of an ideal magnetorotational instability³⁹. As the condition in terms of the rest-mass density, we define a remnant core and a remnant torus as fluid elements with ρ ≥ 10¹³ and <10¹³ g cm⁻³, respectively. Furthermore, we exclude the core region above 10^14.5 g cm⁻³ in the estimate of the quality factor because in such a region, the radial gradient of the angular velocity is positive, and it is not subject to magnetorotational instability³⁹, as plotted in the top panel of Extended Data Fig. 2. Also, we introduce a cutoff density of 10⁷ g cm⁻³ for the torus to suppress the overestimation of the quality factor in the low-density polar region. The middle and bottom panels of the figure show that the fastest-growing mode of the magnetorotational instability is well resolved, both in the core and torus, throughout the simulation. Consequently, magneto-turbulence develops and is sustained. However, the magnetorotational instability, particularly in the core throughout the entire simulation and in the torus for t − t_merger ≲ 80–90 ms, cannot be resolved when Δx_finest = 200 m. Thus, the turbulence is not fully produced in such a low-resolution run.

One caveat is that the neutrino viscosity and drag could affect the magnetorotational instability as diffusive and damping processes⁶⁵. The former (latter) becomes relevant when the neutrino mean free path is shorter (longer) than the wavelength of the magnetorotational instability. Given a profile of the merger remnant, such as the density ρ, angular velocity Ω, temperature T and hypothetical magnetic field strength ${B}_{{{{\rm{hyp}}}}}^{z}$, we solve the two branches of the dispersion relations to quantify the neutrino viscosity and drag effects⁶⁵: For the neutrino viscosity:

$$\begin{array}{r}{\left[\left(\tilde{\sigma }+{\tilde{k}}^{2}{\tilde{\nu }}_{\nu }\right)\tilde{\sigma }+{\tilde{k}}^{2}\right]}^{2}+{\tilde{\kappa }}^{2}\left[{\tilde{\sigma }}^{2}+{\tilde{k}}^{2}\right]-4{\tilde{k}}^{2}=0,\end{array}$$

(6)

and for the neutrino drag:

$$\begin{array}{r}{\left[\left(\tilde{\sigma }+{\tilde{\Gamma }}_{\nu }\right)\tilde{\sigma }+{\tilde{k}}^{2}\right]}^{2}+{\tilde{\kappa }}^{2}\left[{\tilde{\sigma }}^{2}+{\tilde{k}}^{2}\right]-4{\tilde{k}}^{2}=0,\end{array}$$

(7)

where $\tilde{\sigma }=\sigma /\varOmega ,\tilde{k}=k{v}_\mathrm{A}/\varOmega,{\tilde{\kappa }}^{2}={\kappa }^{2}/{\varOmega }^{2},{\tilde{\nu }}_{\nu }={\nu }_{\nu }\varOmega /{v}_\mathrm{A}^{2}$ and ${\tilde{\varGamma }}_{\nu }={\varGamma }_{\nu }/\varOmega$. σ and k are the growth rate and wavenumber of the unstable mode of the magnetorotational instability. ${v}_\mathrm{A}={B}_{\rm{hyp}}^{z}/\sqrt{4\uppi \rho }$ is the Alfvén wave speed. κ² is the epicyclic frequency. ν_ν and Γ_ν are the neutrino viscosity and drag damping rate, respectively. Reference ⁶⁵ provides fitting formulae for them as functions of the rest-mass density and temperature:

$${\nu }_{\nu }=1.2\times 1{0}^{10}{\left(\frac{\rho}{10^{13}\,{\rm{g}}\,{\rm{cm}}^{-3}}\right)}^{-2}{\left(\frac{T}{10\,{\rm{MeV}}}\right)}^{2}\,{\rm{cm}}^{2}\,\mathrm{s}^{-1},$$

(8)

$${\varGamma}_{\nu}=6\times10^{3}{\left(\frac{T}{10\,{\rm{MeV}}}\right)}^{6}\,{\rm{s}}^{-1}.$$

(9)

Also, the neutrino mean free path l_ν is fitted by

$${l}_{\nu }=2\times 1{0}^{3}{\left(\frac{\rho}{10^{13}\,{\rm{g}}\,{\rm{cm}}^{-3}}\right)}^{-1}{\left(\frac{T}{10\,{\rm{MeV}}}\right)}^{-2}\,{\rm{cm}}.$$

(10)

Extended Data Fig. 3 plots the growth rate of the magnetorotational instability for a given profile of a remnant of a massive neutron star and a hypothetical value of ${B}_{{{{\rm{hyp}}}}}^{z}$. We use the profiles on the orbital plane at t − t_merger ≈ 10 ms. The purple curve denotes the boundary where the condition ${\tilde{\nu }}_{\nu }$ or ${\tilde{\Gamma }}_{\nu }\approx 1$ is met⁶⁵. Inside the boundary, the neutrino viscosity or the neutrino drag substantially suppresses the growth rate. Outside it, the growth rate is essentially the same as for the ideal magnetorotational instability. Above it, we plot the azimuthally averaged magnetic field strength ${\bar{B}}_{{{{\rm{sim}}}}}^{z}$ obtained in the simulation. Because of the efficient amplification due to the Kelvin–Helmholtz instability just after the merger, the neutrino viscosity and drag effects are irrelevant in the entire region of the merger remnant.

Generation of mean poloidal magnetic flux due to the magnetorotational instability

To validate our interpretation that the αΩ dynamo is responsible for generating the large-scale magnetic field, we evaluate the mean poloidal magnetic flux on a certain sphere of radius r:

$${\varPhi}_{{\bar{B}}_{r}}(r)=2\int\nolimits_{0}^{\pi /2}{\bar{B}}_{r}{r}^{2}\sin \theta\,\mathrm{d}\theta ,$$

(11)

where ${\bar{B}}_{r}$ denotes the radial mean field in the spherical-polar coordinate.

Extended Data Fig. 4 plots the evolution of the mean poloidal magnetic flux. We selected 8 and 20 km as representative radii for the magnetorotational instability inactive and active regions (see the top panel of Extended Data Fig. 2). On the one hand, in the simulation with Δx_finest = 12.5 m, the mean poloidal magnetic flux on the sphere with r = 8 km exhibits intensive time variability for t − t_merger ≲ 20 ms, which reflects the oscillations of the remnant core, that is, compression and decompression. The poloidal flux gradually decreases, presumably due to the magnetic reconnection, which is imprinted as the decrease of the total poloidal field energy for 10 ≲ t − t_merger ≲ 20 ms in Fig. 1. Nonetheless, it relaxes to a roughly constant value and stays there until t − t_merger ≲ 120 ms. This reflects that magnetorotational instability is inactive in this region, that is, the mean poloidal magnetic flux is not generated. On the other hand, it is evident that the mean poloidal flux on the sphere with r = 20 km is generated around t − t_merger ≈ 25 ms in the high-resolution simulation. Such efficient generation of the mean poloidal flux is not observed in the simulation with Δx_finest = 200 m for t − t_merger ≲ 80 ms, as plotted with the green dashed curve in the figure (but see below for details about the slight generation of the mean poloidal flux for 40 ≲ t − t_merger ≲ 60 ms). For t − t_merger ≳ 100 ms, the mean poloidal magnetic flux is generated even in the low-resolution simulation because the magnetorotational instability starts to be partially, not fully, resolved in the low-density region (see the middle and bottom panels of Extended Data Fig. 2). Therefore, we conclude that the magnetorotational instability is responsible for generating the mean poloidal magnetic flux.

Extended Data Fig. 5 plots meridional profiles of the mean radial magnetic field ${\bar{B}}_{r}$ at selected time slices for the simulation with Δx_finest = 12.5 (200) m in the left (right) column (see also the visualization: http://www2.yukawa.kyoto-u.ac.jp/~kenta.kiuchi/anime/SAKURA/movie_Mean_Poloidal_Flux.mp4). It is evident that the mean poloidal field is generated in the region where the magnetorotational instability is active (ρ ≤ 10^14.5 g cm⁻³), particularly, in the low-latitude region 0.24 ≲ θ ≲ 0.5π. Then, it propagates towards the high-latitude region and generates the jet. We also point out that at a radius of 9–10 km, the polar region has a weak magnetic field. This indicates that the poloidal field below this region does not contribute towards the launch of a jet and stays buried throughout the simulation.

The low-resolution simulation also shows some amplification of the mean poloidal field in the low-density region with ρ ≤ 10¹² g cm⁻³, where the magnetorotational instability is resolved. (see the visualizations for the mean poloidal flux and magnetorotational instability quality factor in the simulation with Δx_finest = 200 m: http://www2.yukawa.kyoto-u.ac.jp/~kenta.kiuchi/anime/SAKURA/movie_Mean_Poloidal_Flux_Low.mp4 and http://www2.yukawa.kyoto-u.ac.jp/~kenta.kiuchi/anime/SAKURA/movie_Mean_MRI_Qfac_Low.mp4). It shows qualitatively similar behaviour to the simulation with Δx_finest = 12.5 m, but the efficiency for generating the mean poloidal magnetic field is much lower, which leads to a weaker Poynting flux luminosity (see the ‘Detailed properties of the Poynting-flux-dominated outflow and magnetically driven post-merger ejecta’ section).

αΩ dynamo

To understand the dynamo process in the simulation, we used the mean field theory with an axisymmetric average. We consider an axisymmetric average because it corresponds to the symmetry of the differential rotation in the system. In the mean field dynamo theory, we assume that the velocity field U and the magnetic field B are, respectively, composed of the mean velocity field $\bar{{{{\bf{U}}}}}$ and velocity fluctuations u so that ${{{\bf{U}}}}=\bar{{{{\bf{U}}}}}+{{{\bf{u}}}}$ and of the mean magnetic field $\bar{{{{\bf{B}}}}}$ and the magnetic fluctuations b. We then average the induction equation, which gives equation (1), where $\bar{{{{\bf{{{{\mathcal{E}}}}}}}}}=\overline{{{{\bf{u}}}}\times {{{\bf{b}}}}}$ is the electromotive force due to the fluctuations. To close the system, the electromotive force is often expressed as a function of the mean magnetic field, and its spatial derivatives as described by equation (2). α_ij and β_ij are, respectively, tensors expressing the contributions of the mean magnetic field $\bar{{{{\bf{B}}}}}$ and its derivatives $\bar{{{{\bf{J}}}}}=\nabla \times \bar{{{{\bf{B}}}}}$ to the electromotive force. These tensors should not depend on $\bar{{{{\bf{B}}}}}$. First, we focus on how the mean poloidal field is generated and then look at the toroidal field.

To show how the mean poloidal field is generated by the alpha or beta effect, we assume that the mean velocity field is composed of the rotation speed $\bar{\bf{U}}=R\varOmega {\bf{e}}_{\phi }$ and project the averaged induction equation (1) in the radial direction in cylindrical coordinates, which gives:

$${\partial }_{t}{\bar{B}}_{R}=-{\partial }_{z}\left[{\left(\bar{{{{\bf{U}}}}}\times \bar{{{{\bf{B}}}}}\right)}_{\phi }+{\bar{{{{\bf{{{{\mathcal{E}}}}}}}}}}_{\phi }\right]=-{\partial }_{z}{\bar{{{{\mathcal{E}}}}}}_{\phi }=-{\partial }_{z}\left({\alpha }_{\phi j}{B}_{j}+{\beta }_{\phi j}{\left(\nabla \times \bar{{{{\bf{B}}}}}\right)}_{j}\right).$$

(12)

Similarly, by projecting the averaged induction equation in the vertical direction, the generation of the mean vertical field ${\bar{B}}_{z}$ is due to the azimuthal component of the electromotive force ${\bar{{{{\mathcal{E}}}}}}_{\phi }$. For the generation of the mean toroidal field, the Ω effect (that is the winding of the magnetic field by differential rotation) is also important and must be compared to the contributions of both the radial and vertical components of the electromotive force ${\bar{{{{\mathcal{E}}}}}}_{R}$ and ${\bar{{{{\mathcal{E}}}}}}_{z}$. To estimate which of the mean magnetic field ${\bar{B}}_{i}$ or the derivatives of the mean magnetic field ${\bar{J}}_{i}$ contributes the most to the electromotive force ${\bar{{{{\mathcal{E}}}}}}_{i}$, we compute the Pearson correlation coefficients between the quantities ${\bar{{{{\mathcal{E}}}}}}_{i}$ and ${Y}_{j}={\bar{B}}_{j}$ or ${\bar{J}}_{j}$ with the following formula:

$${\mathcal{C}}_\mathrm{P}({\bar{{{{\mathcal{E}}}}}}_{i},{Y}_{j})=\frac{{\int}_{t}({\bar{{{{\mathcal{E}}}}}}_{i}-{\langle {\bar{\mathcal{E}}}_{i}\rangle }_{t})({Y}_{j}-{\langle {Y}_{j}\rangle }_{t})\,\mathrm{d}t}{\sqrt{({\int}_{t}{({\bar{{{{\mathcal{E}}}}}}_{i}-{\langle {\bar{\mathcal{E}}}_{i}\rangle }_{t})}^{2}\,\mathrm{d}t)}\sqrt{({\int}_{t}{({Y}_{j}-{\langle {Y}_{j}\rangle }_{t})}^{2}\,\mathrm{d}t)}},$$

(13)

where 〈⋅〉_t represents a time average.

In the following sections, we present a complementary analysis of the other contributions to the generation of the mean magnetic field than the α_ϕϕ effect and Ω effect. We, therefore, show the other correlations between the electromotive force and the magnetic field and the estimated values of the tensor component. Several methods can be used for the values of the alpha tensor components. The simplest one is to estimate them from the correlations²³, but this method assumes that one component is dominant. To take into account the off-diagonal contributions, we compute the values of the alpha tensor coefficients in this study by using a singular value decomposition in a least-squares fit of the mean current data and mean field⁶⁶.

Generation of the mean poloidal field ${\bar{B}}_{R}$

As shown in equation (12), the generation of the axisymmetric poloidal field is due to the curl of the toroidal component of the electromotive force in the averaged induction equation. In this section, we show the correlations between the toroidal electromotive force ${\bar{{{{\mathcal{E}}}}}}_{\phi }$ and the magnetic field components ${\bar{B}}_{j}$ and the mean current ${\bar{J}}_{j}$ (top and middle panels of Extended Data Fig. 6). The low correlations with the mean current ${\bar{J}}_{j}$ shows that its contributions to ${\bar{{{{\mathcal{E}}}}}}_{\phi }$ (that is the β_ij tensor) can be neglected. In the same way, the vertical magnetic field contribution ${\bar{B}}_{z}$ can be neglected. Extended Data Fig. 6 also shows the anti-correlation of ${\bar{B}}_{\phi }$ and ${\bar{{{{\mathcal{E}}}}}}_{\phi }$ and also that the radial magnetic field ${\bar{B}}_{R}$ is correlated to ${\bar{{{{\mathcal{E}}}}}}_{\phi }$. This is because the radial field ${\bar{B}}_{R}$ is anti-correlated to the toroidal magnetic field ${\bar{B}}_{\phi }$ due to the Ω effect (Fig. 3). Since the mean toroidal field ${\bar{B}}_{\phi }$ is anti-correlated to the toroidal electromotive force ${\bar{{{{\mathcal{E}}}}}}_{\phi }$, the mean radial field ${\bar{B}}_{R}$ is correlated to it. To confirm that the contribution of ${\bar{B}}_{\phi }$ dominates, we first compute the non-diagonal components of the alpha tensor (bottom panel of Extended Data Fig. 6). We then compare the contribution of ${\bar{B}}_{\phi }$ and ${\bar{B}}_{R}$, that is, respectively, the time-averaged values of ${\alpha }_{\phi \phi }{\bar{B}}_{\phi }$ and ${\alpha }_{\phi R}{\bar{B}}_{R}$ in the first 10 km:

$$\frac{{\langle {\alpha }_{\phi \phi }{\bar{B}}_{\phi }\rangle }_{t}}{{\langle {\alpha }_{\phi R}{\bar{B}}_{R}\rangle }_{t}}=1.87.$$

(14)

Generation of the mean toroidal field

In the main text, we show that the Ω effect is important for the generation of the mean toroidal field. In this section, we check whether the contribution of the α effect through the poloidal components of the electromotive force ${\bar{{{{\mathcal{E}}}}}}_{R}$ and ${\bar{{{{\mathcal{E}}}}}}_{z}$ is substantial, in which case the dynamo is called an α²Ω dynamo instead of an αΩ dynamo. First, we confirm that the correlations between the electromotive force and the mean current ${\bar{J}}_{j}$ (right panels of Extended Data Fig. 7) are lower than with the mean magnetic field ${\bar{B}}_{j}$ (left panels of Extended Data Fig. 7). The contributions from the mean current can, therefore, be neglected. For the mean magnetic field, some components, for example, the mean toroidal field ${\bar{B}}_{\phi }$, are strongly correlated with the radial component of the electromotive force ${\bar{{{{\mathcal{E}}}}}}_{R}$. This indicates that the α effect might be important to the generation of the mean toroidal field. To compare the strength of these two effects, the α effect and the Ω effect, we computed the corresponding α tensor components α_R i and α_z i and estimated the ratio of the two dynamo numbers ${C}_{\alpha }=\max (| {\alpha }_{Ri}| ,| {\alpha }_{zi}| )R/\eta$ for i ∈ [R, ϕ, z] and C_Ω = ΩR²/η, where η is the resistivity, in the turbulent region averaged for one scale height, which gives at R = 30 km:

$$\frac{{C}_{\varOmega }}{{C}_{\alpha }}=\frac{\varOmega R}{\max (| {\alpha }_{Ri}| ,| {\alpha }_{zi}| )}\approx 30.8.$$

(15)

The contribution of the α effect to the generation of the mean toroidal field can, therefore, be reasonably neglected as the Ω effect dominates the generation of the mean toroidal field. The dynamo in our simulation can, therefore, be interpreted as an αΩ dynamo.

Detailed properties of the Poynting-flux-dominated outflow and magnetically driven post-merger ejecta

The link http://www2.yukawa.kyoto-u.ac.jp/~kenta.kiuchi/anime/SAKURA/DD2_135_135_Dynamo.mp4 is a visualization for the rest-mass density (top left), the magnetic field strength (top, second from left), the magnetization parameter (top, second from right), the unboundedness defined by the Bernoulli criterion (top right), the electron fraction (bottom left), the temperature (bottom, second from left), the entropy per baryon (bottom, second from right) and the geodesic criterion (bottom right) on a plane perpendicular to the orbital plane.

The top panel of Extended Data Fig. 8 plots the angular distribution of the luminosity of the Poynting flux, which is calculated as follows:

$${L}_{{{{\rm{Poynting}}}}}(\theta )=-{\int}_{r\approx 500{{{\rm{km}}}}}\alpha {\psi }^{6}{r}^{2}{({{T}^{r}}_{t})}_{\rm{(EM)}}\,\mathrm{d}\phi,$$

(16)

where α and ψ are the lapse function and the conformal factor, respectively. The high luminosity of ~2–8 × 10⁵² erg s⁻¹ angle⁻¹ is confined in a region with θ < 12^∘.

The middle panel plots the jet-opening-angle-corrected luminosity and luminosity of the Poynting flux (green and blue, respectively) as functions of the post-merger time. According to refs. ^67,68,69, the luminosity of the Poynting-flux-dominated outflow, which is driven by the efficient magnetic winding from the remnant of the binary neutron star merger, can be estimated as

$$\begin{array}{l}{L}_{{{{\rm{Poy}}}}} = \\ \int L_{\rm Poynting}(\theta) \sin\theta d\theta \approx10^{51}\,{\rm{erg}}\,{\rm{s}}^{-1}{\left(\frac{{\bar{B}}_{P}}{1{0}^{15}\,{\rm{G}}}\right)}^{2}{\left(\frac{R}{1{0}^{6}\,{\rm{cm}}}\right)}^{3}\left(\frac{\varOmega }{8,000\,{\rm{rad}}\,{\rm{s}}^{-1}}\right),\end{array}$$

(17)

where ${\bar{B}}_{P}$ is the azimuthally averaged poloidal magnetic field. In the current simulation, it is ~10¹⁵ G. Thus, the Poynting flux luminosity found in the simulation is consistent with this estimation. The jet-opening-angle-corrected luminosity is ~10⁵² erg s⁻¹. Thus, if we assume that the conversion efficiency to gamma-ray photons is ~10%, it is compatible with the observed luminosity of short, hard gamma-ray bursts⁴².

The bottom panel plots the ejecta as a function of the post-merger time. The solid curve denotes the result from using the Bernoulli criteria. The coloured region in the inset shows the error of the baryon mass conservation, which is below 10⁻⁷ M_⊙ throughout the simulation.

Convergence study and effect of the initial large-scale magnetic field

Because we chose an initial large-scale magnetic field strength that is much stronger than those observed in binary pulsars³⁶, we need to clarify whether such a strong large-scale field is responsible for launching the jet.

We begin by estimating the magnetic winding timescale originating from the initial magnetic field. Suppose we consider a binary neutron star merger with a highly magnetized end of 10¹¹ G, as observed for pulsars. In that case, the timescale for the magnetic winding to reach saturation is

$${\bar{B}}_{\phi } \approx10^{16.5}\,{{{\rm{G}}}}\left(\frac{{\bar{B}}_{R}}{5\times 1{0}^{11}\,{{{\rm{G}}}}}\right)\left(\frac{\varOmega }{8,000\,{\rm{rad}}\,{\rm{s}}^{-1}}\right)\left(\frac{t}{8\,{{{\rm{s}}}}}\right),$$

(18)

where we assume that the compression due to the merger amplifies the initial poloidal field by a factor of five, which should be proportional to ρ^2/3 because of conservation of the magnetic flux. We also assume a saturation field strength of 10^16.5 G, as suggested by the super-high-resolution simulation (Extended Data Fig. 1). Therefore, the magnetic winding originating from the initial magnetic field should be minor or irrelevant in reality.

However, the magnetic winding timescale is much shorter than those in reality in both simulations with Δx_finest = 12.5 or 200 m because of the assumed strong initial field. The low-resolution mesh is fine enough to resolve the compression and winding but can only partially capture the Kelvin–Helmholtz and the magnetorotational instability. Therefore, there would be no striking difference between the simulations for the jet launching mechanism if such a strong initial poloidal field and subsequent magnetic winding were responsible for it. The middle panel of Extended Data Fig. 8 shows that the Poynting-flux-dominated outflow is launched at t − t_merger ≈ 60 ms and reaches ~10^49–50 erg s⁻¹ in the low-resolution simulation with Δx_finest = 200 m. However, in the same plot, the super-high-resolution simulation with Δx_finest = 12.5 m shows that the Poynting-flux-dominated outflow is launched at t − t_merger ≈ 35 ms. The luminosity reaches ~10⁵¹ erg s⁻¹. The difference is striking. As discussed in Extended Data Fig. 5, the efficiency of the generation of the mean poloidal flux in the super-high-resolution simulation is much higher than that in the low-resolution simulation. This indicates that the efficient αΩ dynamo is responsible for launching the strong jet.

The bottom panel of Extended Data Fig. 8 also suggests that the Lorentz force-driven post-merger ejecta is launched in the run with Δx_finest = 200 m. However, the ejecta mass is about one order of magnitude smaller than in the run with Δx_finest = 12.5 m, showing that the enhanced activity of magnetohydrodynamics effects by the efficient dynamo action plays an important role in ejecting matter²¹.